!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module initexit use GO, only : gol, goErr, goPr implicit none ! --- in/out ------------------------------ private public :: default, exitus, start ! --- const -------------------------------- character(len=*), parameter :: mname = 'initexit' contains subroutine default !--------------------------------------------------------------------- ! set default values of control variables ! ! mh, 27-jun-1989 ! mh, 26-sep-1992 ! v 9.1 ! mk, 21-dec-2002 !--------------------------------------------------------------------- use dims use datetime, only : tau2date use toolbox, only : escape_tm implicit none ! local integer :: i,k,n ! start ! ! set default calendar ! ! (real calendar, iyear0=1980) icalendo=2 ! ! default time steps of basic tasks nstep=0 ndyn=3600*3 nconv=3600*3 ndiff=24*3600*1000 !never happes in one month ntrans=0 ndiag=4*3600 nchem=0 nsrce=24*3600 !nread=3*3600 ! <--- set in main program by Meteo_Init nwrite=-1 ninst=0 !c ! default is restart istart=10 itaui=0 newsrun=.true. call tau2date(itaui,idatei) if ( mod(idatei(4),3) /= 0 ) & call escape_tm('default: GMT start time should be multiple of 3') itaue=itaui call tau2date(itaue,idatee) itaut=0 call tau2date(itaut,idatet) !c ! output for conservation diagnostics !c ! -1: daily, -2: monthly, -3: yearly, !c ! >=0 interval in sec ndiagp1=-2 !c ! output for mean field diagnostics !c ! -1: daily, -2: monthly, -3: yearly, !c ! >=0 interval in sec ndiagp2=-2 !c ! full convection czeta=1. !c ! full vertical diffusion czetak=1. !c ! scaling factor for horizontal diffusion limits=.true. !c ! checking interval ncheck=6 !c ! control for debug output cdebug=.false. okdebug=.false. revert = 1 end subroutine default !-------------------------------------------------------------- ! terminate a model run ! ! icode=0 normal termination ! else abnormal termination ! ! msg character string to be printed on unit kmain ! ! mh, 29-jun-1989 ! ! changed to v 9.1 mh, 26-sep-1992 ! !-------------------------------------------------------------- subroutine exitus use GO , only : gol, goErr, goPr use GO , only : goGetFU use dims , only : nregions use dims , only : nstep, itau use dims , only : cpu0 use dims , only : nxi, xi, nloop_max use dims , only : okdebug, kmain use chem_param , only : ntrace_chem use tracer_data , only : chem_dat use global_data , only : free_fields use io_save, only : savehdf use restart, only : rs_write use user_output, only : User_output_Done !use user_output, only : user_output_done use sources_sinks, only : trace_end use datetime, only : tstamp use ParTools , only : myid, root, Par_Barrier, Par_StopMPI use tracer_data , only : Par_Check_Domain use global_data , only : outdir use advectm_cfl, only : done_cfl #ifdef with_budgets use budget_global, only : done_budget_global #endif ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/exitus' ! --- local ------------------------------ integer :: status integer :: fu integer :: region, gatheri real :: cpu3, gatherr ! --- begin --------------------------------------- do region=1,nregions ! transfer data to tracer domain call Par_Check_Domain( region, 'n', 'tracer' ) end do #ifndef with_prism ! save the model state. ! in case of coupled system, do not use this, as ! restart-files are used. This routine is quite memory-consuming, ! and should not be called if not used. if(.not.rs_write) then call SaveHDF( 'successful completion of run',trim(outdir)//'/save', status ) if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if else write (gol,'("in ",a,"; no writing of save.hdf files since restart files are used")') rname; call goPr write (gol,'("in ",a,"; ....(restart.write : T or with_prism in rc-file)")') rname; call goPr endif #endif !>>> moved to TM5/TM5_Model_Run !! complete user-specified output !CMK only once! !call user_output_done( status ) !if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if !<<< #ifdef with_budgets ! save budgets, print sumary call done_budget_global ( status ) if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if #endif ! finish up sources/sinks: call trace_end( status ) if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if ! clear short-lived tracers: if ( ntrace_chem > 0 ) then do region = 1, nregions deallocate( chem_dat(region)%rm_k ) end do end if call Par_Barrier if ( okdebug ) then write (gol,'("exitus: starting to free fields ",i4)') myid; call goPr end if call free_fields !free memory of this run call Par_Barrier if ( okdebug ) then write (gol,'("exitus: starting time computation",i4)') myid; call goPr end if call Par_Barrier if ( myid == root ) then write (gol,'(" ")'); call goPr write (gol,'("exitus: CFL info:")'); call goPr do region=1, nregions write (gol,'(a,3i4,f10.4)') ' exitus: x: region, nxi, nloop_max, xi', & region, nxi(region,1), nloop_max(region,1), xi(region,1); call goPr write (gol,'(a,3i4,f10.4)') ' exitus: y: region, nxi, nloop_max, xi', & region, nxi(region,2), nloop_max(region,2), xi(region,2); call goPr write (gol,'(a,3i4,f10.4)') ' exitus: z: region, nxi, nloop_max, xi', & region, nxi(region,3), nloop_max(region,3), xi(region,3); call goPr end do end if ! cfl finished ... call Done_CFL call Par_Barrier if ( myid == root ) then write (gol,'(1x)'); call goPr write (gol,'(1x)'); call goPr write (gol,'(1x,39("--"))'); call goPr write (gol,'(1x,39("--"))'); call goPr write (gol,'(1x)'); call goPr write (gol,'("exitus: program has terminated normally.")'); call goPr call tstamp(kmain,itau,'exitus: final time') call cputim(cpu3) cpu3 = cpu3-cpu0 write (gol,'(a,i8,a,f16.2,a,f16.8)') & ' exitus: no of timesteps: ',nstep,' cpu (s):',cpu3, & ' cpu/step: ',cpu3/nstep; call goPr write (gol,'(1x)'); call goPr write (gol,'(1x,39("--"))'); call goPr write (gol,'(1x,39("--"))'); call goPr write (gol,'(1x)'); call goPr end if ! all Pe's from here call Par_barrier if ( okdebug ) then write (gol,'("exitus: calling stopmpi ",i4)') myid; call goPr end if call Par_barrier ! end program: ! o mpi_finalize ! o stop !call Par_StopMPI end subroutine exitus !----------------------------------------------------------------- ! this subroutine performs the initialization of a model run or ! its continuation. ! ! ! structure of control input : ! line 1 : blank or name of alternative control input file ! line 2-3: text for run header/info (2*80 char) ! line 4-... : namelist inputz ! (must be terminated by string ! '&end' in column 2-5) ! subsequent input may be read in sr trace0 (via unit kinput0) ! ! mh, 1-dec-1990 ! ! changed to v 8.5 mh, mpi HH 23-feb-1994 ! changed to zoom version MK !----------------------------------------------------------------------- subroutine start( tread1, tread2, status ) use GO , only : gol, goPr, goErr, goLabel use GO , only : TrcFile, Init, Done, ReadRc use GO , only : TDate, NewDate, IncrDate, AnyDate use GO , only : operator(+), operator(-), rTotal use GO , only : goTranslate use dims , only : im, jm, lm, dxy11, nlat180, nlon360 use dims , only : nregions, region_name use dims , only : ndyn, nchem, nsrce, nconv, nread, ndyn_max use dims , only : ndiag, ndiagp1, ndiagp2, nwrite, ninst, ndiff use dims , only : nstep0, nstep use dims , only : idatei, idatee, idatet, idate, revert use dims , only : itaui, itaue, itaut, itau use dims , only : itaur use dims , only : iyear0, julday0, icalendo use dims , only : newyr, newday, newmonth, newsrun use dims, only : newmonth, idate use dims, only : mlen, sec_day,sec_month,sec_year use dims , only : okdebug, cpu0, kmain use dims , only : istart use dims , only : czeta, czetak, limits, cdebug use chem_param , only : ntracet, ntrace_chem #ifdef with_pycasso use global_data , only : rcfile #else use global_data , only : rcfile, rcfile_runtime #endif use global_data , only : inputdir, outdir use global_data , only : region_dat use tracer_data , only : mass_dat, chem_dat use global_data , only : declare_fields use global_data , only : fcmode, tfcday0 use Meteo , only : Meteo_Setup_Mass, Meteo_Setup_Other use sources_sinks, only : trace0, trace1 #ifdef with_budgets use budget_global, only : Init_budget_global #endif use redgridZoom , only : grid_reduced, InitRedGrid use redgridZoom , only : nred, nredmax, clustsize, imred, imredj, jred use geometry , only : geomtryh, calc_dxy use io_save , only : readhdf, readhdfmmr, savehdf, readhdfmmix use user_output , only : User_output_Init, user_output_step use datetime, only : calc_sm use datetime , only : date2tau, tau2date, julday, tstamp use zoom_tools , only : coarsen_region, update_parent use advect_tools , only : advect_m, m2phlb, m2phlb1, coarsen_ambm #ifdef MPI use mpi_const , only : localComm use mpi_const , only : MPI_INTEGER, MPI_LOGICAL, MPI_CHARACTER use mpi_const , only : MY_REAL use mpi_const , only : ierr #endif use ParTools , only : myid, root, Par_Barrier use ParTools , only : lmloc use advectm_cfl , only : Init_CFL #ifndef without_advection use advectm_cfl , only : Check_CFL #endif use user_input , only : user_input_start #ifdef oasis3 use tm5_prism , only : PRISM_start_date #endif use restart , only : Restart_Read, Restart_Write use restart_remap, only : Restart_Remap_rm #ifdef with_tendencies #ifdef oasis4 use tm5_tendency_eval , only : prism_init_time #endif use tm5_tendency_eval , only : set_init_tend #endif use meteodata, only : sp_dat ! --- in/out -------------------------------------------- type(TDate), intent(out) :: tread1, tread2 integer, intent(out) :: status ! --- const ---------------------------------------- character(len=*), parameter :: rname = mname//'/start' ! --- local --------------------------------------------- integer :: n, region integer :: imr, jmr, lmr character(len=256) :: fname, fdir character(len=32 ) :: key type(TrcFile) :: rcF type(TDate) :: tdyn, tr(2) #ifdef with_pycasso character(len=32) :: stime #else character(len=8) :: ccyymmdd #endif integer :: ccyy, mm, dd ! --- begin ----------------------------------------------- call goLabel( rname ) write (gol,'(a,": begin")') rname; call goPr write (gol,'(a,": default ...")') rname; call goPr ! set default control parameters call default write (gol,'(a,": rcfile settings ...")') rname; call goPr if ( myid == root ) then ! ! read parameters from rcfile ! call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! grid redudction ? call ReadRc( rcF, 'proces.advection.reduced', grid_reduced, status ) IF_NOTOK_RETURN(status=1) ! forecast series ? call ReadRc( rcF, 'time.fc', fcmode, status, default=.false. ) IF_ERROR_RETURN(status=1) ! read forecast day 0 ? if ( fcmode ) then #ifdef with_pycasso ! read day: yyyy-mm-dd call ReadRc( rcF, 'time.fc.day0' , stime, status ) IF_NOTOK_RETURN(status=1) call goTranslate( stime, '/-', ' ', status ) IF_NOTOK_RETURN(status=1) read (stime,*,iostat=status) ccyy, mm, dd if ( status /= 0 ) then write (gol,'("reading ccyy mm dd from : ",a)') trim(stime); call goErr TRACEBACK; call goErr; status=1; return end if #else ! read day: call ReadRc( rcF, 'time.fc.day0', ccyymmdd, status ) IF_NOTOK_RETURN(status=1) ! extract year, month, day: read (ccyymmdd,'(i4.4,2i2.2)',iostat=status) ccyy, mm, dd if ( status /= 0 ) then write (gol,'("reading ccyy mm dd from : ",a)') ccyymmdd; call goErr write (gol,'("in ",a)') rname; call goErr; status=1; return end if #endif ! store day: tfcday0 = NewDate(year=ccyy,month=mm,day=dd) else ! dummy: tfcday0 = AnyDate() end if call ReadRc(rcF, 'time.ndyn_max', ndyn, status ) IF_NOTOK_RETURN(status=1) nconv = ndyn nchem = ndyn nsrce = ndyn ndyn_max = ndyn ! ensure that every 'nread' seconds is at the end of a dynamic time step: call ReadRc( rcF, 'time.ntimestep', nread, status ) IF_NOTOK_RETURN(status=1) ! how to initialize the tracer fields ? call ReadRc( rcF, 'istart', istart, status ) IF_NOTOK_RETURN(status=1) ! input files: call ReadRc( rcF, 'inputdir', inputdir, status ) IF_NOTOK_RETURN(status=1) ! print debug info ? call ReadRc( rcF, 'okdebug', okdebug, status ) IF_NOTOK_RETURN(status=1) #ifdef with_pycasso ! name of output directory: call ReadRc( rcF, 'output.dir', outdir, status ) IF_NOTOK_RETURN(status=1) #else ! default name for output directory; ! might be overwritten later on by name including time range: call ReadRc( rcF, 'outputdir', outdir, status ) IF_NOTOK_RETURN(status=1) #endif #ifdef with_pycasso #else ! close call Done( rcF, status ) IF_NOTOK_RETURN(status=1) #endif !#ifdef with_prism #ifdef oasis4 call prism_init_time #else #ifdef with_pycasso ! start time: call ReadRc( rcF, 'jobstep.timerange.start' , stime, status ) IF_NOTOK_RETURN(status=1) call goTranslate( stime, '/-:', ' ', status ) IF_NOTOK_RETURN(status=1) read (stime,*,iostat=status) idatei if ( status /= 0 ) then write (gol,'("could not read start time from : ",a)') trim(stime); call goErr TRACEBACK; status=1; return end if ! end time: call ReadRc( rcF, 'jobstep.timerange.end' , stime, status ) IF_NOTOK_RETURN(status=1) call goTranslate( stime, '/-:', ' ', status ) IF_NOTOK_RETURN(status=1) read (stime,*,iostat=status) idatee if ( status /= 0 ) then write (gol,'("could not read end time from : ",a)') trim(stime); call goErr TRACEBACK; status=1; return end if ! 'target' time? idatet = idatee ! close: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) #else ! ! read variables from timing rc file ! call Init( rcF, rcfile_runtime, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'year1' , idatei(1), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'month1', idatei(2), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'day1' , idatei(3), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'hour1' , idatei(4), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'minu1' , idatei(5), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'sec1' , idatei(6), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'year2' , idatee(1), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'month2', idatee(2), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'day2' , idatee(3), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'hour2' , idatee(4), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'minu2' , idatee(5), status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'sec2' , idatee(6), status ) IF_NOTOK_RETURN(status=1) idatet = idatee call ReadRc( rcF, 'istart', istart, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'outputdir', outdir, status ) IF_NOTOK_RETURN(status=1) call Done( rcF, status ) IF_NOTOK_RETURN(status=1) #endif #endif end if write (gol,'(a,": broadcast ...")') rname; call goPr #ifdef MPI ! broadcast namelist call MPI_BCAST(istart ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndyn ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndyn_max ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(nconv ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndiag ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(nchem ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(nsrce ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(nread ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(nwrite ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ninst ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndiff ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(icalendo ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(iyear0 ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(idatei ,6, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(idatee ,6, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(idatet ,6, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndiagp1 ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(ndiagp2 ,1, MPI_INTEGER, root ,localComm,ierr) call MPI_BCAST(czeta ,1, MY_REAL, root ,localComm,ierr) call MPI_BCAST(czetak ,1, MY_REAL, root ,localComm,ierr) call MPI_BCAST(limits ,1, MPI_LOGICAL, root ,localComm,ierr) call MPI_BCAST(cdebug ,1, MPI_LOGICAL, root ,localComm,ierr) call MPI_BCAST(okdebug ,1, MPI_LOGICAL, root ,localComm,ierr) call MPI_BCAST( inputdir, len(inputdir), MPI_CHARACTER, root, localComm, ierr ) call MPI_BCAST( outdir, len( outdir), MPI_CHARACTER, root, localComm, ierr ) call MPI_BCAST(grid_reduced, 1 , MPI_LOGICAL , root ,localComm, ierr ) #endif write (gol,'(a,": times ...")') rname; call goPr !c !--------------------------! !c ! initialize new model run ! !c !--------------------------! !c ! initialize time parameters !c ! and set calendar !c !--------------------------! ! ! itau runs from beginning of year -1 allows running from 1-1-yyyy ! cmk aug/2003 iyear0 = idatei(1)-1 julday0=julday(1,1,iyear0) if (icalendo.eq.2) then julday0=julday(1,1,iyear0) end if call date2tau(idatei,itaui) itau=itaui call tau2date(itau,idate) call declare_fields !allocate memory for this run call date2tau(idatee,itaue) call date2tau(idatet,itaut) ! set time flags newyr=.true. newday=.true. newmonth=.true. newsrun=.true. ! step counter nstep0=0 nstep=0 #ifdef with_prism #ifdef oasis3 ! store initial time for prism coupling via oasis3 ! (time is defined as seconds from begin) PRISM_start_date = idatei ! #endif #endif write (gol,'(a,": init cfl ...")') rname; call goPr ! initialise CFL: CALL init_cfl !initialise CFL preventing routine write (gol,'(a,": determine children etc ...")') rname; call goPr ! cmk: DO THIS FIRST: CHILDREN ARE NEEDED FOR REDGRID TEST! call determine_children_etc( status ) IF_NOTOK_RETURN(status=1) write (gol,'(a,": reduced grid ...")') rname; call goPr ! set up model geometry... if ( grid_reduced ) then do n = 1, nregions if ( myid == root ) then call initredgrid( n, im(n), jm(n), status ) IF_NOTOK_RETURN(status=1) end if #ifdef MPI call MPI_BCAST( status, 1, MPI_INTEGER, root, localComm, ierr ) #endif IF_NOTOK_RETURN(status=1) end do #ifdef MPI call MPI_BCAST( nred , nregions , MPI_INTEGER, root, localComm, ierr) call MPI_BCAST( clustsize, nredmax*nregions, MPI_INTEGER, root, localComm, ierr) call MPI_BCAST( imred , 180 , MPI_INTEGER, root, localComm, ierr) call MPI_BCAST( imredj , nredmax*nregions, MPI_INTEGER, root, localComm, ierr) call MPI_BCAST( jred , nredmax*nregions, MPI_INTEGER, root, localComm, ierr) #endif end if write (gol,'(a,": split order ...")') rname; call goPr call define_splitorderzoom write (gol,'(a,": blank zoom regions ...")') rname; call goPr do region = 1,nregions call blank_zoom_region(region,region_dat(region)) end do !! display zoom regions: !if ( myid == root ) call print_zoom(region_dat) write (gol,'(a,": horizontal geometry ...")') rname; call goPr ! grid cell area in 1x1 grid call calc_dxy(dxy11,nlat180) ! grid definition do region = 1, nregions call GeomtryH( region ) end do write (gol,'(a,": initial meteo ...")') rname; call goPr ! setup from start time to end of interval [k*nread,(k+1)*nread] tread1 = NewDate(time6=idate) tread2 = NewDate(time6=(/idate(1:3),0,0,0/)) + IncrDate(sec=floor(idate(4)*3600.0/nread)*nread+nread) ! n is the number of dynamic intervals within the ! time interval for which the meteo has been setup: n = ceiling( rTotal(tread2-tread1,'sec') / real(ndyn) ) ndyn = nint( rTotal(tread2-tread1,'sec') / n ) ! setup pressure and mass fields ! o do not check pressure implied by advection call Meteo_Setup_Mass( tread1, tread2, status, isfirst=.true. ) IF_NOTOK_RETURN(status=1) #ifndef without_advection ! determine dynamic timestep ndyn for this interval [tread1,tread2] ; ! the initial number of time steps n is increased until no cfl ! violations occure call Check_CFL( tread1, tread2, n, status ) IF_NOTOK_RETURN(status=1) #endif ! * setup meteo for dynamic step tdyn+[0,ndyn] ! current time (begin of dynamics step) tdyn = NewDate( time6=idate ) ! time range of dynamic step: tr(1) = tdyn tr(2) = tdyn + IncrDate( sec=ndyn ) !! convert pu/pv to am/bm/cm, eventually time interpolated !call Setup_MassFlow( tr, status ) !if (status/=0) call escape_tm('ERROR in tracer') ! setup (interpolate?) other meteo: call Meteo_Setup_Other( tr(1), tr(2), status, isfirst=.true. ) IF_NOTOK_RETURN(status=1) ! ! ** init tracer fiels **************************** ! write (gol,'(a,": init tracer fields ...")') rname; call goPr ! allocate non transported tracers: if ( newsrun ) then do region = 1, nregions if ( ntrace_chem > 0 ) then imr = im(region) ; jmr = jm(region) ; lmr = lm(region) allocate( chem_dat(region)%rm_k(imr,jmr,lmloc,ntracet+1:ntracet+ntrace_chem) ) chem_dat(region)%rm_k(:,:,:,:) = 0.0 else nullify( chem_dat(region)%rm_k ) end if end do end if ! get initial tracer fields... (m is also available...) IF(revert ==1)THEN select case (istart) case(1) ! initial tracer fields are set ! to zero do region = 1, nregions #ifdef MPI mass_dat(region)%rm_k = 0.0 #endif mass_dat(region)%rm_t = 0.0 #ifdef slopes #ifdef MPI mass_dat(region)%rxm_k=0. mass_dat(region)%rym_k=0. mass_dat(region)%rzm_k=0. #endif mass_dat(region)%rxm_t=0. mass_dat(region)%rym_t=0. mass_dat(region)%rzm_t=0. #ifdef secmom #ifdef MPI mass_dat(region)%rxxm_k=0. mass_dat(region)%rxym_k=0. mass_dat(region)%rxzm_k=0. mass_dat(region)%ryym_k=0. mass_dat(region)%ryzm_k=0. mass_dat(region)%rzzm_k=0. #endif mass_dat(region)%rxxm_t=0. mass_dat(region)%rxym_t=0. mass_dat(region)%rxzm_t=0. mass_dat(region)%ryym_t=0. mass_dat(region)%ryzm_t=0. mass_dat(region)%rzzm_t=0. #endif #endif end do ! ! initial tracer fields are calculated in trace1 ! case ( 2 ) call trace1 ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! initial tracer fields are obtained ! from a save file of a previous ! model run. m is also obtained! ! case ( 3 ) ! ! On istart == 3 (restart) we check if the savefile names were specified ! Luckily, this asks for problems, because we also use istart=3 for ! rebsubmitting. So we hack in bash scripts with an enormous concentration ! of errors. Therefore, I rebuild istart=30 agagin. ! ! ! Case 3 [pls, 17/9/2010]: ! ! Prior to pycasso : look for files into "savedir" key. This ! required a save.cmd and some post-processing. ! ! With pycasso: look into "output.dir" (note that save files are ! not link/copied/.. to the savedir during post processing anymore: ! we leave that kind of work to the user -too many scenarii-). ! Corollary: this option should NOT be set by the user, only by a ! script when resubmitting a job that do not use restart. And the ! savedir key has been made obsolete. Users, who want to start with ! save files, must use istart=30 or 31 to specify fully qualified ! save filenames. #ifdef with_pycasso ! open runtime rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! directory with save files: call ReadRc( rcF, 'prev.output.dir', fdir, status ) ! close runtime file call Done( rcF, status ) IF_NOTOK_RETURN(status=1) #else ! open rcfile: call Init( rcF, rcfile_runtime, status ) IF_NOTOK_RETURN(status=1) ! directory with save files: call ReadRc( rcF, 'savedir', fdir, status ) ! close rcfile: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) #endif ! open regular rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! loop over regions: do region = 1, nregions ! name of standard save file : save_2000101000_glb3x2.hdf write (fname,'(a,"/save_",i4.4,3i2.2,"_",a,".hdf")') & trim(fdir), idate(1:4), trim(region_name(region)) ! read from save file: call ReadHDF( region, fname ) end do ! close regular rc file call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! idem, but also able to convert grids and levels ! case ( 30 ) ! ! On istart == 30, we use save-files defined in the rc-file, with ! the following format for the key(s): ! ! start.30. : ! ! Example: ! ! start.30.glb6x4 : /data/save_files/mystuff_glb6x4.hdf ! ! open regular rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! loop over regions: do region = 1, nregions ! name of standard save file : save_2000101000_glb3x2.hdf ! try if filename for istart.30 was specified in rcfile: write (key,'("start.30.",a)') trim(region_name(region)) call ReadRc( rcF, key, fname, status, default='file_name_empty' ) IF_NOTOK_RETURN (status = 1) write (gol,*) 'Using save file names from rc-file start.30.* values'; call goPr ! read from save file: call ReadHDF( region, fname ) end do ! close regular rc file call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! idem, but also able to convert grids and levels ! case ( 31 ) call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) do region=1,nregions ! name of save file for this region: call ReadRc( rcF, 'start.31.'//trim(region_name(region)), fname, status ) IF_NOTOK_RETURN(status=1) ! read all tracers, convert grid: call ReadTracer_saved( region, fname, status ) IF_NOTOK_RETURN(status=1) ! overwrite with TM4 fields ? call ReadRc( rcF, 'start.31.'//trim(region_name(region))//'.TM4', fname, status, 'none' ) IF_ERROR_RETURN(status=1) ! key found ? then read: if ( fname /= 'none' ) then call ReadTracer_saved( region, fname, status, tm4=.true. ) IF_NOTOK_RETURN(status=1) end if end do call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! read from nc restart file: ! case ( 33 ) ! read tracer masses and tendencies; ! surface pressure is read in meteo call Restart_Read( status, tracer_mass=.true., tendencies=.true. , air_mass=.true.) IF_NOTOK_RETURN(status=1) ! NOTE: Do not update parents ! ! A restart file is written after all processes in a time step are finished, ! and therefore the parents should be updated alread. ! If the 'update_parent' routine is called again here, it causes tiny ! differences between a long run and two short ones. ! ! initial tracer fields are obtained ! from a save file of a previous ! case ( 4 ) !mk now HDF !WP! read tracer mixing ratios from saveoldfile call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) do region=1,nregions call ReadRc( rcF, 'start.4.'//trim(region_name(region)), fname, status ) IF_NOTOK_RETURN(status=1) call readhdfmmr( region, fname ) end do call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! Initial tracer fields are obtained from a mmix file. ! ! This is typically the choice for combining different versions ! or extending the number of tracers ! The compounds are searched by name. If not in the mmix file ! the field is initialized as zero ! case ( 5 ) call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) do region=1,nregions call ReadRc( rcF, 'start.5.'//trim(region_name(region)), fname, status ) IF_NOTOK_RETURN(status=1) call readhdfmmix(region,fname) end do call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do case(6) ! start up field from a customized co2 field on 1x1 degree, then coarsened !call start_co2 ! instart 61: start up field interpolated from lower vertical resolution fields ! saved in a restart file, then interpolated to higher resolution. m, sp, pblh ! are all handled by standard TM5 startup/meteo.F90. case( 61 ) call restart_remap_rm(status) IF_NOTOK_RETURN(status=1) call restart_write(status=status,key="_remapped") IF_NOTOK_RETURN(status=1) case(62) ! start up co2 fields (rm) from restart write (gol,'(a,": read tracer masses from restart file ...")') rname; call goPr call Restart_Read( status, tracer_mass=.true., tendencies=.false., & air_mass=.false., surface_pressure=.false.,pressure=.false., & override_istart=.true.) IF_NOTOK_RETURN(status=1) ! ! User defined ! case ( 9 ) call user_input_start( status ) IF_NOTOK_RETURN(status=1) ! transfer information from fine--->coarse... write (gol,'(a,": update parents ...")') rname; call goPr do region = nregions,2,-1 call update_parent( region ) end do ! ! error ... ! case default write (gol,'("unsupported istart : ",i2)') istart; call goErr call goErr; status=1; return end select end if ! forward run #ifdef with_tendencies call set_init_tend (rcfile,status) IF_NOTOK_RETURN(status=1) #endif write (gol,'(a,": cpu time counting ...")') rname; call goPr call cputim(cpu0) call Par_Barrier write (gol,'(a,": trace0 ...")') rname; call goPr call trace0( status ) IF_NOTOK_RETURN(status=1) call Par_Barrier #ifdef with_budgets write (gol,'(a,": budgets ...global:")') rname; call goPr call init_budget_global ( status ) IF_NOTOK_RETURN(status=1) #endif !>>> moved to TM5/TM_model_run !write (gol,'(a,": user output ...")') rname; call goPr ! !! initialise user-specified output (!CMK only once!) !call user_output_init( status ) !IF_NOTOK_RETURN(status=1) ! !! write output for start time: !itaur = itau ! <--- seems to be necessary for user_output ... !do region = 1, nregions ! call user_output_step( region, status ) ! IF_NOTOK_RETURN(status=1) !end do !<<< write (gol,'(a,": model start up complete")') rname; call goPr ! ok call goLabel(); status=0 end subroutine start ! *** subroutine ReadTracer_saved( region, fname, status, tm4 ) use file_hdf , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData use grid , only : TllGridInfo, TLevelInfo, Init, Done, Fill3D use dims , only : im, jm, lm use chem_param , only : ntrace, ntracet, fscale, names, ntrace_chem use partools , only : myid, root use partools , only : previous_par, tracer_active, ntracetloc, tracer_loc use partools , only : Par_Barrier #ifdef MPI use mpi_comm , only : scatter_after_read_k_3d #endif use tracer_data , only : mass_dat, chem_dat use meteo , only : lli, levi, sp_dat, m_dat ! --- in/out ------------------------------- integer, intent(in) :: region character(len=*), intent(in) :: fname integer, intent(out) :: status logical, intent(in), optional :: tm4 ! --- const ------------------------------- character(len=*), parameter :: rname = mname//'/ReadTracer_saved' ! 31 (from 60) layer defintition for TM4 save files: real, parameter :: at_tm4_60_31(0:31) = (/ & 0., 7.367743, 467.333588, 1385.912598, 2887.696533, & 4941.778320, 6144.314941, 8802.356445, 10209.500977, 11632.758789, & 14411.124023, 15706.447266, 16899.468750, 17961.357422, 18864.750000, & 19584.330078, 20097.402344, 20384.480469, 20222.205078, 19755.109375, & 19027.695313, 18045.183594, 16819.474609, 15379.805664, 13775.325195, & 12077.446289, 10376.126953, 8765.053711, 6018.019531, 3960.291504, & 2082.273926, 0./) real, parameter :: bt_tm4_60_31(0:31) = (/ & 1.00000000, 0.99401945, 0.96764523, 0.93194032, 0.87965691, & 0.81125343, 0.77159661, 0.68326861, 0.63554746, 0.58616841, & 0.48477158, 0.43396294, 0.38389215, 0.33515489, 0.28832296, & 0.24393314, 0.20247590, 0.16438432, 0.09967469, 0.07353383, & 0.05169041, 0.03412116, 0.02067788, 0.01114291, 0.00508112, & 0.00181516, 0.00046139, 0.00007582, 0.00000000, 0.00000000, & 0.00000000, 0.00000000/) ! 34 (from 91) layer defintition for TM4 save files: real, parameter :: at_tm4_91_34(0:34) = (/ & 0.000000, 6.575628, 336.772369, 1297.656128, 3010.146973, & 5422.802734, 8356.252930, 11543.166992, 14665.645508, 17385.595703, & 19348.775391, 20319.011719, 20348.916016, 19919.796875, 19184.544922, & 18191.029297, 16990.623047, 15638.053711, 14192.009766, 12713.897461, & 11262.484375, 9873.560547, 8564.624023, 7341.469727, 6199.839355, & 4663.776367, 3358.425781, 2292.155518, 1463.163940, 857.945801, & 450.685791, 204.637451, 76.167656, 21.413612, 0.000000/) real, parameter :: bt_tm4_91_34(0:34) = (/ & 1.000000, 0.994204, 0.973466, 0.935157, 0.875518, & 0.795385, 0.698224, 0.589317, 0.475016, 0.362203, & 0.259554, 0.176091, 0.112979, 0.080777, 0.055474, & 0.036227, 0.022189, 0.012508, 0.006322, 0.002765, & 0.001000, 0.000279, 0.000055, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/) ! 34 (from standard) layer defintition for TM4 save files: real, parameter :: at_tm4_std_34(0:34) = (/ & 0., 7.367743, 210.393890, 855.361755, & 2063.779785, 3850.913330, 6144.314941, 8802.356445, & 11632.758789, 14411.124023, 16899.468750, 17961.357422, & 18864.750000, 19584.330078, 20097.402344, 20384.480469, & 20429.863281, 20222.205078, 19755.109375, 19027.695313, & 18045.183594, 16819.474609, 15379.805664, 13775.325195, & 12077.446289, 10376.126953, 8765.053711, 7306.631348, & 6018.019531, 3960.291504, 1680.640259, 713.218079, & 298.495789, 95.636963, 0./) real, parameter :: bt_tm4_std_34(0:34) = (/ & 1.00000000, 0.99401945, 0.97966272, 0.95182151, & 0.90788388, 0.84737492, 0.77159661, 0.68326861, & 0.58616841, 0.48477158, 0.38389215, 0.33515489, & 0.28832296, 0.24393314, 0.20247590, 0.16438432, & 0.13002251, 0.09967469, 0.07353383, 0.05169041, & 0.03412116, 0.02067788, 0.01114291, 0.00508112, & 0.00181516, 0.00046139, 0.00007582, 0.00000000, & 0.00000000, 0.00000000, 0.00000000, 0.00000000, & 0.00000000, 0.00000000, 0.00000000/) ! --- local ------------------------------- integer :: imr, jmr, lmr real, dimension(:,:,:), pointer :: sp real, dimension(:,:,:), pointer :: m real, dimension(:,:,:,:), pointer :: rm type(THdfFile) :: hdf type(TSds) :: sds_rm, sds_rmc integer :: n, n_loc integer :: ntrace_in, ntracet_in integer :: n_in, nl character(len=8), allocatable :: names_in(:) integer :: xbeg_in, ybeg_in real :: dx_in, dy_in integer :: im_in, jm_in, lm_in real, allocatable :: at_in(:), bt_in(:) type(TllGridInfo) :: lli_in type(TLevelInfo) :: levi_in real, allocatable :: rm_in(:,:,:,:) real, allocatable :: rmt(:,:,:) integer :: l logical :: with_tm4 ! --- begin ------------------------------- ! should be parallel over tracers ... if ( previous_par(region) /= 'tracer' ) then write (gol,'("expected parallel over tracers:")'); call goErr write (gol,'(" region : ",i3)') region; call goErr write (gol,'(" previous_par : ",a)') trim(previous_par(region)); call goErr TRACEBACK; call goErr; status=1; return end if ! grid sizes imr = im(region) jmr = jm(region) lmr = lm(region) ! set pointers sp => sp_dat(region)%data m => m_dat(region)%data rm => mass_dat(region)%rm_t ! special input ? with_tm4 = .false. if ( present(tm4) ) with_tm4 = tm4 ! open hdf file (all pe) call Init( hdf, fname, 'read', status ) IF_NOTOK_RETURN(status=1) ! read list with tracer names in file: call ReadAttribute( hdf, 'ntrace', ntrace_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'ntracet', ntracet_in, status ) IF_NOTOK_RETURN(status=1) allocate( names_in(ntrace_in) ) call ReadAttribute( hdf, 'names', names_in, status ) IF_NOTOK_RETURN(status=1) ! read level definition: call ReadAttribute( hdf, 'lm', lm_in, status ) IF_NOTOK_RETURN(status=1) allocate( at_in(0:lm_in) ) allocate( bt_in(0:lm_in) ) if ( with_tm4 ) then select case ( lm_in ) case ( 31 ) at_in = at_tm4_60_31 bt_in = bt_tm4_60_31 case ( 34 ) at_in = at_tm4_std_34 bt_in = bt_tm4_std_34 case default write (gol,'("input from TM4 save file implemented for 60-31 and 91-34levels only ...")'); call goErr TRACEBACK; status=1; return end select else call ReadAttribute( hdf, 'at', at_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'bt', bt_in, status ) IF_NOTOK_RETURN(status=1) end if call Init( levi_in, lm_in, at_in, bt_in, status ) IF_NOTOK_RETURN(status=1) deallocate( at_in ) deallocate( bt_in ) ! read grid definition: call ReadAttribute( hdf, 'im', im_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'jm', jm_in, status ) IF_NOTOK_RETURN(status=1) if ( with_tm4 ) then xbeg_in = -180.0 dx_in = 360.0/im_in ybeg_in = -90.0 dy_in = 180.0/jm_in else call ReadAttribute( hdf, 'xbeg', xbeg_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'dx', dx_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'ybeg', ybeg_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'dy', dy_in, status ) IF_NOTOK_RETURN(status=1) end if call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, & ybeg_in+0.5*dy_in, dy_in, jm_in, status ) IF_NOTOK_RETURN(status=1) ! temporary storage: allocate( rm_in(im_in,jm_in,lm_in,1) ) allocate( rmt(imr,jmr,lmr) ) ! open tracer masses: call Init( sds_rm, hdf, 'rm', status ) IF_NOTOK_RETURN(status=1) if ( .not. with_tm4 ) then call Init( sds_rmc, hdf, 'rmc', status ) if(status .ne. 0) write (gol,'("WARNING - No short lived tracer array found in save file...continueing execution")'); call goPr end if ! loop over all tracers: do n = 1, ntrace ! search input record: n_in = -1 do nl = 1, ntrace_in if ( names(n) == names_in(nl) ) then n_in = nl exit end if end do ! not found ? if ( n_in < 0 ) then write (gol,'("WARNING - tracer not found in save file:")'); call goPr write (gol,'("WARNING - file : ",a)') trim(fname); call goPr write (gol,'("WARNING - tracer : ",i3," ",a)') n, trim(names(n)); call goPr cycle end if ! transported or short lived ? if ( n <= ntracet ) then ! not on this proces ? if ( .not. tracer_active(n) ) cycle ! read 3D tracer field from rm or rmc: if ( n_in <= ntracet_in ) then call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) else call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) end if ! expand polar cell for TM4 files ... if ( with_tm4 ) then do l = 1, lm_in rm_in(:, 1,l,1) = rm_in(1, 1,l,1) / im_in rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in end do end if ! local index: n_loc = tracer_loc(n) ! convert grids and levels: call Fill3D( lli(region), levi, 'n', sp(1:imr,1:jmr,1), rm(1:imr,1:jmr,1:lmr,n_loc), & lli_in, levi_in, rm_in(:,:,:,1), 'sum', status ) IF_NOTOK_RETURN(status=1) ! info write (gol,*) 'initialised tracer ', n, names(n), '; max mass, max vmr : ', & maxval(rm(1:imr,1:jmr,1:lmr,n_loc)), & maxval(rm(1:imr,1:jmr,1:lmr,n_loc)/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr else ! short lived tracer if ( myid == root ) then ! read 3D tracer field: if ( with_tm4 ) then ! tm4 save files have short lived tracers in rm call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) ! expand polar cell ... do l = 1, lm_in rm_in(:, 1,l,1) = rm_in(1, 1,l,1) / im_in rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in end do else ! read 3D tracer field from rm or rmc: if ( n_in <= ntracet_in ) then call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) else call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) end if end if ! convert grids and levels: call Fill3D( lli(region), levi, 'n', sp(1:imr,1:jmr,1), rmt, & lli_in, levi_in, rm_in(:,:,:,1), 'sum', status ) IF_NOTOK_RETURN(status=1) ! info write (gol,*) 'initialised tracer ', n, names(n), '; max mass, max vmr : ', & maxval(rmt), & maxval(rmt/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr else ! for safety .. rmt = 0.0 end if ! distribute over levels: #ifdef MPI call scatter_after_read_k_3d( rmt, imr,jmr,lmr, 0,0,0, & chem_dat(region)%rm_k(:,:,:,n), root ) #else chem_dat(region)%rm_k(:,:,:,n) = rmt #endif end if end do ! clear deallocate( names_in ) deallocate( rm_in ) deallocate( rmt ) ! done with grid definitions: call Done( lli_in, status ) IF_NOTOK_RETURN(status=1) call Done( levi_in, status ) IF_NOTOK_RETURN(status=1) ! close data sets: call Done( sds_rm, status ) IF_NOTOK_RETURN(status=1) if ( .not. with_tm4 ) then call Done( sds_rmc, status ) end if ! close hdf file (all pe) call Done( hdf, status ) IF_NOTOK_RETURN(status=1) ! wait ... call Par_Barrier ! ok status = 0 end subroutine ReadTracer_saved ! ! *********************************************************************** ! cpu time ! *********************************************************************** ! subroutine cputim(time ) implicit none real,intent(out) :: time integer :: clock,clockrate call system_clock(clock,clockrate) time=clock*(1.0/clockrate) !time=second() end subroutine cputim ! the following routines are called at start-up to calculate some ! useful information for the model, like array adresses and ! the lay-out of the defined regions. ! ! determines children and ibeg,iend,jbeg,jend,lbeg,lend ! the subroutine is normally called once per run of TM5 ! written by patrick berkvens and mike botchev, march-june 1999 ! subroutine determine_children_etc( status ) use dims, only : nregions use dims, only : parent, children use dims, only : im, jm, lm use dims, only : xbeg, ybeg, zbeg use dims, only : xend, yend, zend use dims, only : xref, yref, zref use dims, only : ibeg, jbeg, lbeg use dims, only : iend, jend, lend use dims, only : isr, jsr use dims, only : ier, jer use dims, only : dx, dy, dz use dims, only : touch_np, touch_sp use dims, only : xcyc use dims, only : okdebug ! --- in/out --------------------------------- integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/determine_children_etc' ! --- local --------------------------------- integer :: icells, jcells, lcells, jj, region integer :: xref_, yref_, zref_ ! --- begin --------------------------------- if ( okdebug ) print *,'determine_children_etc: ' children = 0 do region=2,nregions ! increase children's counter of the parent: children(parent(region),0) = children(parent(region),0) + 1 jj = children(parent(region),0) ! this child is region: children(parent(region),jj) = region end do ! for example: !children(0,0) = 1 ! means: there is one child in the region 0 !children(0,1) = 1 ! means: this child is region 1 !children(1,0) = 1 ! means: there is one child in region 1 !children(1,1) = 2 ! means: this child is region 1 !children(2,0) = 0 ! no children in region 2 do region=2,nregions ibeg(region) = nint((xbeg(region)-xbeg(parent(region))) & / (dx/xref(parent(region)))) + 1 jbeg(region) = nint((ybeg(region)-ybeg(parent(region))) & / (dy/yref(parent(region)))) + 1 lbeg(region) = nint((zbeg(region)-zbeg(parent(region))) & / (dz/zref(parent(region)))) + 1 iend(region) = nint((xend(region)-xbeg(parent(region))) & / (dx/xref(parent(region)))) jend(region) = nint((yend(region)-ybeg(parent(region))) & / (dy/yref(parent(region)))) lend(region) = nint((zend(region)-zbeg(parent(region))) & / (dz/zref(parent(region)))) ! test xref_ = xref(region)/xref(parent(region)) yref_ = yref(region)/yref(parent(region)) zref_ = zref(region)/zref(parent(region)) ! compute icells - how many cells of parent in x-direction a region occupies if (iend(region)>ibeg(region)) then icells = iend(region)-ibeg(region)+1 else icells = im(parent(region)) - ibeg(region) + 1 + iend(region) end if if (icells*xref_/=im(region)) then write (gol,*) 'check xbeg/xend, im and ibeg/iend for region ',region; call goErr write (gol,*) ' xbeg,xend,im,ibeg,iend:', & xbeg(region),xend(region),im(region),ibeg(region),iend(region); call goErr TRACEBACK; call goErr; status=1; return end if jcells = jend(region)-jbeg(region)+1 if (jcells*yref_/=jm(region)) then write (gol,*) 'check ybeg/yend, jm and jbeg/jend for region ',region; call goErr write (gol,*) ' ybeg,yend,jm,jbeg,jend:', & ybeg(region),yend(region),jm(region),jbeg(region),jend(region); call goErr TRACEBACK; call goErr; status=1; return end if lcells = lend(region)-lbeg(region)+1 if (lcells*zref_/=lm(region)) then write (gol,*) 'check zbeg/zend, lm and lbeg/lend for region ',region; call goErr write (gol,*) ' zbeg,zend,lm,lbeg,lend:', & zbeg(region),zend(region),lm(region),lbeg(region),lend(region); call goErr TRACEBACK; call goErr; status=1; return end if end do if ( okdebug ) then write (gol,*) ' ibeg(2:nregions): ',ibeg; call goPr write (gol,*) ' iend(2:nregions): ',iend; call goPr write (gol,*) ' jbeg(2:nregions): ',jbeg; call goPr write (gol,*) ' jend(2:nregions): ',jend; call goPr write (gol,*) ' lbeg(2:nregions): ',lbeg; call goPr write (gol,*) ' lend(2:nregions): ',lend; call goPr write (gol,*) ' parents: ',parent; call goPr do region=1,nregions write (gol,*) ' children(',region, & ',0:nregions)=',children(region,:); call goPr end do end if ! determine isr,ier,jsr,jer: the scope of the cells to be ! processed per region. ! processes like vertical mixing, chemistry, sources/sinks ! are processed on this scope: ! do i=isr(region),ier(region) ! do j=jsr(region),jer(region) isr(1) = 1 ier(1) = im(1) jsr(1) = 1 jer(1) = jm(1) do region = 2, nregions xref_ = xref(region)/xref(parent(region)) yref_ = yref(region)/yref(parent(region)) isr(region) = xref_ + 1 ier(region) = im(region)-xref_ jsr(region) = yref_+1 jer(region) = jm(region)-yref_ if(xcyc(region)/=0) then isr(region) = 1 ier(region) = im(region) end if if(touch_sp(region)==1) jsr(region) = 1 if(touch_np(region)==1) jer(region) = jm(region) end do ! ok status = 0 end subroutine determine_children_etc subroutine define_splitorderzoom ! ! splitorderzoom is splitorder specified for each region ! the subroutine is normally called once per run of TM5 ! ATTN: subroutine was written based on assumption that a region can not ! have number less than number of its parent: region > parent(region) ! written by patrick berkvens and mike botchev, march-june 1999 ! use dims implicit none ! local integer :: i, j, j0, region, tref_ !character,dimension(3):: reverse if ( okdebug ) then write (gol,*) 'define_splitorderzoom: '; call goPr end if splitorderzoom = ' ' splitorderzoom(1,1:nsplitsteps) = splitorder if ( zoom_mode == 1 ) then ! zoom_mode==1 means: ! x y z | z y x ! x y z z y x | z y x x y z ! ... do region=2,nregions tref_ = tref(region)/tref(parent(region)) if ( (tref_ > 4 ) .or. ( tref_ == 3 ) ) then write (gol,*) 'define_splitorderzoom: ERROR region = ',region; call goErr write (gol,*) 'define_splitorderzoom: ', 'wrong value for tref(region): ',tref(region); call goErr write (gol,*) 'define_splitorderzoom: ', 'should be 1,2 or 4 times bigger than tref of its parent'; call goErr write (gol,*) 'define_splitorderzoom: ', 'use another value for parameter zoom_mode'; call goErr stop end if j0 = 1 ! step counter for parent j = 1 ! step counter for region do while ( splitorderzoom(parent(region),j0) /= ' ' ) ! use step triple of parent tref_ times: do i=1,tref_ if (mod(i,2)==1) then ! copy step triple of parent: !WP! changed 2 into n_oper splitorderzoom(region,j:j+(n_operators-1)) = & splitorderzoom(parent(region),j0:j0+(n_operators-1)) else ! step triple of parent in reverse order: ! WP! changed 2 into n_op splitorderzoom(region,j:j+(n_operators-1)) = & reverse(splitorderzoom(parent(region), & j0:j0+(n_operators-1))) end if !WP! changed 3 into n_ope... j = j + (n_operators) end do ! next step triple of parent !WP! changed 3 into n_oper.. j0 = j0 + (n_operators) end do end do else if ( zoom_mode == 2 ) then ! zoom_mode==2 means: ! x y z | z y x ! x y z x y z | z y x z y x ! ... do region=2,nregions tref_ = tref(region)/tref(parent(region)) j0 = 1 ! step counter for parent j = 1 ! step counter for region do while (splitorderzoom(parent(region),j0)/=' ') ! replicate step triple of parent tref_ times: do i=1,tref_ splitorderzoom(region,j:j+(n_operators-1)) = & splitorderzoom(parent(region),j0:j0+(n_operators-1)) j = j + (n_operators) end do ! next step of parent j0 = j0 + (n_operators) end do end do else write (gol,*) 'define_splitorderzoom: wrong value for zoom_mode ',zoom_mode; call goErr stop end if contains function reverse(str) result(inv_str) !WP! changed to n_operators.... character,dimension(n_operators):: str,inv_str integer ::i do i=1,(n_operators) inv_str(i)=str(n_operators+1-i) end do end function reverse end subroutine define_splitorderzoom ! *** subroutine print_zoom(region_dat) use dims ! --- in/out ---------------------------------- type(region_data), intent(in) :: region_dat(nregions) ! --- local --------------------------------- integer :: region, j integer :: n ! --- begin --------------------------------- write (gol,'("print_zoom")'); call goPr write (gol,'("print_zoom: nregions : ",i6)') nregions; call goPr write (gol,'("print_zoom: xbeg : ",20i6)') xbeg; call goPr write (gol,'("print_zoom: xend : ",20i6)') xend; call goPr write (gol,'("print_zoom: ybeg : ",20i6)') ybeg; call goPr write (gol,'("print_zoom: yend : ",20i6)') yend; call goPr write (gol,'("print_zoom: zbeg : ",20i6)') zbeg; call goPr write (gol,'("print_zoom: zend : ",20i6)') zend; call goPr write (gol,'("print_zoom: xref : ",20i6)') xref; call goPr write (gol,'("print_zoom: yref : ",20i6)') yref; call goPr write (gol,'("print_zoom: zref : ",20i6)') zref; call goPr write (gol,'("print_zoom: im : ",20i6)') im; call goPr write (gol,'("print_zoom: jm : ",20i6)') jm; call goPr write (gol,'("print_zoom: lm : ",20i6)') lm; call goPr write (gol,'("print_zoom: isr : ",20i6)') isr; call goPr write (gol,'("print_zoom: ier : ",20i6)') ier; call goPr write (gol,'("print_zoom: jsr : ",20i6)') jsr; call goPr write (gol,'("print_zoom: jer : ",20i6)') jer; call goPr write (gol,'("print_zoom: ibeg(2:nregions) : ",20i6)') ibeg; call goPr write (gol,'("print_zoom: iend(2:nregions) : ",20i6)') iend; call goPr write (gol,'("print_zoom: jbeg(2:nregions) : ",20i6)') jbeg; call goPr write (gol,'("print_zoom: jend(2:nregions) : ",20i6)') jend; call goPr write (gol,'("print_zoom: lbeg(2:nregions) : ",20i6)') lbeg; call goPr write (gol,'("print_zoom: lend(2:nregions) : ",20i6)') lend; call goPr write (gol,'("print_zoom: parents : ",20i6)') parent; call goPr do region=1,nregions write (gol,'("print_zoom: children(",i2,",0:nregions)=",20i3)') region, children(region,:); call goPr end do do region=1,nregions n = count(splitorderzoom(region,:)/=' ') write (gol,'("print_zoom: splitorderzoom(",i2,",:)=",100(a1))') & region, splitorderzoom(region,1:n); call goPr end do do region = 1,nregions write (gol,'("print_zoom: Zoomed ... region ",i3)') region ; call goPr do j=jm(region),1,-1 write (gol,'("print_zoom: ",i3,a3,200i1)') j,'---',region_dat(region)%zoomed(:,j); call goPr end do write (gol,'("print_zoom: Edge ... region ",i3)') region ; call goPr do j=jm(region),1,-1 write (gol,'("print_zoom: ",i3,a3,200i1)') j,'---',region_dat(region)%edge(:,j); call goPr end do end do end subroutine print_zoom ! *** subroutine blank_zoom_region(region,region_dat) ! ! determine which cells in a region are calculated in an overlying ! zoom region. Only the core of a zoom region is flagged. The ! interface cells are not flagged, since they are treated coarse ! ref: Krol et al, IAMAS paper, 2001 ! ! Example region1: ! Zoomed I2 represent the core-zoom region ! 1111111111111111111 ! 1111111111111111111 ! 1111222111111111111 ! 1111222111111111111 ! 1111222111111111111 ! 1111222111111111111 ! 1111111111111111111 ! 1111111111111111111 ! Edge !indicates the interface cells ! 0000000000000000000 ! 0003222300000000000 ! 0001000100000000000 ! 0001000100000000000 ! 0001000100000000000 ! 0001000100000000000 ! 0003222300000000000 ! 0000000000000000000 ! 3 means x+y, 2 means y, and 1 means only x interface cell. ! use dims implicit none ! input/output integer, intent(in) :: region type(region_data),intent(inout) :: region_dat ! local integer :: i, j, n, n1, xref_, yref_ integer,dimension(:,:),pointer :: zoomed, edge zoomed => region_dat%zoomed edge => region_dat%edge zoomed = region edge = 0 do n=1,children(region,0) !loop over the children n1 = children(region,n) !child region do i=1,im(region) if( (i <= ibeg(n1) .or. i >= iend(n1)) .and. (xcyc(n1) /= 1) ) cycle do j=1,jm(region) if( (j <= jbeg(n1)) .and. (touch_sp(n1)/=1) ) cycle if( (j >= jend(n1)) .and. (touch_np(n1)/=1) ) cycle zoomed(i,j) = n1 end do end do if ( xcyc(n1) /= 1 ) then edge(ibeg(n1),jbeg(n1):jend(n1)) = 1 edge(iend(n1),jbeg(n1):jend(n1)) = 1 end if if ( touch_sp(n1) /= 1 ) then edge(ibeg(n1):iend(n1),jbeg(n1)) = & edge(ibeg(n1):iend(n1),jbeg(n1)) + 2 end if if ( touch_np(n1) /= 1 ) then edge(ibeg(n1):iend(n1),jend(n1)) = & edge(ibeg(n1):iend(n1),jend(n1)) + 2 end if end do ! put the edge of the region itself to -1..... if ( parent(region) /= 0 ) then xref_ = xref(region)/xref(parent(region)) yref_ = yref(region)/yref(parent(region)) if(xcyc(region)/=1) then edge(1:xref_,:) = -1 edge(im(region)-xref_+1:im(region),:) = -1 end if if(touch_sp(region)/=1) edge(:,1:yref_) = -1 if(touch_np(region)/=1) edge(:,jm(region)-yref_+1:jm(region)) = -1 end if nullify(zoomed) nullify(edge) end subroutine blank_zoom_region end module initexit