! #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" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: MODELINTEGRATION ! ! !DESCRIPTION: handles time ordering of model processes !\\ !\\ ! !INTERFACE: ! MODULE MODELINTEGRATION ! ! !USES: ! use GO, only : gol, goPr, goErr, goLabel IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Proces_Init, Proces_Done, Proces_Update public :: Proces_Region ! ! !PRIVATE DATA MEMBERS: ! character(len=1) :: output_after_step character(len=*), parameter :: mname = 'ModelIntegration' ! module name ! ! timer handles integer :: itim_advectx, itim_advecty, itim_advectz integer :: itim_vertical, itim_chemistry, itim_source_sink ! ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PROCES_INIT ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE PROCES_INIT( status ) ! ! !USES: ! use GO, only : GO_Timer_Def use GO, only : TrcFile, Init, Done, ReadRc use user_output, only : User_output_Init use Binas, only : p_global use dims, only : nregions, nregions_all, iglbsfc use global_data, only : rcfile use Meteodata, only : sp_dat, phlb_dat, m_dat use Meteodata, only : sp1_dat, sp2_dat, tsp_dat, spm_dat use Meteodata, only : mfu_dat, mfv_dat, mfw_dat use Meteodata, only : temper_dat, humid_dat use Meteodata, only : entu_dat, entd_dat, detu_dat, detd_dat use Meteodata, only : lwc_dat, iwc_dat, cc_dat, cco_dat, ccu_dat use Meteodata, only : gph_dat, omega_dat use Meteodata, only : oro_dat, lsmask_dat use Meteodata, only : Set use Meteo, only : Meteo_Alloc #ifndef without_advection use Advect, only : Advect_Init #endif #ifndef without_convection use Convection, only : Convection_Init #ifndef without_diffusion use Diffusion, only : Diffusion_Init #endif #endif #ifndef without_chemistry use Chemistry, only : Chemistry_Init #endif use Sources_Sinks, only : Sources_Sinks_Init #ifndef without_dry_deposition use Dry_Deposition, only : Dry_Deposition_Init #endif #ifndef without_wet_deposition use Wet_Deposition, only : Wet_Deposition_Init #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Proces_Init' type(TrcFile) :: rcF integer :: n ! --- BEGIN --------------------------------- call goLabel(rname) ! ! (1) meteo always used ! do n = 1, nregions ! current (surface) pressure and mass: call Set( sp_dat(n), status, used=.true. ) call Set( m_dat(n), status, used=.true. ) call Set( phlb_dat(n), status, used=.true. ) ! surface pressures, tendency, and average (mid) call Set( sp1_dat(n), status, used=.true. ) call Set( sp2_dat(n), status, used=.true. ) call Set( tsp_dat(n), status, used=.true. ) call Set( spm_dat(n), status, used=.true. ) end do ! ! (2) Init processes (these inits should select the appropriate meteo to use) ! #ifndef without_advection call Advect_Init( status ) IF_NOTOK_RETURN(status=1) #endif #ifndef without_convection call Convection_Init( status ) IF_NOTOK_RETURN(status=1) #ifndef without_diffusion call Diffusion_Init( status ) IF_NOTOK_RETURN(status=1) #endif #endif #ifndef without_chemistry call Chemistry_Init( status ) IF_NOTOK_RETURN(status=1) #endif call Sources_Sinks_Init( status ) IF_NOTOK_RETURN(status=1) #ifndef without_dry_deposition call Dry_Deposition_Init( status ) IF_NOTOK_RETURN(status=1) #endif #ifndef without_wet_deposition call Wet_Deposition_Init( status ) IF_NOTOK_RETURN(status=1) #endif call user_output_init( status ) IF_NOTOK_RETURN(status=1) ! ! (3) Check for (and set-to-use) required meteo ! do n = 1, nregions_all ! convec requires gph and omega: if ( entu_dat(n)%used ) then call Set( gph_dat(n), status, used=.true. ) call Set( omega_dat(n), status, used=.true. ) end if ! omega (Pa/s) is computed form mfw (kg/s): if ( omega_dat(n)%used ) then call Set( mfw_dat(n), status, used=.true. ) end if ! gph is computed from oro/sp/temper/humid if ( gph_dat(n)%used ) then call Set( oro_dat(n), status, used=.true. ) call Set( sp_dat(n), status, used=.true. ) call Set( temper_dat(n), status, used=.true. ) call Set( humid_dat(n), status, used=.true. ) end if ! sp is interpolated in time between sp1 and sp2: if ( sp_dat(n)%used ) then call Set( sp1_dat(n), status, used=.true. ) call Set( sp2_dat(n), status, used=.true. ) end if ! cloud covers and overhead/underfeet cloud covers if ( cc_dat(n)%used ) then call Set( cco_dat(n), status, used=.true. ) call Set( ccu_dat(n), status, used=.true. ) end if end do ! !>>> Please use obscure meteo fields only if you really need them ... ! ! If too many fields are used by default, even a simple meteo production run ! ! of an extra surface field will request fields that are not used anyway ... ! ! ! write (gol,'("WARNING - ")'); call goPr ! write (gol,'("WARNING - global orography not in use by default anymore;")'); call goPr ! write (gol,'("WARNING - see the comment in ",a)') rname; call goPr ! write (gol,'("WARNING - ")'); call goPr ! ! ! ! always store also the orography on 1x1: not here - set it in process init that really need it ! ! call Set( oro_dat(iglbsfc), status, used=.true. ) ! ! ! !<<< ! ! (4) Allocate used meteo ! call Meteo_Alloc( status ) IF_NOTOK_RETURN(status=1) ! ! (5) default values for surface pressure (might be required for mass fields used in init routines) ! do n = 1, nregions_all if ( sp_dat(n)%used ) sp_dat(n)%data = p_global if ( spm_dat(n)%used ) spm_dat(n)%data = p_global end do ! ! (6) When do we sample/accumulate output? ! ! - x,y,z,c,v,e : always after this process ! - a : after all steps (testing, not recommended) ! - d : the 'old' way (after all, outside do_steps) ! - Recommended value: 'v' (default) ! call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.after.step', output_after_step, status, default = 'v' ) IF_ERROR_RETURN(status=1) if (output_after_step == 'd') then write(gol, *) ' **********************************************' ; call goPr write(gol, *) ' Output will be collected in the OLD way' ; call goPr write(gol, *) ' consider to include e.g. output.after.step : v' ; call goPr write(gol, *) ' **********************************************' ; call goPr else write(gol, *) ' ****************************************' ; call goPr write(gol, *) ' Output will be collected after step : ',output_after_step ; call goPr write(gol, *) ' ****************************************' ; call goPr endif call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! ! (7) timing ! call GO_Timer_Def( itim_advectx , 'advectx' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_advecty , 'advecty' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_advectz , 'advectz' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_vertical , 'vertical' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_chemistry , 'chemistry' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_source_sink, 'source_sink' , status ) IF_NOTOK_RETURN(status=1) ! ! (8) done ! call goLabel() status = 0 END SUBROUTINE PROCES_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PROCES_DONE ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE PROCES_DONE( status ) ! ! !USES: ! #ifndef without_chemistry use Chemistry , only : Chemistry_Done #endif #ifndef without_advection use Advect , only : Advect_Done #endif #ifndef without_convection use Convection , only : Convection_Done #ifndef without_diffusion use Diffusion , only : Diffusion_Done #endif #endif use Sources_Sinks , only : Sources_Sinks_Done #ifndef without_dry_deposition use Dry_Deposition, only : Dry_Deposition_Done #endif #ifndef without_wet_deposition use Wet_Deposition, only : Wet_Deposition_Done #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Proces_Done' ! --- BEGIN --------------------------------- #ifndef without_advection call Advect_Done( status ) IF_NOTOK_RETURN(status=1) #endif #ifndef without_convection call Convection_Done( status ) IF_NOTOK_RETURN(status=1) #ifndef without_diffusion call Diffusion_Done( status ) IF_NOTOK_RETURN(status=1) #endif #endif #ifndef without_chemistry call Chemistry_Done( status ) IF_NOTOK_RETURN(status=1) #endif call Sources_Sinks_Done( status ) IF_NOTOK_RETURN(status=1) #ifndef without_dry_deposition call Dry_Deposition_Done( status ) IF_NOTOK_RETURN(status=1) #endif #ifndef without_wet_deposition call Wet_Deposition_Done( status ) IF_NOTOK_RETURN(status=1) #endif status = 0 END SUBROUTINE PROCES_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PROCES_UPDATE ! ! !DESCRIPTION: update process fields, for example because meteo has been changed !\\ !\\ ! !INTERFACE: ! SUBROUTINE PROCES_UPDATE( status ) ! ! !USES: ! use dims , only : nregions #ifndef without_diffusion use Diffusion , only : Calc_Kzz, diffusion_write, read_diffusion, write_diffusion #endif use sources_sinks , only : ss_after_read_meteo_update #ifndef without_dry_deposition use dry_deposition, only : dd_surface_fields #endif #ifndef without_wet_deposition use wet_deposition, only : calculate_lsp_scav #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Proces_Update' integer :: region ! --- BEGIN --------------------------------- #ifndef without_diffusion ! Read or recalculate diffusion coefficients if necessary ! If diffusion_write, then read existing DKG file or ! write it if it does not already exist. if (diffusion_write) then call read_diffusion( status ) IF_ERROR_RETURN(status=1) ! if error reading ! else if dkg file could not be found --> calculate and write if (status < 0) then call Calc_Kzz(status) IF_NOTOK_RETURN(status=1) call write_diffusion( status ) IF_NOTOK_RETURN(status=1) end if else call Calc_Kzz(status) IF_NOTOK_RETURN(status=1) end if #endif #ifndef without_dry_deposition ! calculate the dry_deposition fields from stuff in fsurf file (will coarsen...) call dd_surface_fields( status ) IF_NOTOK_RETURN(status=1) #endif #ifndef without_wet_deposition ! pre-calculate the lsp scavenging rloss1, rloss2, rloss3 do region = 1, nregions call calculate_lsp_scav(region) end do #endif ! update sources and sinks that depends on meteo call ss_after_read_meteo_update( status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE PROCES_UPDATE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PROCES_REGION ! ! !DESCRIPTION: for a given region, performs tref(region)/tref(parent(region)) steps. ! ! !example (cmk) ! XYZ ECCEZYX ! xyz eccezyx cezyx xyz ec ! xyzeccezyx ecxyzzyxec cezyxxyzec zyxeccexyz ! in this case the operations e and c should be executed in the ! interface region ! but not in the core of the zoom region (otherwise double counting) ! This results in the most simple algorithm, because you should leave ! the DO_STEPS routine AFTER ! (1) finishing all the steps... ! (2) XYZ ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE PROCES_REGION( region, tr, status ) ! ! !USES: ! use GO, only : TDate, NewDate use dims, only : tref, revert, ndyn, itaur, parent use dims, only : statusr => status use dims, only : n_operators use datetime, only : tau2date use user_output, only : user_output_step use advect_tools, only : m2phlb1 ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region type(TDate) :: tr(2) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/proces_region' integer :: i, tref_, dtime integer :: time6(6) type(TDate) :: tri(2) ! --- BEGIN ---------------------------------- ! determine refinement factor with respect to the parent tref_ = tref(region)/tref(parent(region)) ! main time step for all processes. note that all timesteps are default/2 dtime = revert*ndyn/(2*tref(region)) do i = 1, tref_ ! time step: itaur(region)+[0,dtime] call tau2date( itaur(region), time6 ) tri(1) = NewDate( time6=time6 ) call tau2date( itaur(region)+dtime, time6 ) tri(2) = NewDate( time6=time6 ) call do_steps( region, tri, status ) IF_NOTOK_RETURN(status=1) ! do the remaining steps if necessary if ( mod(statusr(region),n_operators) /= 0 ) then call do_steps( region, tri, status ) IF_NOTOK_RETURN(status=1) end if ! update region time..... itaur(region) = itaur(region) + dtime ! Accumulate or Sample various output (note: used to be the default output_after_step) if (output_after_step == 'd') then call user_output_step( region, status ) IF_NOTOK_RETURN(status=1) endif end do ! compute pressure grid from (changed) air mass: call m2phlb1 ! ok status = 0 END SUBROUTINE PROCES_REGION !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DO_STEPS ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE DO_STEPS( region, tr, status ) ! ! !USES: ! use GO, only : TDate use GO, only : GO_Timer_Start, GO_Timer_End use dims, only : statusr => status use dims, only : okdebug, splitorderzoom, n_operators #ifndef without_chemistry use chemistry, only : chemie #endif use sources_sinks, only : sources_sinks_apply #ifdef with_budgets use budget_global, only : budget_transportg #endif use chem_param, only : names, ntracet #ifndef without_convection use convection, only : convec #endif #ifndef without_advection use advectx, only : advectxzoom use advecty, only : advectyzoom use advectz, only : dynamw #endif #ifndef without_wet_deposition use wet_deposition, only : lspscav #endif use user_output, only : user_output_step !DBG ! DEBUG STATEMENT !DBG use ParTools, only : isRoot !DBG use tracer_data, only : Par_Check_mass, tracer_print ! debug !DBG ! use restart, only : Restart_Write ! debug ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region type(TDate) :: tr(2) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/do_steps' integer :: i123, reg, rgi character :: tobedone character(len=1) :: next_step, prev_step !cmk integer :: start_x, stop_x, start_y, stop_y, n !CMKTEMP ! --- BEGIN ----------------------------------- !example (cmk) ! XYZ EC CEX ! xyzeccezyx cexyzzyxec ! in this case the operations e and c should be executed in ! the interface region but not in the core of the zoom region ! (otherwise double counting) ! This results in the most simple algorithm, because you should leave ! the DO_STEPS routine AFTER a Z OR after finishing all steps... ! note that the parent is responsible for the interface with the children ! this means that the interface may not be 'update' in case of (e.g.) ! xyz a ! xyzaazyx ! in this case, the edges of the child are not updated with the proces a. ! this has consequences for the mmix and save output.... ! THE INTERFACE IS PART OF THE PARENT, NOT THE CHILD..... do i123=1,n_operators next_step = splitorderzoom(region,statusr(region)+1) if ( statusr(region) /= 0 ) then prev_step = splitorderzoom(region,statusr(region)) else prev_step = ' ' end if !DBG ! DEBUG STATEMENT !DBG if ( okdebug ) then !DBG if (isRoot) then !DBG ! print current step with capital letter (X,Y or Z) !DBG do reg=1,region !DBG tobedone = ' ' !DBG if ( reg == region ) tobedone = upper(next_step) !DBG write(gol,*) ' do_steps: ',reg,': ', splitorderzoom(reg,1:statusr(reg)),tobedone !DBG call goPr !DBG end do !DBG end if !DBG !DBG call tracer_print(region, "do_steps befor "//next_step, status) !DBG IF_NOTOK_RETURN(status=1) !DBG call Par_Check_mass(region, "bef_"//next_step, debug=.true.) !DBG end if SELECT CASE(next_step) CASE ( 'v' ) call GO_Timer_Start( itim_vertical, status ) IF_NOTOK_RETURN(status=1) #ifndef without_convection #ifdef with_budgets call budget_transportg(region,0,'conv',prev_step) #endif call convec(region, status ) IF_NOTOK_RETURN(status=1) #ifdef with_budgets call budget_transportg(region,1,'conv',prev_step) #endif #endif call GO_Timer_End( itim_vertical, status ) IF_NOTOK_RETURN(status=1) CASE ( 'c' ) call GO_Timer_Start( itim_chemistry, status ) IF_NOTOK_RETURN(status=1) #ifndef without_chemistry CALL CHEMIE( region, tr, status ) IF_NOTOK_RETURN(status=1) #endif call GO_Timer_End( itim_chemistry, status ) IF_NOTOK_RETURN(status=1) CASE( 's' ) call GO_Timer_Start( itim_source_sink, status ) IF_NOTOK_RETURN(status=1) call sources_sinks_apply( region, tr, status ) IF_NOTOK_RETURN(status=1) #ifndef without_wet_deposition ! wet removal by lsp (dry deposition is applied in chemistry) call lspscav( region ) #endif call GO_Timer_End( itim_source_sink, status ) IF_NOTOK_RETURN(status=1) CASE( 'x' ) call GO_Timer_Start( itim_advectx, status ) IF_NOTOK_RETURN(status=1) #ifndef without_advection #ifdef with_budgets call budget_transportg(region,0,'advx',prev_step) #endif CALL ADVECTXZOOM(region, status) IF_NOTOK_RETURN(status=1) #ifdef with_budgets call budget_transportg(region,1,'advx',prev_step) #endif #endif /* ADVECTX*/ call GO_Timer_End( itim_advectx, status ) IF_NOTOK_RETURN(status=1) CASE( 'y' ) call GO_Timer_Start( itim_advecty, status ) IF_NOTOK_RETURN(status=1) #ifndef without_advection #ifdef with_budgets call budget_transportg(region,0,'advy',prev_step) #endif CALL ADVECTYZOOM(region) #ifdef with_budgets call budget_transportg(region,1,'advy',prev_step) #endif #endif call GO_Timer_End( itim_advecty, status ) IF_NOTOK_RETURN(status=1) CASE( 'z' ) call GO_Timer_Start( itim_advectz, status ) IF_NOTOK_RETURN(status=1) #ifndef without_advection #ifdef with_budgets call budget_transportg(region,0,'advz',prev_step) #endif CALL DYNAMW #ifdef with_budgets call budget_transportg(region,1,'advz',prev_step) #endif #endif call GO_Timer_End( itim_advectz, status ) IF_NOTOK_RETURN(status=1) CASE default write(gol,*)'do_steps: strange value in splitorderzoom: ', next_step ; call goErr write(gol,*)'do_steps: (must be c, e, v, x, y or z)' ; call goErr TRACEBACK status=1; return END SELECT !DBG ! DEBUG STATEMENT !DBG if (okdebug) then !DBG call tracer_print(region, "do_steps after "//next_step, status) !DBG IF_NOTOK_RETURN(status=1) !DBG call Par_Check_mass(region, "aft_"//next_step, debug=.true.) !DBG ! call Restart_Write( status, key="aft_"//next_step ) !DBG ! IF_NOTOK_RETURN(status=1) !DBG end if statusr(region) = statusr(region)+1 ! e.g. xyzvec...........cevzyx ! exit after xyz or vec of cevzyx.... if ( mod(statusr(region),n_operators) == 0 ) then exit end if ! flexible sample scheme ('a' means after all steps), d = old 'default', ! output_after_step = v,c,e,x,y,z,d,a if ((output_after_step == next_step) .or. (output_after_step == 'a')) then call user_output_step( region, status ) IF_NOTOK_RETURN(status=1) endif !exit after 'yz' sequence to proces xyz children if ( next_step == 'z' ) then prev_step = splitorderzoom(region,statusr(region)-1) if ( prev_step == 'y' ) exit end if end do ! ok status = 0 END SUBROUTINE DO_STEPS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: UPPER ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! CHARACTER FUNCTION UPPER(xyz) ! ! !INPUT PARAMETERS: ! character(1),intent(in) :: xyz ! !EOP !------------------------------------------------------------------------ !BOC if (xyz=='x') then upper = 'X' else if (xyz=='y') then upper = 'Y' else if (xyz=='z') then upper = 'Z' else if (xyz=='c') then upper = 'C' else if (xyz=='s') then upper = 'S' else if (xyz=='v') then upper = 'V' else upper = '_' end if END FUNCTION UPPER !EOC END MODULE MODELINTEGRATION