************************************************************************ program stratmesh ************************************************************************ * * * Generates three-dimensional TOUGH2 mesh based on a rectangular * * surface grid with regular or irregular grid spacing. * * Multiple hydrostratigraphic units of varibale thickness can be * * specified. A material name is assigned to each unit. * * Each unit can be subdivided into several layers of elements. * * Boundary elements can be generated for each face of the model. * * * * Elevations of hydrostratigraphic units have to be provided at each * * point of the surface mesh. Future updates of STRATMESH may include * * some interpolation routine to calculate elevations at points where * * no direct measurements are available. * * * * Stefan Finsterle, Ron Falta, and Rex Hodges * * Clemson University, South Carolina, April 8, 1996 * * * * Formatted input has to be provided as follows: * * * * STRAT (Format: a5; block identifier for STRATMESH) * * * * NTYPE,NO,DEL,CONST (Format: a2,3x,i5,2e10.4) * * - NTYPE : set equal to NX, or NY for specifying grid increments * * in X or Y direction (surface mesh) * * - NO : number of grid increments (maximum: 99) * * - DEL : constant grid increment for NO grid blocks, if set to a * * non-zero value. * * - CONST : a constant added to the coordinates * * * * DEL(I),I=1,NO (Format: 8E10.4; optional, DEL=0.0 or blank only) * * - DEL(I) : a set of grid increments in the direction specified by * * NTYPE * * * * NS,NO (Format: a2,3x,i5) * * - NS : set equal to NS for specifying hydrostratigraphic units * * - NO : number of hydrostratigraphic units (maximum: 61) * * * * ROCK,NSUB (Format: a5,i5, NO records following keyword NS) * * - ROCK : Material name assigned to hydrostratigraphic unit * * - NSUB : Number of element layers per hydrostratigraphic unit * * * * ZS (Format: a2; block identifier) * * X,Y,ZS(I,J,K),K=1,NS+1 * * (Format: x(8e10.4,/) * * - X : X-coordinate, has to match coordinate of surface mesh * * - Y : Y-coordinate, has to match coordinate of surface mesh * * - ZS : Elevation of interfaces between hydrostratigraphic units * * There are NX*NY such records. * * For each grid block (I,J) provide NS+1 elevations, * * defining the boundaries between the NS hydrostrati- * * graphic units. Units can be specified either from top * * to bottom or bottom to top. * * * * BC (Format: a2; block identifier) * * IBOUND(I),I=1,6,NALLCON,TANMAX,DAMAX * * (Format: 6i1,i4,2e10.4) * * - IBOUND : Specifies treatment of boundary elements. * * Provide 6 integers corresponding to the 6 faces of the * * model, i.e.: * * I=1:LEFT; 2:RIGHT; 3:FRONT; 4:BACK; 5:TOP; 6:BOTTOM * * IBOUND(I)=0 --> no special treatment of elements * * IBOUND(I)=1 --> elements will be inactive (neg. volume) * * IBOUND(I)=2 --> large volume will be assigned to element * * IBOUND(I)=3 --> only one boundary element will be * * generated representing the entire face * * A fixed rock type ('LEFT ','RIGHT','FRONT','BACK ', * * 'TOP ', and 'BOTTO') will be assigned to the * * corresponding boundary elements. * * - NALLCON : By default, connections between boundary elements are * * omitted. Connections between boundary elements and the * * model domain have a very small nodal distance. * * If NALLCON is not euqal to 0, all connections (including * * those between boundary elements) are printed, and nodal * * distances are not reduced to small values. * * - TANMAX : A warning meassage will be printed if the tangent of a * * connection increases TANMAX (default: 0.2) * * - ZDIFFMAX: A warning message will be printed if the ratio bewteen * * adjacent element thicknesses is greater than ZDIFFMAX * * (default: 2.0) * * * * ENDFI (Format: a5; closes STRATMESH data block) * * * ************************************************************************ C C --- STRATMESH main program C character filename*80 write(*,7000) 7000 format(/,20x,'-----------------------', & /,20x,'| S T R A T M E S H |', & /,20x,'-----------------------',/, & /,' Generates three-dimensional TOUGH2/T2VOC mesh.', & /,' The surface mesh is rectangular', & ' with regular or irregular grid spacing.', & /,' Multiple hydrostratigraphic units of variable', & ' thicknesses are generated.', & /,' A material name is assigned to each unit.', & /,' Each unit may be subdivided into several layers', & ' of elements.',/, & /,' S. Finsterle, R. Falta, and R. Hodges', & /,' Clemson University', & /,' South Carolina', & /,' April 22, 1996',/) C C --- Read name of input file C write(*,'(a,$)') ' STRATMESH input file : ' read(*,'(a)') filename open(unit=10,file=filename,status='old',err=9000) lf=index(filename,'.') if (lf.eq.0) then call lenos(filename,lf) else lf=lf-1 endif C C --- Open output files C write(*,7001) filename(:lf),filename(:lf),filename(:lf) 7001 format(' STRATMESH output file : ',a,'.out',/, & ' TOUGH2 mesh file : ',a,'.mes',/, & ' TECPLOT file : ',a,'.tec') open(unit=11,file=filename(:lf)//'.out',status='unknown') open(unit=12,file=filename(:lf)//'.mes',status='unknown') open(unit=13,file=filename(:lf)//'.tec',status='unknown') C C --- Write header of STRATMESH output file C write(11,7002) filename(:lf),filename(:lf),filename(:lf) 7002 format(' STRATMESH',/, & ' ---------',/,/, & ' Input file : ',a,/, & ' Mesh file : ',a,'.mes',/, & ' Tecplot file : ',a,'.tec',/) C C --- Call mesh generator C call doitnow write(*,7003) 7003 format(/,' STRATMESH terminated.',/) stop C C --- Error message C 9000 continue call lenos(filename,lf) write(*,7004) filename(:lf) 7004 format(/,' ERROR. File ',a,' does not exist.') stop end C --- End of main program STRATMESH ************************************************************************ subroutine doitnow ************************************************************************ C C --- Set maximum problem size here C C maxx : maximum number of elements in X-direction ( < 99) C maxy : maximum number of elements in Y-direction ( < 99) C maxz : maximum number of elements in Z-direction ( < 61) C maxs : maximum number of hydrostratigraphic units ( < 61) C parameter (maxx=99,maxy=99,maxz=60,maxs=30) C C --- ----------------^^------^^------^^------^^------------------------ C C --- Set constants here C C big : Multiplication factor for volumes of boundary elements C small: Multiplication factor for connections to boundary elements C dtm : Default maximum tangent of connections C dzdm : Default maximum ratio of element thicknesses C parameter (big=1.e+20,small=1.e-08) parameter (dtm=2.e-01,dzdm =2.e+00) C C --- ---------------^^^^^^-------^^^^^^-------------------------------- C parameter (mkey=6) parameter (eps=1.e-10) character*5 ckey,word,rock,bcrock,rtype character*5 celem,celem1,celem00,celem11 character abc*1,nw*2,cinact*8,ctopbot*13 dimension dx(maxx),dy(maxy),xs(maxx),ys(maxy),elev(maxs) dimension ddx(2),ddy(2),ddz(2) dimension nsub(maxs),ibound(6),rock(maxs),zs(maxx,maxy,maxs) dimension izs(maxx,maxy) dimension nw(mkey),cinact(3),abc(61),bcrock(6),ctopbot(2) data nx,ny,nz,ns,iloop,jzs,nelbou,nelmod,ncon,nconbou/10*0/ data ibound/6*0/,izsok,nallcon,nline/3*0/ data itopbot/1/ data totvol,xtot,ytot/3*0./ data nw/'EN','NX','NY','NS','ZS','BC'/ data ctopbot/'top to bottom','bottom to top'/ data bcrock/'LEFT ','RIGHT','FRONT','BACK ','TOP ','BOTTO'/ data cinact/'inactive','very big','huge '/ data abc/'1','2','3','4','5','6','7','8','9', & 'A','B','C','D','E','F','G','H','I','J','K','L','M', & 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z', & 'a','b','c','d','e','f','g','h','i','j','k','l','m', & 'n','o','p','q','r','s','t','u','v','w','x','y','z'/ C C --- Check maximum problem size C if (maxx.gt.99) call error(4,nline, & 'Max. number of elements in X-direction is 99') if (maxy.gt.99) call error(4,nline, & 'Max. number of elements in Y-direction is 99') if (maxz.gt.61) call error(4,nline, & 'Max. number of elements in Z-direction is 61') if (maxs.gt.61) call error(4,nline, & 'Max. number of hydrostratigraphic units is 61') C C --- Find STRATMESH block name C 1000 read(10,'(a)',end=9000) word nline=nline+1 if (word.ne.'STRAT') goto 1000 C C --- Read keywords C 2000 continue read(10,7000,iostat=ios,end=9001) ckey,no,del,const nline=nline+1 if (ios.ne.0) call error(1,nline,'CKEY,NO,DEL,CONST') 7000 format(a5,i5,2e10.4) do 1001 i=1,mkey if (ckey(:2).eq.nw(i)) & goto (2001,2002,2003,2004,2005,2006),i 1001 continue if (ckey.ne.' ') write(11,7001) ckey 7001 format(' WARNING. Have read unknown block label ''',a,'''.') goto 2000 C C --- Read X-gridding C 2002 continue nx1=nx+1 nx=nx+no if (nx.gt.maxx) call error(3,nline,'maxx') xconst=const if (del.ne.0.) then do 1002 i=nx1,nx dx(i)=del 1002 continue else read(10,7002,iostat=ios) (dx(i),i=nx1,nx) nline=nline+(nx-nx1)/8+1 if (ios.ne.0) call error(1,nline,'DX(I)') 7002 format(8e10.4) endif goto 2000 C C --- Read Y-gridding C 2003 continue ny1=ny+1 ny=ny+no if (ny.gt.maxy) call error(3,nline,'maxy') yconst=const if (del.ne.0.) then do 1003 i=ny1,ny dy(i)=del 1003 continue else read(10,7002,iostat=ios) (dy(i),i=ny1,ny) nline=nline+(ny-ny1)/8+1 if (ios.ne.0) call error(1,nline,'DY(I)') endif goto 2000 C C --- Read hydrostratigraphic units C 2004 continue ns1=ns+1 ns=ns+no if (ns.ge.maxs) call error(3,nline,'maxs') do 1004 i=ns1,ns read(10,7003,iostat=ios) rock(i),nsub(i) nline=nline+1 7003 format(a5,i5) if (ios.ne.0) call error(1,nline,'ROCK,NSUB') nz=nz+nsub(i) if (nz.gt.maxz) call error(3,nline,'maxz') 1004 continue goto 2000 C C --- Read elevations of interfaces between hydrostratigraphic units C 2005 continue if (nx.eq.0.or.ny.eq.0) & call error(4,nline,'Define NX and NY before NS') C C --- Calculate coordinates of surface mesh first C xs(1)=xconst+dx(1)/2. do 1005 i=2,nx xs(i)=xs(i-1)+(dx(i)+dx(i-1))/2. 1005 continue ys(1)=yconst+dy(1)/2. do 1006 i=2,ny ys(i)=ys(i-1)+(dy(i)+dy(i-1))/2. 1006 continue iel=0 jzs=jzs+1 C C --- Read nx*ny elevation sets C kkk=0 2007 continue read(10,7004,iostat=ios) x,y,(elev(kk),kk=1,ns+1) nline=nline+(ns+1)/8+1 7004 format(8e10.4) if (ios.ne.0) then backspace(10) nline=nline-1 goto 2009 endif kkk=kkk+1 do 1007 i=1,nx if (abs(x-xs(i)).lt.eps) goto 2008 1007 continue write (11,7005) 'X',kkk,x,y 7005 format(' ERROR. ',a1,'-coordinate of elevation set no.',i4, & ' is incorrect.',/, & ' Coordinates specified:',E15.6,' /',E15.6,/) goto 2007 2008 continue do 1008 j=1,ny if (abs(y-ys(j)).lt.eps) goto 2010 1008 continue write (11,7005) 'Y',kkk,x,y goto 2007 2010 continue iel=iel+1 izs(i,j)=izs(i,j)+1 do 1009 l=1,ns+1 zs(i,j,l)=elev(l) if (l.gt.2.and.kkk.eq.1) then if (zs(i,j,l-1)-zs(i,j,l).lt.0.) itopbot=2 if (((zs(i,j,l-1)-zs(i,j,l)).lt.0.and.itopbot.eq.1).or. & ((zs(i,j,l-1)-zs(i,j,l)).gt.0.and.itopbot.eq.2)) then write(11,7006) iel,i,j,xs(i),ys(j),l-1, & zs(i,j,l-1),l,zs(i,j,l) 7006 format(' ERROR. Inconsistency detected in elevations.',/,/, & ' Cell number (I,J) :',i6,5x,'(',i10,' ,',i10,')',/, & ' Coordinates (X,Y) :',11x,'(',e10.4,' ,',e10.4,')',/,/, & ' Sequence of elevations: unit elevation',/, & 2(24x,i6,e15.6,/),/, & ' Consistently use either increasing', & ' or decreasing elevations.') write(*,'(/,a)') ' ERROR in STRATMESH input file' stop endif endif 1009 continue goto 2007 2009 continue do 1010 j=1,ny do 1011 i=1,nx if (izs(i,j).ne.1) then izsok=izsok+1 write(11,7007) izs(i,j),i,j,xs(i),ys(j) 7007 format(/,' ERROR.',i3, & ' elevation set(s) specified at grid point (', & i2,'/',i2,').',/,' Coordinates:',E15.6,' /',E15.6) endif 1011 continue 1010 continue goto 2000 2006 continue C C --- Read which faces of the model should be prepared for assigning C boundary conditions C read(10,'(6i1,i4,2e10.4)',iostat=ios) & (ibound(i),i=1,6),nallcon,tanmax,zdiffmax nline=nline+1 if (ios.ne.0) call error(1,nline, & 'IBOUND(I),I=1..6,NALLCON,TANMAX,ZDIFFMAX') if (tanmax.le.0.) tanmax=dtm if (zdiffmax.lt.1.) zdiffmax=dzdm goto 2000 2001 continue C C --- Go through loop 3 times. First time to generate element information C for model domain, second time for inactive boundary elements, and C third time for connections C iloop=iloop+1 if (iloop.eq.1) then C C --- Print header information C write(13,7008) nx,ny,nz 7008 format('VARIABLES = X',19x,'Y "Z or ROCKTYPE"',/, & 'ZONE t="3D-MESH" i =',i3,' j =',i3,' k =',i3) write(12,7009) nx,ny,nz 2013 continue 7009 format('ELEME NX=',i4,' NY=',i4,' NZ=',i4) write(11,7010) nx,ny,nz,nx*ny*nz,ns,ctopbot(itopbot) 7010 format(/,' Number of elements in x-direction :',i6,/, & ' Number of elements in y-direction :',i6,/, & ' Number of elements in z-direction :',i6,/, & ' Total number of elements :',i6,/, & ' Number of hydrostratigraphic units:',i6,2x, & '(from ',a13,')',/,/, & ' Unit Material Subdivisions') do 1012 i=1,ns write(11,7011) i,rock(i),nsub(i) 7011 format(i5,9x,a5,i14) 1012 continue write(11,'(/)') kloop=2 C C --- Check boundary element information C do 1013 i=1,6 if (ibound(i).gt.3) & call error(4,nline,'Unknown option for IBOUND') if (ibound(i).ne.0) then kloop=1 write(11,7012) bcrock(i),cinact(ibound(i)) 7012 format(' Elements representing ',a5,' face of model are ',a8) endif 1013 continue C C --- Check whether layers are defined from top to bottom or C bottom to top. Switch definition of top and bottom face accordingly. C if (itopbot.eq.2) then idum=ibound(5) ibound(5)=ibound(6) ibound(6)=idum bcrock(5)='BOTTO' bcrock(6)='TOP ' endif write(11,7013) 'X',(dx(i),i=1,nx) write(11,7013) 'Y',(dy(i),i=1,ny) 7013 format(/,/,' Grid spacing in ',a1,'-direction',/, & ' ---------------------------',/,/, & 100(5e13.6,/)) if (nz.eq.0) then C C --- Calculate coordinates of surface mesh C xs(1)=xconst+dx(1)/2. do 1105 i=2,nx xs(i)=xs(i-1)+(dx(i)+dx(i-1))/2. 1105 continue ys(1)=yconst+dy(1)/2. do 1106 i=2,ny ys(i)=ys(i-1)+(dy(i)+dy(i-1))/2. 1106 continue endif write(11,7014) 'X',(xs(i),i=1,nx) write(11,7014) 'Y',(ys(i),i=1,ny) 7014 format(/,/,1x,a1,'-coordinates',/, & ' -------------',/,/,100(5e13.6,/)) write(11,7115) nx,ny,nx*ny 7115 format(/,/,' Surface mesh has',i3,' *',i3,' =',i5, & ' pairs of coordinates:',/,/, & 4x,'#',4x,'I',4x,'J',19x,'X',19x,'Y') itot=0 do 1107 i=1,nx do 1108 j=1,ny itot=itot+1 write(11,'(3i5,2e20.6)') itot,i,j,xs(i),ys(j) 1108 continue 1107 continue if (izsok.ne.0) call error(1,nline,'X,Y,ELEV(K),K=1,NS+1)') if (nz.eq.0) then call error(4,nline,'Surface mesh calculated only') stop endif write(11,7015) 7015 format(/,/,' # ELEME ROCK X Y', & ' Z DX DY DZ') else if (iloop.eq.3) then write(12,'(/,a5)') 'CONNE' endif nzbot=0 iel=0 C C --- Main loop over hydrostratigraphic units starts here C do 1014 l=1,ns nztop=nzbot+1 nzbot=nztop+nsub(l)-1 kk=0 C C --- Loop over element layers starts here C do 1015 k=nztop,nzbot kk=kk+1 C C --- Loop over Y-gridding starts here C do 1016 j=1,ny y=ys(j) ddy(1)=dy(j)/2. if (j.lt.ny) ddy(2)=dy(j+1)/2. C C --- Innermost loop over X-gridding starts here C do 1017 i=1,nx celem(1:1)=abc(k) write(celem(2:3),'(i2)') j write(celem(4:5),'(i2)') i x=xs(i) ddx(1)=dx(i)/2. if (i.lt.nx) ddx(2)=dx(i+1)/2. dzk=(zs(i,j,l)-zs(i,j,l+1))/nsub(l) if (kk.lt.nsub(l)) then dzk1=dzk else if (l.lt.ns) then dzk1=(zs(i,j,l+1)-zs(i,j,l+2))/nsub(l+1) endif z=zs(i,j,l)-dzk/2.-dzk*(kk-1) ddz(1)=abs(dzk)/2. if (k.lt.nz) ddz(2)=abs(dzk1)/2. iel=iel+1 if (iloop.eq.1) then write(11,7016) iel,celem,rock(l), & x,y,z,dx(i),dy(j),abs(dzk) write(13,'(3e20.10)') x,y,z endif 7016 format(i5,' <',a5,'> ',a5,6f10.3) if (iloop.le.2) then C C --- Generate block ELEME C aij=0. if (k.eq.1) aij=aij+dx(i)*dy(j) if (k.eq.nz) aij=aij+dx(i)*dy(j) volume=dx(i)*dy(j)*abs(dzk) iprint=0 if (kloop.eq.1) then C C --- Check for boundary elements C do 1018 m=1,6 if (ibound(m).gt.0) then if((m.eq.1.and.i.eq.1).or. & (m.eq.2.and.i.eq.nx).or. & (m.eq.3.and.j.eq.1).or. & (m.eq.4.and.j.eq.ny).or. & (m.eq.5.and.k.eq.1).or. & (m.eq.6.and.k.eq.nz)) then rtype=bcrock(m) if (ibound(m).eq.1) then volume=-volume if (iloop.eq.2) iprint=1 else if (ibound(m).eq.2) then volume=volume*big if (iloop.eq.1) iprint=1 else if (ibound(m).eq.3) then ibound(m)=99 volume=big celem=bcrock(m) celem(4:5)=' 0' if (iloop.eq.1) iprint=1 else if (ibound(m).eq.99) then iprint=0 endif C C --- Boundary element C if (iprint.eq.1) nelbou=nelbou+1 goto 2012 endif endif 1018 continue endif C C --- Regular element C if (iloop.eq.1) then iprint=1 nelmod=nelmod+1 totvol=totvol+volume rtype=rock(l) endif 2012 continue C C --- Write element information C if (iprint.eq.1) & write(12,7017) celem,rtype,volume,aij,x,y,z 7017 format(a5,10x,a5,2e10.4,10x,3e10.4) else C C --- Generate block CONNE C if (kloop.eq.1) then C C --- Check for connections to boundary elements C do 1019 m=1,6 if (ibound(m).eq.0) goto 1019 celem00='!' celem11='!' if (m.eq.1.and.i.eq.1) then ddf1=small ddf2=1. iprintx=1 iprinty=0 iprintz=0 if (ibound(m).eq.99) celem00='LEF 0' else if (m.eq.2.and.i.eq.nx-1) then ddf1=1. ddf2=small iprintx=1 iprinty=0 iprintz=0 if (ibound(m).eq.99) celem11='RIG 0' else if (m.eq.3.and.j.eq.1) then ddf1=small ddf2=1. iprintx=0 iprinty=1 iprintz=0 if (ibound(m).eq.99) celem00='FRO 0' else if (m.eq.4.and.j.eq.ny-1) then ddf1=1. ddf2=small iprintx=0 iprinty=1 iprintz=0 if (ibound(m).eq.99) celem11='BAC 0' else if (m.eq.5.and.k.eq.1) then ddf1=small ddf2=1. iprintx=0 iprinty=0 iprintz=1 if (ibound(m).eq.99) then celem00=bcrock(m) celem00(4:5)=' 0' endif else if (m.eq.6.and.k.eq.nz-1) then ddf1=1. ddf2=small iprintx=0 iprinty=0 iprintz=1 if (ibound(m).eq.99) then celem00=bcrock(m) celem00(4:5)=' 0' endif else goto 1019 endif C C --- Connection to boundary element C nconbou=nconbou+1 goto 2011 1019 continue endif C C --- Regular connection C iprintx=1 iprinty=1 iprintz=1 ddf1=1. ddf2=1. celem00='!' celem11='!' 2011 continue C C --- Force printout of all connections between boundary elements C if (nallcon.ne.0) then iprintx=1 iprinty=1 iprintz=1 ddf1=1. ddf2=1. endif ncon=ncon+iprintx+iprinty+iprintz C C --- Print connections in X-direction C if (i.lt.nx.and.iprintx.eq.1) then dzki1=(zs(i+1,j,l)-zs(i+1,j,l+1))/nsub(l) zi1=zs(i+1,j,l)-dzki1/2.-dzki1*(kk-1) ayz=dy(j)*(abs(dzk)+abs(dzki1))/2. beta=(z-zi1)/sqrt((ddx(1)+ddx(2))**2+(zi1-z)**2) celem1=celem write(celem1(4:5),'(i2)') i+1 if (celem00(1:1).ne.'!') celem=celem00 if (celem11(1:1).ne.'!') celem1=celem11 write(12,7018) celem,celem1,1, & ddx(1)*ddf1,ddx(2)*ddf2,ayz,beta C C --- Check tangent and thickness ratio in X-direction C if (abs(beta).gt.tanmax) write(11,7019) & celem,celem1,celem,z,celem1, & zi1,ddx(1)+ddx(2),(z-zi1)/(ddx(1)+ddx(2)) arat=abs(dzk/dzki1) if (arat.lt.1.) arat=1./arat if (arat.gt.zdiffmax) write(11,7020) & celem,dzk,celem1,dzki1,arat endif C C --- Print connections in Y-direction C if (j.lt.ny.and.iprinty.eq.1) then dzkj1=(zs(i,j+1,l)-zs(i,j+1,l+1))/nsub(l) zj1=zs(i,j+1,l)-dzkj1/2.-dzkj1*(kk-1) axz=dx(i)*(abs(dzk)+abs(dzkj1))/2. beta=(z-zj1)/sqrt((ddy(1)+ddy(2))**2+(zj1-z)**2) celem1=celem write(celem1(2:3),'(i2)') j+1 if (celem00(1:1).ne.'!') celem=celem00 if (celem11(1:1).ne.'!') celem1=celem11 write(12,7018) celem,celem1,2, & ddy(1)*ddf1,ddy(2)*ddf2,axz,beta C C --- Check tangent and thickness ratio in Y-direction C if (abs(beta).gt.tanmax) write(11,7019) & celem,celem1,celem,z,celem1, & zj1,ddy(1)+ddy(2),(z-zj1)/(ddy(1)+ddy(2)) arat=abs(dzk/dzkj1) if (arat.lt.1.) arat=1./arat if (arat.gt.zdiffmax) write(11,7020) & celem,dzk,celem1,dzkj1,arat endif C C --- Print connections in Z-direction C if (k.lt.nz.and.iprintz.eq.1) then axy=dx(i)*dy(j) beta=dzk/abs(dzk) celem1=celem write(celem1(1:1),'(a1)') abc(k+1) if (celem00(1:1).ne.'!') celem=celem00 if (celem11(1:1).ne.'!') celem1=celem11 write(12,7018) celem,celem1,3, & ddz(1)*ddf1,ddz(2)*ddf2,axy,beta endif 7018 format(2a5,19x,i1,3e10.4,f10.7) 7019 format(/,' WARNING. Steep slope of connection between element ', & a5,' and ',a5,' !',/, & ' Elevation of element <',a5,'> :',e15.5,/, & ' Elevation of element <',a5,'> :',e15.5,/, & ' Horzontal distance between elements :',e15.5,/, & ' Tangent :',e15.5) 7020 format(/,' WARNING. Large contrast in element thickness!',/, & ' Thickness of element <',a5,'> :',e15.5,/, & ' Thickness of element <',a5,'> :',e15.5,/, & ' Ratio :',e15.5) endif 1017 continue 1016 continue 1015 continue 1014 continue if (iloop.eq.1) iloop=kloop if (iloop.lt.3) goto 2001 C C --- End of three loops; produce mesh statistics. C do 1020 i=1,nx xtot=xtot+dx(i) 1020 continue do 1021 j=1,ny ytot=ytot+dy(j) 1021 continue write(11,7021) nelmod+nelbou,nelmod,nelbou,ncon,ncon-nconbou, & nconbou,xtot*ytot,totvol 7021 format(/, &' Total number of elements :',i15,/, &' Total number of elements within model domain :',i15,/, &' Total number of boundary elements :',i15,/, &' Total number of connections :',i15,/, &' Total number of connections within model domain :',i15,/, &' Total number of connections to boundary elements :',i15,/, &' Surface area of model :',e15.8,/, &' Total volume of model :',e15.8) C C --- Generate more output for TECPLOT C C --- Vertical slices in X-direction C do 1022 i=1,nx write(13,7022) 'X',xs(i),ny,nz 7022 format('ZONE t = "Vertical slice at ',a1,' =',e12.6,'"', & ' i =',i3,' j =',i3) nzbot=0 do 1023 l=1,ns nztop=nzbot+1 nzbot=nztop+nsub(l)-1 kk=0 do 1024 k=nztop,nzbot kk=kk+1 do 1025 j=1,ny dzk=(zs(i,j,l)-zs(i,j,l+1))/nsub(l) z=zs(i,j,l)-dzk/2.-dzk*(kk-1) write(13,'(2e20.10,i20)') ys(j),z,l 1025 continue 1024 continue 1023 continue 1022 continue C C --- Vertical slices in Y-direction C do 1026 j=1,ny write(13,7022) 'Y',ys(j),nx,nz nzbot=0 do 1027 l=1,ns nztop=nzbot+1 nzbot=nztop+nsub(l)-1 kk=0 do 1028 k=nztop,nzbot kk=kk+1 do 1029 i=1,nx dzk=(zs(i,j,l)-zs(i,j,l+1))/nsub(l) z=zs(i,j,l)-dzk/2.-dzk*(kk-1) write(13,'(2e20.10,i20)') xs(i),z,l 1029 continue 1028 continue 1027 continue 1026 continue return C C --- Error messages C 9000 continue call error(2,nline,'STRAT') 9001 continue call error(2,nline,'ENDFI') end C --- End of doitnow ************************************************************************ subroutine error(imes,nline,message) ************************************************************************ C C --- Write error messages C character*(*) message call lenos(message,lm) if (imes.eq.1) then write(11,7000) nline,message(:lm) 7000 format(/,' ERROR near line:',i4,' while reading ',a,'!') else if (imes.eq.2) then write(11,7001) message(:lm) 7001 format(/,' ERROR. Keyword ',a,' not found.') else if (imes.eq.3) then write(11,7003) message(:lm) 7003 format(/,' ERROR. Increase parameter ',a,' and recompile.') else if (imes.eq.4) then write(11,7004) message(:lm) 7004 format(/,' ERROR. ',a,'.') endif write(*,9000) 9000 format(/,' ERROR in STRATMESH input file.', & ' Check STRATMESH output file for details.',/) stop end C --- End of error ************************************************************************ subroutine lenos(string,lstr) ************************************************************************ C C --- Determine lenth of string without space C character*(*) string lstr=0 ilenstr=len(string) if (ilenstr.eq.0) return do 1000 i=ilenstr,1,-1 if (string(i:i).ne.' ') then lstr=i return endif 1000 continue end C --- End of lenos