c Method Of Splitting Tsunami (MOST) model c (c) Vasily Titov 1998 c Finite difference algorithm computes tsunami wave (long wave) c evolution with splitting method (method of fractional steps) c using nested grids with variable space resolution. c*********************************************************************** c SOFTWARE LICENSE AND DISCLAIMER c c Scope of License. Subject to all the terms and conditions of this c license, DOC grants USER the royalty-free, nonexclusive, c nontransferable, and worldwide rights to reproduce, modify, and c distribute ComMIT and MOST Model Software, herein referred to as the c Product. c c Conditions and Limitations of Use Warranties. Neither the U.S. c Government, nor any agency or employee thereof, makes any warranties, c expressed or implied, with respect to the Product provided under this c License, including but not limited to the implied warranties or c merchantability and fitness for any particular purpose. c c Liability. In no event shall the U.S. Government, nor any agency or c employee thereof, be liable for any direct, indirect, or consequential c damages flowing from the use of the Product provided under this c License. c c Non-Assignment. Neither this License nor any rights granted hereunder c are transferable or assignable without the explicit prior written c consent of DOC. c c Names and Logos. USER shall not substitute its name or logo for the c name or logo of DOC, or any of its agencies, in identification of the c Product. c c Export of technology. USER shall comply with all U.S. laws and c regulations restricting the export of the Product to other countries. c c Governing Law. This License shall be governed by the laws of United c States as interpreted and applied by the Federal courts in the c District of Columbia. c c Term of License. This License shall remain in effect as long as USER c uses the Product in accordance with Paragraphs 1 and 2. c*********************************************************************** c Array descriptions: c u* - x-velocities (along longitude) c v* - y-velocities (along latitude) c q* - total wave height (depth + amplitude) c d* - depth c i* - initial bottom deformation c x* - values on vertical boundaries c y* - values on horizontal boundaries c in* -coordinates of corners of no.$ nested grid in * array c in*(1,$)=x1,in*(2,$)=x2,in*(3,$)=y1,in*(4,$)=y2 c Constants: c m* , n* -> x-, y-dimensions of arrays c Additional notes c Last modiified : Arun Chawla 3/20/2006, DCB 3/20/2006 c########################################################################### c This version of the MOST model uses 3 nested grids for computation # c and initial conditions are provided from a source file. Except for # c the parameter file and the three grid files all I/O is done ussing # c the netcdf file format. # c # c The parameter file is named most3_facts_nc.in and its format is # c >> MIN WAVE -> minimum amplitude to start computation # c >> MIN DEPTH -> minimum depth for offshore computation # c >> DRY DEPTH -> dry land depth for inundation # c >> FRIC COEFF -> manning's coefficient for bottom friction # c >> RUNUP COEFF -> switch to determine if there is runup in # c grids A and B (0 = no runup 1 = runup) # c >> BLOWUP -> max value for eta (m), if exceeded run terminates# c >> TSTEP -> input time step in seconds # c >> NTEPS -> input number of steps # c >> NA -> Multiple of input time step to Compute "A" # c arrays # c >> NB -> Multiple of input time step to Compute "B" # c arrays # c >> NOUT -> Multiple of input time step for output # c (has to be a multiple of NA and NB) # c >> NSTART -> Starting time step # c >> NSKIP -> Nodes to skip in saving to output # c (NSKIP = 1, if saving all the nodes to output) # c >> GRID A -> File name for coarsest grid # c >> GRID B -> File name for the next nested grid # c >> GRID C -> File name for the finest nested grid # c >> DIR SOURCE -> Directory path where the source files are # c >> DIR OUT -> Directory path for output files # c # c # c The command line arguments for running the model are # c $$> most3_facts_nc [dir] # c where # c -> common prefix assigned to all output files # c -> common prefix for all source files # c -> directory where the 3 grids and par files # c are kept (they have to be together) # c [dir] -> optional argument if set and not "false" # c then input parameter file read from # c if not set or false read from "site"_run2d/ # c # c The code can be compiled using the attached make file # c no nest/narr/in one way nesting, no talkback, friction, froude cap # c########################################################################### program MOST include 'netcdf.inc' integer*4 nx_nc,ny_nc,icnt #ifdef SMALLMEM parameter (nx0 = 400,ny0 = 400,n_t0=4500) parameter (nx0g = 3000, ny0g = 3000) parameter (nxa1 = 2000, nya1 = 2000) parameter (nxb1 = 2000, nyb1 = 2000) parameter (nxc1 = 2000, nyc1 = 2000) parameter (ndum = 2000) parameter (ntsmax = 3000) #else parameter (nx0 = 1000,ny0 = 1000,n_t0=4500) parameter (nx0g = 3000, ny0g = 3000) parameter (nxa1 = 3000, nya1 = 3000) parameter (nxb1 = 3000, nyb1 = 3000) parameter (nxc1 = 3000, nyc1 = 3000) parameter (ndum = 3000) parameter (ntsmax = 3000) #endif c parameters for initial condition interpolation parameter (lw = 8*(nxa1+nya1), liw = 2*(nxa1+nya1)) c..... 0-level arrays real*8 ua0(nx0,ny0),va0(nx0,ny0),qa0(nx0,ny0),da0(nx0,ny0) real*8 ua00(nx0,ny0),va00(nx0,ny0),qa00(nx0,ny0) real*8 ua01(nx0,ny0),va01(nx0,ny0),qa01(nx0,ny0) real*8 h1a0(nx0), h2a0(ny0), t0(n_t0), w(lw) real*8 x1ra0(ny0,3,1),x2la0(ny0,3,1) & ,y1ra0(nx0,3,1),y2la0(nx0,3,1) real*8 maxa0, taper integer itcalc integer ina0(4,3),outa0(4,3),crnra0(4,3),iw(liw),intpol(2) c..... deformation arrays integer*4 dx0,dy0 integer i,j real*8 h1a0g(nx0g),h2a0g(ny0g) real*8 ddeform(nx0g,ny0g) real*8 defthr,defmax logical deformislocal c..... 1-level arrays real*8 deformA(nxa1,nya1) real*8 ua(nxa1,nya1),va(nxa1,nya1),qa(nxa1,nya1),da(nxa1,nya1) c & ,ia(nxa1,nya1) real*8 h1a(nxa1), h2a(nya1) real*8 x1ra(nya1,3,1),x2la(nya1,3,1) & ,y1ra(nxa1,3,1),y2la(nxa1,3,1) real*8 maxa,mina integer mxai,mxaj,mnai,mnaj real*4 maxha(nxa1,nya1), minha(nxa1,nya1), maxspa(nxa1,nya1) integer ina(4,3),outa(4,3),crnra(4,3) c..... 2-level arrays real*8 deformB(nxb1,nyb1) real*8 ub1(nxb1,nyb1),vb1(nxb1,nyb1),qb1(nxb1,nyb1),db1(nxb1,nyb1) c & ,ib1(nxb1,nyb1) real*8 h1b1(nxb1), h2b1(nyb1),PI,RE,GA,dt,t,grnd,dry,fi,fr real*8 x1rb1(nyb1,3,1),x1lb1(nya1,3,1),x2rb1(nya1,3,1) & ,x2lb1(nyb1,3,1),y1rb1(nxb1,3,1),y1lb1(nxa1,3,1) & ,y2rb1(nxa1,3,1),y2lb1(nxb1,3,1) real*8 maxb,minb integer mxbi,mxbj,mnbi,mnbj real*4 maxhb(nxb1,nyb1), minhb(nxb1,nyb1), maxspb(nxb1,nyb1) integer inb1(4,3),outb1(4,3),crnrb1(4,3) c..... 3-level arrays real*8 deformC(nxc1,nyc1) real*8 uc1(nxc1,nyc1),vc1(nxc1,nyc1),qc1(nxc1,nyc1),dc1(nxc1,nyc1) c & ,ic1(nxc1,nyc1) real*8 h1c1(nxc1), h2c1(nyc1),minwave logical nowave, first real*8 x1rc1(nyc1,3,1),x1lc1(nyb1,3,1),x2rc1(nyb1,3,1) & ,x2lc1(nyc1,3,1),y1rc1(nxc1,3,1),y1lc1(nxb1,3,1) & ,y2rc1(nxb1,3,1),y2lc1(nxc1,3,1) real*8 maxc,minc integer mxci,mxcj,mnci,mncj real*4 maxhc(nxc1,nyc1), minhc(nxc1,nyc1), maxspc(nxc1,nyc1) C uncomment to get maximum values at each timestep c real*8 maxia(nxa1,nya1),maxib(nxb1,nyb1),maxic(nxc1,nyc1) c integer*4 dax c..... Dummy arrays for timestep computations c..... (dimension should be larger than all dims) real*8 uw(ndum),vw(ndum),qw(ndum),dw(ndum),u1(ndum),q1(ndum) & ,v1(ndum),h(ndum) c..... netCDF output arrays character*120 fname,fname1,fname2,fname3,dumb,prefix,site character*120 indir,output,tmp,durlin,durlout,defprefix character ga_hst*250,ga_ttl*250,ds*5 integer*4 ndims,nvar,ncunl_id,ngatts,inlen(12),imxchr,ga_len,drstr c.....error status return integer*4 iret c.....netCDF id for output to all 3 grids integer*4 axisids(3) integer*4 ncid1C,ncid2C,ncid3C integer*4 ncid1B,ncid2B,ncid3B integer*4 ncid1A,ncid2A,ncid3A integer*4 ncid1in,ncid2in,ncid3in integer*4 TIME_idC, HA_idC, UA_idC, VA_idC integer*4 TIME_idB, HA_idB, UA_idB, VA_idB integer*4 TIME_idA, HA_idA, UA_idA, VA_idA integer*4 MAXH_idC, MAXU_idC, MAXV_idC integer*4 MAXH_idB, MAXU_idB, MAXV_idB integer*4 MAXH_idA, MAXU_idA, MAXV_idA integer*4 TIME_idin, HA_idin, UA_idin, VA_idin integer*4 DEF_id, GLON_id, GLAT_id c..... netCDF data variables real*8 lat(nyc1),lon(nxc1),ct(1) real*4 cUA(nxc1,nyc1,1),cVA(nxc1,nyc1,1),cHA(nxc1,nyc1,1) real*4 bUA(nxb1,nyb1,1),bVA(nxb1,nyb1,1),bHA(nxb1,nyb1,1) real*4 aUA(nxa1,nya1,1),aVA(nxa1,nya1,1),aHA(nxa1,nya1,1) real*4 rUA(nx0,ny0,1),rVA(nx0,ny0,1),rHA(nx0,ny0,1) c..... program runtime information character*10 date(3) integer values(8), time1, time2, time real eee, etime, eeet(2) c.....error checking variables integer nargs, CheckVar real maxerror, capmax c.....runup switch integer irunup c..... timeseries locations for output integer*4 numts, tsi(ntsmax),tsj(ntsmax),siftid,tstimeid,tsid integer*4 mxaid,mnaid,spaid,mxbid,mnbid,spbid,mxcid,mncid,spcid integer*4 floodid character*80 fnsift real*8 tslat(ntsmax), tslon(ntsmax) integer*4 tsgrid(30),suppressA,suppressB,suppressC,suppressS data suppressA,suppressB,suppressC,suppressS/1,1,1,0/ c..... switch when finding image-ordered bathymetry files logical latimgA,latimgB,latimgC c..... num points inundated in C-grid integer*4 numinun integer*4 inunc(nxc1,nyc1) logical inunnotset integer inunthreshold common /run/ t,nstep,npics,grnd,dry,fi,fr common /const/ PI,RE,GA common /nc/ ns1,ne1,ns2,ne2,istep,jstep,scl_m2cm common /arr/ qa,ua,va,qb1,ub1,vb1,qc1,uc1,vc1 data mma,nna/nxa1,nya1/,mmb1,nnb1/nxb1,nyb1/,mmc1,nnc1/nxc1,nyc1/ data mma0,nna0/nx0,ny0/ data PI/3.14159265/,RE/1.569612e-07/,GA/9.81D0/,scl_m2cm/100./ C uncomment to get maximum values at each timestep in surfer format c data maxia/4000000*0.0/ c data maxib/4000000*0.0/ c data maxic/4000000*0.0/ c data dax/1/ c..... interpolation switches for initial conditions data intpol/3,3/ data imxchr/250/ c..... Restart real*8 tstop integer n_restart character*80 rstrt common /restart/ terminate logical terminate external sig_handler terminate =.false. inunnotset = .true. inunthreshold = 10 ! TODO: add to parameter file c.....track total model run time (wall clock time) time1 = time() call date_and_time (date(1),date(2),date(3),values) #ifdef PGI call signal(15, sig_handler, -1) ! for Portland Group #elif IFORT call signal(15, sig_handler, -1) ! for Intel ifort #elif XLF include 'fexcp.h' call signal(15, sig_handler) ! for Intel xlf #else call signal(15, sig_handler) ! for gfortran #endif c..... fast initialization of maximum value variable c..... (data-statement initializations are slow to compile) do i=1,nxc1 do j=1,nyc1 maxhc(i,j) = -9.9999998e+33 minhc(i,j) = 1.0e+34 maxspc(i,j) = -9.9999998e+33 qc1(i,j) = 0.0 inunc(i,j) = 0 end do end do do i=1,nxb1 do j=1,nyb1 maxhb(i,j) = -9.9999998e+33 minhb(i,j) = 1.0e+34 maxspb(i,j) = -9.9999998e+33 end do end do do i=1,nxa1 do j=1,nya1 maxha(i,j) = -9.9999998e+33 minha(i,j) = 1.0e+34 maxspa(i,j) = -9.9999998e+33 end do end do latimgA = .false. latimgB = .false. latimgC = .false. c.....track total model run time, date, revision 102 format ('MOST v2.5 ',a15,' ',i4.4,'-',i2.2,'-',i2.2,' ', & i2,':',i2.2,':',i2.2) c..... Don't edit the Rev: number below, svn does it for you c..... revision - simply make some change to this file, commit, update and make write (ga_hst,102) '$Rev: 2584 $ ', & values(1),values(2),values(3), & values(5),values(6),values(7) nargs = iargc() ds = ' ' if (nargs.NE.3 .and. nargs.NE.4) then write(*,104) 104 format('Usage ./most3_facts_nc ', & ' ', & '[directory structure check]') ga_len = len_trim(ga_hst) write(*,*) ga_hst(1:ga_len) stop else if (nargs.EQ.4) then c.......use operational directory structure if arg(4) = false call getarg(4, ds) if (ds.ne.'false'.and.ds.ne.'FALSE') then drstr = 1 else drstr = 0 endif else c.......use operational directory structure if nargs = 3 drstr = 0 endif call getarg(1, site) nc = lnblnk(site) call getarg(2, prefix) nc3 = lnblnk(prefix) call getarg(3, indir) nc2 = lnblnk(indir) if (nc2.eq.0) then nc2 = len(indir) write(3,*) ' directory is ',nc2,' characters too long!!' endif output = indir(1:nc2)//"output_"//site(1:nc)//".lis" c.....open text output file, clear old contents, reopen it open(unit=3,file=output,status='unknown',form='formatted') write(3,*) ' ' close(unit=3, status='keep') open(unit=3,file=output,status='unknown',form='formatted') c..... Restart filename (rstrt) rstrt=indir(1:nc2)//site(1:nc)//"_restart.nc" c.....track total model run time, date, revision ga_len = len_trim(ga_hst) write (3,*) ga_hst(1:ga_len) write(3,803) 803 format(' Site: ',$) write(3,*) site(1:nc) write(3,804) 804 format(' Input prefix: ',$) write(3,*) prefix(1:nc3) write(3,806) 806 format(' Input Directory: ',$) write(3,*) indir(1:nc2) c.....Get parameter inputs file c.....if 3 arguments drstr = false else 4th arg set: false or true c fname = ' ' if (drstr.EQ.1) then c.......input parameter file - read directly from run directory fname = indir(1:nc2)//'most3_facts_nc.in' write(3,*) 'Read Computational parameters from run directory: ', & fname(1:nc2+17) c.......create deformation prefix deform defprefix = 'deform ' else c.......input parameter file - read from site directory fname = indir(1:nc2)//site(1:nc)//'_run2d/most3_facts_nc.in' write(3,*) 'Read Computational parameters from site directory: ' & ,fname(1:nc+24+nc2) c.......create deformation prefix with site name defprefix = site(1:nc)//'_deform ' endif open(unit=2,file=fname,status='old',form='formatted',err=808) goto 809 808 write(3,*) 'ERORR: input param file not found: ', fname write(3,*) 'Halting execution.' stop 809 continue c..... Computational parameters write(3,801) 801 format(' Minimum amplitude of input offshore wave (m):',$) read (2,*) minwave write(3,*) minwave minwave=minwave*scl_m2cm write(3,811) 811 format(' Input minimum depth for offshore (m): ',$) read (2,*) grnd write(3,*) grnd write(3,8111) 8111 format (' Input "dry land" depth for inundation (m): ',$) read (2,*) dry write(3,*) dry write(3,8112) 8112 format (' Input friction coefficient (n**2): ',$) read (2,*) fr write(3,*) fr write(3,8113) 8113 format (' Input runup switch (0 - runup only in gridC, 1 - runup &in all grids): ',$) read(2,*)irunup write(3,*)irunup if (irunup.NE.0.AND.irunup.NE.1) then write(3,*)'ERROR: Incorrect runup flag value ',irunup write(3,*) 'Halting execution.' stop endif write(3,8114) 8114 format(' Max allowed eta (m): ',$) read (2,*) capmax write(3,*) capmax write(3,812) 812 format(' Input time step (sec): ',$) read (2,*) dt write(3,*) dt write(3,813) 813 format(' Input amount of steps: ',$) read (2,*) mfin write(3,*) mfin write(3,8131) 8131 format(' Compute "A" arrays every n-th time step, n=',$) read (2,*) ka write (3,*) ka write(3,8132) 8132 format(' Compute "B" arrays every n-th time step, n=',$) read (2,*) kb write (3,*) kb write(3,814) 814 format(' Input number of steps between snapshots (should be a mul &tiple of A,B and C time steps) : ',$) read (2,*) ksf write(3,*) ksf if (ksf/ka*ka.NE.ksf.OR.ksf/kb*kb.NE.ksf) then write(3,*)'ERROR: Output step ',ksf,' needs to be a multiple of & ',ka,' and ',kb write(3,*) 'Halting execution.' stop endif write(3,815) 815 format(' ...Starting from: ',$) read (2,*) nstart write(3,*) nstart write(3,816) 816 format(' ...Saving grid every n-th node, n=',$) read (2,*) npics write(3,*) npics c..... Read Bathymetry and initialize variables write (3,*) 'Reading Bathymetry' write (3,*) ' 1-ST LEVEL:' call BathRead(da,qa,h1a,h2a,ma,na,mma,nna,grnd,indir, & latimgA) write (3,*) ' 2-ND LEVEL:' call BathRead(db1,qb1,h1b1,h2b1,mb1,nb1,mmb1,nnb1,grnd,indir, & latimgB) write (3,*) ' 3-RD LEVEL:' call BathRead(dc1,qc1,h1c1,h2c1,mc1,nc1,mmc1,nnc1,dry,indir, & latimgC) c..... Read input/output directories write(3,600) 600 format(' DODS URL:',$) read (2,'(a)') durlin write(3,*) durlin(1:lnblnk(durlin)) write(3,805) 805 format(' output directory: ',$) read (2,'(a)') durlout write(3,*) durlout(1:lnblnk(durlout)) c..... If the terminate file exists, delete it open(unit=17,file=durlout(1:lnblnk(durlout))//'shootmenow', & status='old',err=601) write(3,*) 'WARNING: Found kill request file, deleting it' close(unit=17,status='delete') 601 continue c..... Read output suppression switches write(3,*) ' Produce/Suppress netCDF output for grids' write(3,850) 850 format(' A-grid B-grid C-grid Sift-file (default, 1 1 1 0): ',$) read(unit=2,fmt=*,end=855,err=855) suppressA,suppressB,suppressC, & suppressS write(3,*) suppressA,suppressB,suppressC,suppressS c..... Read Timeseries locations, if any numts = 0 read (unit=2, fmt=*, end=853, err=853) numts do i=1,numts read (unit=2,fmt=*,end=853,err=853) tsgrid(i),tsi(i),tsj(i) end do goto 854 853 write(3,*) 'Cannot read timeseries locations, not outputting ts' numts = 0 854 write(3,*) 'number of timeseries locations: ', numts do i=1,numts write(3,*) 'grid: ',tsgrid(i),' i: ',tsi(i),' j: ',tsj(i) enddo 855 continue c..... close model parameters file close(unit=2,status='keep') write(3,*) ' Done with input file ',fname call indx(qa,da,qb1,db1,h1a,h2a,h1b1,h2b1,ina,outa,crnra & ,x1lb1,x2rb1,y1lb1,y2rb1,ma,na,mb1,nb1 & ,mma,nna,mmb1,nnb1) call indx(qb1,db1,qc1,dc1,h1b1,h2b1,h1c1,h2c1,inb1,outb1,crnrb1 & ,x1lc1,x2rc1,y1lc1,y2rc1,mb1,nb1,mc1,nc1,mmb1,nnb1 & ,mmc1,nnc1) c..................... FACTS input files........................ nc4 = lnblnk(durlin) nc2 = lnblnk(prefix) if ((nc4+nc2) .ge. 120) then write(3,*) ' DODS URL name too long ',nc4+nc2,' characters' endif fname = durlin(1:nc4)//prefix(1:nc2)//' ' nc3 = lnblnk(fname) write(3,*) ' Input FACTS files:' write(3,610) 610 format(' zonal U: ',$) fname2 = fname(1:nc3)//'u.nc' write(3,*) fname2(1:nc3+4) if (CheckVar(fname2,ncid2in,'UA',UA_idin).eq.0) goto 9999 write(3,620) 620 format(' meridial V: ',$) fname3 = fname(1:nc3)//'v.nc' write(3,*) fname3(1:nc3+4) if (CheckVar(fname3,ncid3in,'VA',VA_idin).eq.0) goto 9999 write(3,630) 630 format(' amplitudes H: ',$) fname1 = fname(1:nc3)//'h.nc' write(3,*) fname1(1:nc3+4) if (CheckVar(fname1,ncid1in,'HA',HA_idin).eq.0) goto 9999 c get run command & add to history call getarg(0, tmp) nc4 = lnblnk(tmp) nc3 = lnblnk(indir) ga_len = len_trim(ga_hst) ga_hst = ga_hst(1:ga_len)//'\n'//tmp(1:nc4)//' '//site(1:nc)//' ' + //prefix(1:nc2)//' '//indir(1:nc3)//' '//ds c ga_len = len_trim(ga_hst) c write(3,*) ga_len,' history: ','*'//ga_hst(1:ga_len)//'*' c get title c iret = nf_inq_attlen(ncid1in,NF_GLOBAL,'title',ga_len) c iret = nf_get_att_text(ncid1in,NF_GLOBAL,'title',ga_ttl) write(ga_ttl,103) 103 format('MOST tsunami inundation model output ') ga_len = len_trim(ga_ttl) if (ga_len .gt. imxchr) then write(3,*) ' TITLE longer than',imxchr,' character *** truncated' endif c write(3,*) ' length:',ga_len,' title: ',ga_ttl(1:ga_len) call check_err(nf_inq(ncid1in,ndims,nvar,ngatts,ncunl_id)) c write(3,'("nf_inq(",5i3,")")') ncid1in,ndims,nvar,ngatts,ncunl_id do i=1,ndims call check_err(nf_inq_dim(ncid1in,i,dumb,inlen(i))) c write(3,*) ' nf_inq_dim(',ncid1in,i,dumb(1:lnblnk(dumb)), c + inlen(i),')' if ((dumb(1:3) .eq. 'TIM').or.(dumb(1:3) .eq. 'tim')) then nt0 = inlen(i) call check_err(nf_inq_varid(ncid1in,dumb(1:lnblnk(dumb)), + TIME_idin)) c write(3,*) ' nf_inq_varid(',ncid1in,dumb(1:lnblnk(dumb)), c + TIME_idin,')' else if ((dumb(1:3) .eq. 'LON').or.(dumb(1:3) .eq. 'lon')) then ma0 = inlen(i) call check_err(nf_inq_varid(ncid1in,dumb(1:lnblnk(dumb)), + LON_idin)) c write(3,*) ' nf_inq_varid(',ncid1in,dumb(1:lnblnk(dumb)), c + LON_idin,')' else if ((dumb(1:3) .eq. 'LAT').or.(dumb(1:3) .eq. 'lat')) then na0 = inlen(i) call check_err(nf_inq_varid(ncid1in,dumb(1:lnblnk(dumb)), + LAT_idin)) c write(3,*) ' nf_inq_varid(',ncid1in,dumb(1:lnblnk(dumb)), c + LAT_idin,')' endif end do c..... Read 0-grid lat/lon/time arrays call freadnc(ncid1in,TIME_idin,LON_idin,LAT_idin, & nt0,ma0,na0,h1a0,h2a0,t0) c..... Check that A-grid is within 0-grid if (h1a(1).lt.h1a0(1) .OR. h1a(ma).gt.h1a0(ma0) .OR. & h2a(1).lt.h2a0(1) .OR. h2a(na).gt.h2a0(na0)) then write(3,*) 'ERROR: 0-grid (prop) does not contain the A-grid' write(3,*) '0-grid extents, minLon: ',h1a0(1), & ' maxLon: ', h1a0(ma0) write(3,*) 'minLat: ',h2a0(1), ' maxLat: ',h2a0(na0) write(3,*) 'A-grid extents, minLon: ',h1a(1), & ' maxLon: ', h1a(ma) write(3,*) 'minLat: ',h2a(1), ' maxLat: ',h2a(na) write(3,*) 'Halting execution.' stop endif c..... Deformation input call readdeform(ncid1in,h1a0g,h2a0g,dx0,dy0,nx0g,ny0g, & ma,na,h1a,h2a,ddeform,defmax) write(3,*) 'defmax: ', defmax c..... Threshold for deforming bathymetry (0.001 m recommended) defthr=0.001 if (defmax.gt.defthr) then write(3,*) ' Deformation is over threshold, deforming bathies' c..... regrid deform onto a-grid call rgrd2(dx0,dy0,nx0g,ny0g,h1a0g,h2a0g,ddeform & ,ma,na,mma,nna,h1a,h2a,deformA & ,intpol,w,lw,iw,liw,ier) if (ier .ne. 0) then write(3,*) ' Error: regridding def onto A-grid: ', ier call exit(-1) endif c..... regrid deform onto b-grid call rgrd2(dx0,dy0,nx0g,ny0g,h1a0g,h2a0g,ddeform & ,mb1,nb1,mmb1,nnb1,h1b1,h2b1,deformB & ,intpol,w,lw,iw,liw,ier) if (ier .ne. 0) then write(3,*) ' Error: regridding def onto B-grid: ', ier call exit(-1) endif c..... regrid deform onto c-grid call rgrd2(dx0,dy0,nx0g,ny0g,h1a0g,h2a0g,ddeform & ,mc1,nc1,mmc1,nnc1,h1c1,h2c1,deformC & ,intpol,w,lw,iw,liw,ier) if (ier .ne. 0) then write(3,*) ' Error: regridding def onto C-grid: ', ier call exit(-1) endif c..... create output files nc2 = lnblnk(defprefix) fname = defprefix(1:nc2)//'a.dat' do j=1,na do i=1,ma da(i,j)=da(i,j)+deformA(i,j) end do end do write(3,*) ' Writing deformed A-grid bathy file: ',fname call bath_write(fname,ma,na,mma,nna,h1a,h2a,da) fname = defprefix(1:nc2)//'b.dat' do j=1,nb1 do i=1,mb1 db1(i,j)=db1(i,j)+deformB(i,j) end do end do write(3,*) ' Writing deformed B-grid bathy file: ',fname call bath_write(fname,mb1,nb1,mmb1,nnb1,h1b1,h2b1,db1) fname = defprefix(1:nc2)//'c.dat' do j=1,nc1 do i=1,mc1 dc1(i,j)=dc1(i,j)+deformC(i,j) end do end do write(3,*) ' Writing deformed C-grid bathy file: ',fname call bath_write(fname,mc1,nc1,mmc1,nnc1,h1c1,h2c1,dc1) write(3,*) ' ' 666 format (3000F8.2) else write(3,*) ' no local deformation detected' write(3,*) ' ' endif c..... if error in local deformation checking, continue from here 640 continue write (3,*) 'size of input array:', ma0,na0,nt0 write (3,*) ' Longitude:',h1a0(1),' to',h1a0(ma0) write (3,*) ' Latitude:',h2a0(1),' to',h2a0(na0) write (3,*) ' Time:',t0(1),' to',t0(nt0) call flush(3) call indx_off(h1a0,h2a0,h1a,h2a,ina0,outa0,crnra0 & ,ma0,na0,ma,na) c..... The initial number of inundated points (should be zero) call countinun(qc1,dc1,inunc,numinun,mc1,nc1,nxc1,nyc1,dry) write(3,*) 'Number of inundated points in c-grid: ', numinun do i=1,ma0 do j=1,na0 da0(i,j)=0. end do end do dt0=(t0(2)-t0(1)) t=t0(1) it0=1 cm=0.01 taper=1.0 nowave=.TRUE. c..... adding switch for first time first = .TRUE. call flush(3) c.....next lines until 'end4restart' are for restart call read4restart(rstrt,'ua',n_restart,tstop,icnt, & nxa1,nya1,ma,na,ua) if(n_restart.gt.0) then call read4restart(rstrt,'va',n_restart,tstop,icnt, & nxa1,nya1,ma,na,va) call read4restart(rstrt,'qa',n_restart,tstop,icnt, & nxa1,nya1,ma,na,qa) call read4restart(rstrt,'ub1',n_restart,tstop,icnt, & nxb1,nyb1,mb1,nb1,ub1) call read4restart(rstrt,'vb1',n_restart,tstop,icnt, & nxb1,nyb1,mb1,nb1,vb1) call read4restart(rstrt,'qb1',n_restart,tstop,icnt, & nxb1,nyb1,mb1,nb1,qb1) call read4restart(rstrt,'uc1',n_restart,tstop,icnt, & nxc1,nyc1,mc1,nc1,uc1) call read4restart(rstrt,'vc1',n_restart,tstop,icnt, & nxc1,nyc1,mc1,nc1,vc1) call read4restart(rstrt,'qc1',n_restart,tstop,icnt, & nxc1,nyc1,mc1,nc1,qc1) if(irunup.eq.0) then call bounds(qb1,ub1,vb1,db1,h1b1,h2b1,h1c1,h2c1,x1rc1,x2lc1 & ,y1rc1,y2lc1,crnrb1,inb1,mb1,nb1,mc1,nc1,mmc1,nnc1,mmb1,nnb1) call bounds(qa,ua,va,da,h1a,h2a,h1b1,h2b1,x1rb1,x2lb1,y1rb1 & ,y2lb1,crnra,ina,ma,na,mb1,nb1,mmb1,nnb1,mma,nna) else call bounds_nr(qa,ua,va,da,h1a,h2a,h1b1,h2b1,x1rb1,x2lb1,y1rb1 & ,y2lb1,crnra,ina,ma,na,mb1,nb1,mmb1,nnb1,mma,nna,dry) call bounds_nr(qb1,ub1,vb1,db1,h1b1,h2b1,h1c1,h2c1,x1rc1,x2lc1 & ,y1rc1,y2lc1,crnrb1,inb1,mb1,nb1,mc1,nc1,mmc1,nnc1,mmb1,nnb1 & ,dry) endif endif c.....end4restart c...........Initialize NetCDF output files for all 3 grids.............. nx_nc = (mc1-1)/npics+1 ny_nc = (nc1-1)/npics+1 nx_nb = (mb1-1)/npics+1 ny_nb = (nb1-1)/npics+1 nx_na = (ma-1)/npics+1 ny_na = (na-1)/npics+1 write(3,*) 'NetCDF array size for grid C: ',nx_nc,ny_nc write(3,*) 'NetCDF array size for grid B: ',nx_nb,ny_nb write(3,*) 'NetCDF array size for grid A: ',nx_na,ny_na c..... create output files nc4 = lnblnk(durlout) prefix = durlout(1:nc4)//site(1:nc)//'_runup ' nc2 = lnblnk(prefix) 818 format('creating netCDF output file: ',A) if (suppressC.ne.0) then fname = prefix(1:nc2)//'_ha.nc' nvar = 1 call fgennc(h1c1,h2c1,fname,nvar,n_restart,ncid1C,TIME_idC, & HA_idC, MAXH_idC, nx_nc,ny_nc,mc1,nc1,ga_hst,ga_ttl,latimgC) write(3,818) fname(1:nc2+6) nvar = 2 fname = prefix(1:nc2)//'_ua.nc' call fgennc(h1c1,h2c1,fname,nvar,n_restart,ncid2C,TIME_idC, & UA_idC, MAXU_idC,nx_nc,ny_nc,mc1,nc1,ga_hst,ga_ttl,latimgC) write(3,818) fname(1:nc2+6) nvar = 3 fname = prefix(1:nc2)//'_va.nc' call fgennc(h1c1,h2c1,fname,nvar,n_restart,ncid3C,TIME_idC, & VA_idC, MAXV_idC, nx_nc,ny_nc,mc1,nc1,ga_hst,ga_ttl,latimgC) write(3,818) fname(1:nc2+6) endif if (suppressB.ne.0) then prefix = durlout(1:nc4)//site(1:nc)//'_runupB ' nc2 = lnblnk(prefix) fname = prefix(1:nc2)//'_ha.nc' nvar = 1 call fgennc(h1b1,h2b1,fname,nvar,n_restart,ncid1B,TIME_idB, & HA_idB, MAXH_idB, nx_nb,ny_nb,mb1,nb1,ga_hst,ga_ttl,latimgB) write(3,818) fname(1:nc2+6) nvar = 2 fname = prefix(1:nc2)//'_ua.nc' call fgennc(h1b1,h2b1,fname,nvar,n_restart,ncid2B,TIME_idB, & UA_idB, MAXU_idB, nx_nb,ny_nb,mb1,nb1,ga_hst,ga_ttl,latimgB) write(3,818) fname(1:nc2+6) nvar = 3 fname = prefix(1:nc2)//'_va.nc' call fgennc(h1b1,h2b1,fname,nvar,n_restart,ncid3B,TIME_idB, & VA_idB, MAXV_idB, nx_nb,ny_nb,mb1,nb1,ga_hst,ga_ttl,latimgB) write(3,818) fname(1:nc2+6) end if if (suppressA.ne.0) then prefix = durlout(1:nc4)//site(1:nc)//'_runupA ' nc2 = lnblnk(prefix) fname = prefix(1:nc2)//'_ha.nc' nvar = 1 call fgennc(h1a,h2a,fname,nvar,n_restart,ncid1A,TIME_idA, & HA_idA, MAXH_idA, nx_na,ny_na,ma,na,ga_hst,ga_ttl,latimgA) write(3,818) fname(1:nc2+6) nvar = 2 fname = prefix(1:nc2)//'_ua.nc' call fgennc(h1a,h2a,fname,nvar,n_restart,ncid2A,TIME_idA, & UA_idA, MAXU_idA, nx_na,ny_na,ma,na,ga_hst,ga_ttl,latimgA) write(3,818) fname(1:nc2+6) nvar = 3 fname = prefix(1:nc2)//'_va.nc' call fgennc(h1a,h2a,fname,nvar,n_restart,ncid3A,TIME_idA, & VA_idA, MAXV_idA, nx_na,ny_na,ma,na,ga_hst,ga_ttl,latimgA) write(3,818) fname(1:nc2+6) end if if (suppressS.ne.0) then prefix = durlout(1:nc4)//site(1:nc)//'_sift ' nc2 = lnblnk(prefix) fnsift = prefix(1:nc2)//'.nc' if (n_restart.gt.0 .AND. & nf_open(fnsift, NF_WRITE, siftid).eq.NF_NOERR) then write(3,*) 'reading from: ', fnsift(1:lnblnk(fnsift)) call flush(3) call siftread(siftid,tstimeid,tsid, & mxaid,mnaid,spaid,mxbid,mnbid,spbid, & mxcid,mncid,spcid, & maxha,minha,maxspa,ma,na,nxa1,nya1, & maxhb,minhb,maxspb,mb1,nb1,nxb1,nyb1, & maxhc,minhc,maxspc,mc1,nc1,nxc1,nyc1,numts,floodid, & latimgA,latimgB,latimgC) else write(3,818) fnsift(1:lnblnk(fnsift)) call flush(3) call siftcreate(fnsift, & h1a,h2a,ma,na,nxa1,nya1, & h1b1,h2b1,mb1,nb1,nxb1,nyb1, & h1c1,h2c1,mc1,nc1,nxc1,nyc1, & numts,tsgrid,tsi,tsj,siftid,tstimeid,tsid, & mxaid,mnaid,spaid,mxbid,mnbid,spbid,mxcid,mncid,spcid, & floodid,inunthreshold,latimgA,latimgB,latimgC) end if end if write(3,*) ' netCDF initialization complete' call flush(3) c..... New algorithm for starting calculations: c..... threshold = greater of: c..... 1.0% of max amp of incoming wave .OR. minwave c..... Calculate max amp of incoming wave on A-Grid boundary: maxa0 = 0.0 do it0=1,nt0 call readrecs(ncid1in,HA_idin,it0,rHA,qa01,ma0,na0,mma0,nna0) do i=1,ma0 do j=1,na0 if (qa01(i,j).GT.maxa0.AND. & i.GE.ina0(1,1).AND.i.LE.ina0(2,1).AND. & j.GE.ina0(3,1).AND.j.LE.ina0(4,1)) then maxa0=qa01(i,j) itcalc=it0 end if end do end do end do write(3,*) 'incoming wave max amp on A-Grid boundary: ',maxa0 maxa0=0.01*maxa0 if (maxa0 .GT. minwave) then minwave=maxa0 write(3,*) 'setting threshold to 1% of max amp: ',minwave else write(3,*) 'leaving threshold at: ',minwave,' cm' end if call flush(3) it0=1 do i=1,ma0 do j=1,na0 qa01(i,j)=0. end do end do c..... Time Cycle maxerror = 0.0 c.... icnt = 0 ! get from restart file do 888 nstep = 1,mfin ! begins time cycle c..... Reading input wave do while ((nowave.OR.t0(it0).LE.t+dt0).AND.(it0.LT.nt0)) do i=1,ma0 do j=1,na0 ua00(i,j)=ua01(i,j) va00(i,j)=va01(i,j) qa00(i,j)=qa01(i,j) if (nowave) then if (abs(qa01(i,j)).GT.minwave.AND. & i.GE.ina0(1,1).AND.i.LE.ina0(2,1).AND. & j.GE.ina0(3,1).AND.j.LE.ina0(4,1)) then nowave=.FALSE. t=t0(it0) write(3,*) & 'input',it0,' wave detected at',t0(it0),' amp:',qa01(i,j),'cm' &,' at ',h1a0(i),' , ',h2a0(j) call flush(3) end if end if end do end do c..... interpolating initial conditions if (.NOT.nowave.AND.first) then first = .FALSE. c..... if t > tstop (restart stop time) if (t.GT.tstop) then call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,ua01 & ,ma,na,mma,nna,h1a,h2a,ua & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,va01 & ,ma,na,mma,nna,h1a,h2a,va & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,qa01 & ,ma,na,mma,nna,h1a,h2a,qa & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,ua01 & ,mb1,nb1,mmb1,nnb1,h1b1,h2b1,ub1 & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,va01 & ,mb1,nb1,mmb1,nnb1,h1b1,h2b1,vb1 & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,qa01 & ,mb1,nb1,mmb1,nnb1,h1b1,h2b1,qb1 & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,ua01 & ,mc1,nc1,mmc1,nnc1,h1c1,h2c1,uc1 & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,va01 & ,mc1,nc1,mmc1,nnc1,h1c1,h2c1,vc1 & ,intpol,w,lw,iw,liw,ier) call rgrd2(ma0,na0,mma0,nna0,h1a0,h2a0,qa01 & ,mc1,nc1,mmc1,nnc1,h1c1,h2c1,qc1 & ,intpol,w,lw,iw,liw,ier) if (ier.NE.0) then write (3,*) 'ERROR: **** Regridding error ', ier goto 9999 endif write(3,*) 'Initial surface is read at t=',t0(it0) call flush(3) do i=1,ma do j=1,na ua(i,j)=ua(i,j)*cm va(i,j)=va(i,j)*cm qa(i,j)=qa(i,j)*cm+da(i,j) end do end do do i=1,mb1 do j=1,nb1 ub1(i,j)=ub1(i,j)*cm vb1(i,j)=vb1(i,j)*cm qb1(i,j)=qb1(i,j)*cm+db1(i,j) end do end do do i=1,mc1 do j=1,nc1 uc1(i,j)=uc1(i,j)*cm vc1(i,j)=vc1(i,j)*cm qc1(i,j)=qc1(i,j)*cm+dc1(i,j) end do end do endif ! end of if added for restart cA Adjusting bathymetry for grid B and grid A runs if(irunup.eq.1) then do i = 1,ma do j = 1,na if (da(i,j).LE.grnd) then da(i,j) = -500.0 qa(i,j) = -500.0 endif end do end do do i = 1,mb1 do j = 1,nb1 if (db1(i,j).LE.grnd) then db1(i,j) = -500.0 qb1(i,j) = -500.0 endif end do end do end if end if ! .NOT.nowave.AND.first call readrecs(ncid2in,UA_idin,it0,rUA,ua01,ma0,na0,mma0,nna0) call readrecs(ncid3in,VA_idin,it0,rVA,va01,ma0,na0,mma0,nna0) call readrecs(ncid1in,HA_idin,it0,rHA,qa01,ma0,na0,mma0,nna0) it0=it0+1 c if (it0.EQ.nt0) cm=0.0 CC uncomment for more verbose .lis file CC if (t0(it0).GT.t+dt0) then CC write (3,*) 'input',it0-1,' is read at t=',t0(it0-1) CC end if end do ! Input wave is read t=t+dt if (nowave) then write(3,*) 'WARNING: input wave smaller than trigger threshold' write(3,*) ' No model output produced.' goto 9999 endif c..... if t < tstop (restart time), dont compute if (t.le.tstop) goto 888 c write(*,*)'At step = ',nstep call TimeStep(qc1,uc1,vc1,dc1,qw,uw,vw,dw,q1,u1,v1,d1,h,ndum & ,x1rc1,x2lc1,y1rc1,y2lc1,h1c1,h2c1 & ,mc1,nc1,mmc1,nnc1,dt) c.... Compute maxh(i,j) and maxspeed(i,j) for _sift.nc output 645 format('WARNING: num inundated pts GT ',i3,' at timestep ',i5) if (suppressS .ne. 0) then call MaxSlab(qc1,dc1,maxhc,minhc,mc1,nc1,nxc1,nyc1,dry) call MaxSpeed(uc1,vc1,qc1,maxspc,mc1,nc1,nxc1,nyc1,dry) c..... Check for inundation, set flag if (nstep .gt. 3) then ! skip first few steps call countinun(qc1,dc1,inunc,numinun,mc1,nc1,nxc1,nyc1,dry) if (inunnotset .and. numinun .gt. inunthreshold) then inunnotset = .false. call writeinunflag(siftid,floodid,icnt) write(3,645) inunthreshold, icnt endif endif endif c.... B grid if((nstep-1)/kb*kb.EQ.nstep-1)then if(irunup.eq.0) then call TimeStep_n(qb1,ub1,vb1,db1,qw,uw,vw,dw,q1,u1,v1,d1,h & ,ndum,x1rb1,x2lb1,y1rb1,y2lb1 & ,h1b1,h2b1,mb1,nb1,mmb1,nnb1,dt*kb) call bounds(qb1,ub1,vb1,db1,h1b1,h2b1,h1c1,h2c1,x1rc1,x2lc1 & ,y1rc1,y2lc1,crnrb1,inb1,mb1,nb1,mc1,nc1,mmc1,nnc1,mmb1,nnb1) elseif(irunup.eq.1) then call TimeStep_nr(qb1,ub1,vb1,db1,qw,uw,vw,dw,q1,u1,v1,d1,h & ,ndum,x1rb1,x2lb1,y1rb1,y2lb1 & ,h1b1,h2b1,mb1,nb1,mmb1,nnb1,dt*kb) call bounds_nr(qb1,ub1,vb1,db1,h1b1,h2b1,h1c1,h2c1,x1rc1,x2lc1 & ,y1rc1,y2lc1,crnrb1,inb1,mb1,nb1,mc1,nc1,mmc1,nnc1,mmb1,nnb1 & ,dry) endif if(suppressS .ne. 0) then call MaxSlab(qb1,db1,maxhb,minhb,mb1,nb1,nxb1,nyb1,grnd) call MaxSpeed(ub1,vb1,qb1,maxspb,mb1,nb1,nxb1,nyb1,grnd) endif end if c.... A grid if((nstep-1)/ka*ka.EQ.nstep-1)then c.....Check to see if there are any more forcing timesteps. No? taper if (it0.GE.nt0-18) then taper = taper * 0.9 if (taper.lt. 1.0) then write(3,*) 'WARNING: reached end of forcing file...' write(3,*) ' tapering last forcing frame.' end if end if do i=1,ma0 ! adjust boundary input wave do j=1,na0 ua0(i,j)=(ua00(i,j) & +(ua01(i,j)-ua00(i,j))*(t-t0(it0-2))/dt0)*cm*taper va0(i,j)=(va00(i,j) & +(va01(i,j)-va00(i,j))*(t-t0(it0-2))/dt0)*cm*taper qa0(i,j)=(qa00(i,j) & +(qa01(i,j)-qa00(i,j))*(t-t0(it0-2))/dt0)*cm*taper end do end do call bounds(qa0,ua0,va0,da0,h1a0,h2a0,h1a,h2a,x1ra,x2la,y1ra & ,y2la,crnra0,ina0,ma0,na0,ma,na,mma,nna,mma0,nna0) c write(*,*)'Computing gridA ..' if(irunup.eq.0) then call TimeStep_n(qa,ua,va,da,qw,uw,vw,dw,q1,u1,v1,d1,h,ndum & ,x1ra,x2la,y1ra,y2la & ,h1a,h2a,ma,na,mma,nna,dt*ka) call bounds(qa,ua,va,da,h1a,h2a,h1b1,h2b1,x1rb1,x2lb1,y1rb1 & ,y2lb1,crnra,ina,ma,na,mb1,nb1,mmb1,nnb1,mma,nna) elseif(irunup.eq.1) then call TimeStep_nr(qa,ua,va,da,qw,uw,vw,dw,q1,u1,v1,d1,h,ndum & ,x1ra,x2la,y1ra,y2la & ,h1a,h2a,ma,na,mma,nna,dt*ka) call bounds_nr(qa,ua,va,da,h1a,h2a,h1b1,h2b1,x1rb1,x2lb1,y1rb1 & ,y2lb1,crnra,ina,ma,na,mb1,nb1,mmb1,nnb1,mma,nna,dry) endif if(suppressS .ne. 0) then call MaxSlab(qa,da,maxha,minha,ma,na,nxa1,nya1,grnd) call MaxSpeed(ua,va,qa,maxspa,ma,na,nxa1,nya1,grnd) endif endif c Check for maximum value over limit (blowup) call MaxValue(qc1,dc1,mc1,nc1,mmc1,nnc1,dry,maxc, & minc,mxci,mxcj,mnci,mncj) call MaxValue(qb1,db1,mb1,nb1,mmb1,nnb1,grnd,maxb, & minb,mxbi,mxbj,mnbi,mnbj) call MaxValue(qa,da,ma,na,mma,nna,grnd,maxa, & mina,mxai,mxaj,mnai,mnaj) if (maxc .gt. capmax) then write(3,*) 'ERROR: Maximum exceeded in C grid' write(3,*) 'blowup RUN STOPPED maxerror = ',maxc write(3,*) 'occurred at i = ', mxci, ' x = ',h1c1(mxci) write(3,*) 'occurred at j = ', mxcj, ' y = ',h2c1(mxcj) goto 9999 endif if (maxb .gt. capmax) then write(3,*) 'ERROR: Maximum exceeded in B grid' write(3,*) 'blowup RUN STOPPED maxerror = ',maxb write(3,*) 'occurred at i = ', mxbi, ' x = ',h1b1(mxbi) write(3,*) 'occurred at j = ', mxbj, ' y = ',h2b1(mxbj) goto 9999 endif if (maxa .gt. capmax) then write(3,*) 'ERROR: Maximum exceeded in A grid' write(3,*) 'blowup RUN STOPPED maxerror = ',maxa write(3,*) 'occurred at i = ', mxai, ' x = ',h1a(mxai) write(3,*) 'occurred at j = ', mxaj, ' y = ',h2a(mxaj) goto 9999 endif c..... Write into NetCDF file if (nstep.GE.nstart) then C uncomment to get maximum values at each timestep (in ascii) c call SurfWrite('cg',qc1,dc1,h1c1,h2c1,mc1,nc1,mmc1,nnc1,t, c & 1,dax,maxic,dry) c call SurfWrite('bg',qb1,db1,h1b1,h2b1,mb1,nb1,mmb1,nnb1,t c & 1,dax,maxib,grnd) c call SurfWrite('ag',qa,da,h1a,h2a,ma,na,mma,nna,t c & 1,dax,maxia,grnd) C end of maximum value at each timestep C If the terminate file exists, save restart file, end program open(unit=17,file=durlout(1:lnblnk(durlout))//'shootmenow', & status='old',err=401) write(3,*) 'WARNING: Found kill request file, terminating run' c close(unit=17,status='delete') goto 9990 401 continue if((nstep-nstart)/ksf*ksf.EQ.(nstep-nstart))then icnt = icnt+1 if((nstep-nstart)/(ksf*10)*(ksf*10).EQ.(nstep-nstart)) then write(3,402)icnt,t,t/60.,t/3600.,100*real(nstep)/real(mfin) 402 format('Output time step: ',i8,' time: ',f10.1,' sec, ', & f8.1,' min, ', f4.1,' hrs,', f5.1,' % comp') 403 format(' Max/Min elevation in grid C are: ',f12.8,' /',f12.8) 404 format(' Max/Min elevation in grid B are: ',f12.8,' /',f12.8) 405 format(' Max/Min elevation in grid A are: ',f12.8,' /',f12.8) write(3,403) maxc, minc write(3,404) maxb, minb write(3,405) maxa, mina call flush(3) end if if (suppressC .ne. 0) then call writerecs(ncid1C,TIME_idC,HA_idC,1, & icnt,cHA,qc1,uc1,vc1,dc1 & ,nx_nc,ny_nc,nxc1,nyc1,tsi,tsj,0,tsid,latimgC) call writerecs(ncid2C,TIME_idC,UA_idC,2, & icnt,cUA,qc1,uc1,vc1,dc1 & ,nx_nc,ny_nc,nxc1,nyc1,tsi,tsj,0,tsid,latimgC) call writerecs(ncid3C,TIME_idC,VA_idC,3, & icnt,cVA,qc1,uc1,vc1,dc1 & ,nx_nc,ny_nc,nxc1,nyc1,tsi,tsj,0,tsid,latimgC) endif if (suppressB .ne. 0) then call writerecs(ncid1B,TIME_idB,HA_idB,1, & icnt,bHA,qb1,ub1,vb1,db1 & ,nx_nb,ny_nb,nxb1,nyb1,tsi,tsj,0,tsid,latimgB) call writerecs(ncid2B,TIME_idB,UA_idB,2, & icnt,bUA,qb1,ub1,vb1,db1 & ,nx_nb,ny_nb,nxb1,nyb1,tsi,tsj,0,tsid,latimgB) call writerecs(ncid3B,TIME_idB,VA_idB,3, & icnt,bVA,qb1,ub1,vb1,db1 & ,nx_nb,ny_nb,nxb1,nyb1,tsi,tsj,0,tsid,latimgB) endif if (suppressA .ne. 0) then call writerecs(ncid1A,TIME_idA,HA_idA,1, & icnt,aHA,qa,ua,va,da & ,nx_na,ny_na,nxa1,nya1,tsi,tsj,0,tsid,latimgA) call writerecs(ncid2A,TIME_idA,UA_idA,2, & icnt,aUA,qa,ua,va,da & ,nx_na,ny_na,nxa1,nya1,tsi,tsj,0,tsid,latimgA) call writerecs(ncid3A,TIME_idA,VA_idA,3, & icnt,aVA,qa,ua,va,da & ,nx_na,ny_na,nxa1,nya1,tsi,tsj,0,tsid,latimgA) endif if (suppressS .ne. 0) then call siftsave(siftid,tstimeid,icnt,t, & qa,da,nxa1,nya1,ma,na,qb1,db1,nxb1,nyb1,mb1,nb1, & qc1,dc1,nxc1,nyc1,mc1,nc1, & tsgrid,tsi,tsj,numts,tsid,floodid,inunc,h1c1,h2c1, & latimgA,latimgB,latimgC,dry) call siftsavemax(siftid,mxaid,mnaid,spaid,ma,na,nxa1,nya1, & maxha, minha, maxspa,latimgA) call siftsavemax(siftid,mxbid,mnbid,spbid,mb1,nb1,nxb1,nyb1, & maxhb, minhb, maxspb,latimgB) call siftsavemax(siftid,mxcid,mncid,spcid,mc1,nc1,nxc1,nyc1, & maxhc, minhc, maxspc,latimgC) iret = nf_sync(siftid) call check_err(iret) endif end if end if if (terminate) then write(3,*) ' WARNING: Terminated before run completed...' call flush(3) goto 9990 endif 888 continue ! end of time cycle 9990 continue call countinun(qc1,dc1,inunc,numinun,mc1,nc1,nxc1,nyc1,dry) 455 format(a,i12) write(3,455)' Total number of inundated points: ', numinun c.... save restart c.... (note the antics to keep from exceeding stack limits) write(3,*) ' saving restart file' call save4restart(rstrt,'ua',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mma,nna,ma,na,ua) call save4restart(rstrt,'va',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mma,nna,ma,na,va) call save4restart(rstrt,'qa',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mma,nna,ma,na,qa) call save4restart(rstrt,'ub1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmb1,nnb1,mb1,nb1,ub1) call save4restart(rstrt,'vb1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmb1,nnb1,mb1,nb1,vb1) call save4restart(rstrt,'qb1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmb1,nnb1,mb1,nb1,qb1) call save4restart(rstrt,'uc1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmc1,nnc1,mc1,nc1,uc1) call save4restart(rstrt,'vc1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmc1,nnc1,mc1,nc1,vc1) call save4restart(rstrt,'qc1',n_restart,t,icnt, & mma,nna,ma,na,mmb1,nnb1,mb1,nb1,mmc1,nnc1,mc1,nc1, & mmc1,nnc1,mc1,nc1,qc1) write(3,*) ' done saving restart file' call flush(3) if (suppressS .ne. 0) then write(3,*) ' calculating maximums' call flush(3) call MaxSlab(qc1,dc1,maxhc,minhc,mc1,nc1,nxc1,nyc1,dry) call MaxSpeed(uc1,vc1,qc1,maxspc,mc1,nc1,nxc1,nyc1,dry) call MaxSlab(qb1,db1,maxhb,minhb,mb1,nb1,nxb1,nyb1,grnd) call MaxSpeed(ub1,vb1,qb1,maxspb,mb1,nb1,nxb1,nyb1,grnd) call MaxSlab(qa,da,maxha,minha,ma,na,nxa1,nya1,grnd) call MaxSpeed(ua,va,qa,maxspa,ma,na,nxa1,nya1,grnd) write(3,*) ' saving maximums to _sift.nc' call flush(3) call siftsave(siftid,tstimeid,icnt,t, & qa,da,nxa1,nya1,ma,na,qb1,db1,nxb1,nyb1,mb1,nb1, & qc1,dc1,nxc1,nyc1,mc1,nc1, & tsgrid,tsi,tsj,numts,tsid,floodid,inunc,h1c1,h2c1, & latimgA,latimgB,latimgC,dry) call siftsavemax(siftid,mxaid,mnaid,spaid,ma,na,nxa1,nya1, & maxha, minha, maxspa,latimgA) call siftsavemax(siftid,mxbid,mnbid,spbid,mb1,nb1,nxb1,nyb1, & maxhb, minhb, maxspb,latimgB) call siftsavemax(siftid,mxcid,mncid,spcid,mc1,nc1,nxc1,nyc1, & maxhc, minhc, maxspc,latimgC) iret = nf_sync(siftid) endif 9999 continue iret = nf_close(ncid1in) iret = nf_close(ncid2in) iret = nf_close(ncid3in) iret = nf_close(ncid1C) iret = nf_close(ncid2C) iret = nf_close(ncid3C) iret = nf_close(ncid1B) iret = nf_close(ncid2B) iret = nf_close(ncid3B) iret = nf_close(ncid1A) iret = nf_close(ncid2A) iret = nf_close(ncid3A) iret = nf_close(siftid) eee = etime(eeet) time2 = time() call date_and_time (date(1),date(2),date(3),values) c..... Don't edit the Rev: number below, svn does it for you c..... Revision - simply make some change to this file, commit, update and make write (3,102) '$Rev: 2584 $ ', & values(1),values(2),values(3), & values(5),values(6),values(7) write(3,*) ' elapsed secs:',eee,', user:',eeet(1),', sys:',eeet(2) write(3,*) ' clock time:',time2-time1, & 'secs, ',(time2-time1)/60.,' minutes' write(3,*) ' Run finished' call flush(3) close(unit=3,status='keep') call exit(0) end