!################################################################# ! ! tracer and air mass ! !### 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" ! #define IF_MPI_NOTOK_RETURN(action) if (status/=MPI_SUCCESS) then; TRACEBACK; action; return; end if ! !################################################################# module tracer_data use GO, only : gol, goPr, goErr use GO , only : gol, goErr, goPr use dims , only : nregions use global_types, only : chem_data #ifdef with_tendencies use GO , only : TDate use chem_param , only : ntrace use tm5_tendency, only : plc_ntr, plc_trname use tm5_tendency, only : plc_itr_nox use tm5_tendency, only : plc_npr, plc_prname use tm5_tendency, only : plc_ipr_pchem use tm5_tendency, only : plc_ipr_lchem, plc_ipr_pchem, plc_ipr_pemi , plc_ipr_conc, & plc_ipr_lddep, plc_ipr_lwdep, plc_ipr_tcnvd, plc_ipr_tchem #endif implicit none ! --- in/out ------------------------------------ private public :: mass_data, mass_dat public :: assign_massdat, free_massdat public :: chem_data, chem_dat public :: Par_Check_Domain public :: Par_Check_Mass public :: AdjustTracer public :: tracer_print #ifdef with_mpi public :: Tracer_Get_Here #endif #ifdef with_feedback public :: Tracer_Fill_Slabs #endif #ifdef with_tendencies public :: plc_ntr public :: plc_npr, plc_ipr_pchem, plc_ipr_lchem, plc_ipr_pemi , plc_ipr_conc, & plc_ipr_lddep, plc_ipr_lwdep, plc_ipr_tcnvd, plc_ipr_tchem public :: plc_kg_from_tm public :: plc_data, plc_dat public :: PLC_Init, PLC_Done, PLC_Set, PLC_Reset public :: PLC_AddValue, PLC_AddValue_t public :: plc_reset_period #endif ! --- const -------------------------------------- character(len=*), parameter :: mname = 'tracer_data' ! --- types ------------------------------------ type mass_data ! m_k air mass (kg) (halo = 2) ! p surface pressure (Pa) (halo = 2) ! phlb_k pressure at the grid boundaries (Pa) : ! l=1 means bottom lowest box (halo = 0) ! rm_k tracer masses (kg) (halo = 2) ! rxm_k tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! rym_k tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! rzm_k tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! m_t air mass (kg) (halo = 2) ! phlb_t pressure at the grid boundaries (Pa) : ! l=1 means bottom lowest box (halo = 0) ! rm_t tracer masses (kg) (halo = 2) ! rxm_t tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! rym_t tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! rzm_t tracer slopes (only transported) in kg/(halfgridsize) (halo = 1) ! real,dimension(:,:,:,:),pointer :: rm_k #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm_k real,dimension(:,:,:,:),pointer :: rym_k real,dimension(:,:,:,:),pointer :: rzm_k #ifdef secmom real,dimension(:,:,:,:),pointer :: rxxm_k real,dimension(:,:,:,:),pointer :: rxym_k real,dimension(:,:,:,:),pointer :: rxzm_k real,dimension(:,:,:,:),pointer :: ryym_k real,dimension(:,:,:,:),pointer :: ryzm_k real,dimension(:,:,:,:),pointer :: rzzm_k #endif #endif real,dimension(:,:,:,:),pointer :: rm_t #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm_t real,dimension(:,:,:,:),pointer :: rym_t real,dimension(:,:,:,:),pointer :: rzm_t #ifdef secmom real,dimension(:,:,:,:),pointer :: rxxm_t real,dimension(:,:,:,:),pointer :: rxym_t real,dimension(:,:,:,:),pointer :: rxzm_t real,dimension(:,:,:,:),pointer :: ryym_t real,dimension(:,:,:,:),pointer :: ryzm_t real,dimension(:,:,:,:),pointer :: rzzm_t #endif #endif end type mass_data #ifdef with_tendencies ! 3D production or loss rates, or concentrations type plc_data ! timing type(TDate) :: t ! slab of levels real, pointer :: rm_k(:,:,:) ! temporary storage for budgets if proces is parallel over tracers: logical :: par_t real, pointer :: rm_t(:,:,:) ! unit character(len=32) :: unit end type plc_data #endif ! --- var ------------------------------ ! ! transported tracer masses and air mass: ! type(mass_data), target :: mass_dat(nregions) ! ! non-transported tracer masses: ! o used to store emissions ? ! o allocated in trace0, deallocated in exitus ! chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,ntracet+1:ntracet+ntrace_chem) ! o read in in trace0 ! o applied in source1 ! type(chem_data), target :: chem_dat(nregions) #ifdef with_tendencies ! mapping of indices: integer :: plc_itr_to_tm(plc_ntr) integer :: plc_itr_from_tm(ntrace) ! mixing from tm tracer mass to tendency tracer mass: real :: plc_kg_from_tm(ntrace) ! concentration, production and loss rates ! note: extra proces '0' is used to store concentrations at begin of a step type(plc_data), target, save :: plc_dat(nregions,plc_ntr,0:plc_npr) ! reset period in hours: integer :: plc_reset_period ! hours #endif contains ! ====================================================================== #ifdef with_tendencies subroutine PLC_Init( rcfile, status ) use GO , only : AnyDate use GO , only : TrcFile, Init, Done, ReadRc use dims , only : nregions use dims , only : im, jm, lm use chem_param , only : ntracet use chem_param , only : iO3, iCO, iNO2, iNOx, iCH2O, iSO2 use chem_param , only : xmNOx, xmNO2 use tm5_tendency, only : plc_ipr_conc, plc_ipr_pchem, plc_ipr_lchem use tm5_tendency, only : plc_trname use partools , only : lmloc, tracer_active ! --- in/out ------------------------------- character(len=*), intent(in) :: rcfile integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_Init' ! --- local -------------------------------- type(TrcFile) :: rcF integer :: n integer :: imr, jmr, lmr integer :: itr, itr_tm integer :: ipr ! --- begin ------------------------------- ! ! settings ! call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'tend.reset.period', plc_reset_period, status ) IF_ERROR_RETURN(status=1) call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! ! indices mapping from plc to tm and vs ! ! initialise with dummy values: plc_itr_to_tm = -1 plc_itr_from_tm = -1 ! loop over plc tracers: do itr = 1, plc_ntr ! corresponding global tracer index and conversion factor: plc_kg_from_tm = 1.0 select case ( plc_trname(itr) ) case ( 'o3' ) ; itr_tm = iO3 case ( 'co' ) ; itr_tm = iCO case ( 'no2' ) ; itr_tm = iNO2 case ( 'nox' ) ; itr_tm = iNOx ; plc_kg_from_tm(itr_tm) = xmNO2/xmNOx ! kg N -> kg NO2 case ( 'ch2o' ) ; itr_tm = iCH2O case ( 'so2' ) ; itr_tm = iSO2 case default write (gol,'("unsupported plc tracer : ",a)') trim(plc_trname(itr)); call goErr TRACEBACK; status=1; return end select ! store: plc_itr_to_tm (itr ) = itr_tm plc_itr_from_tm(itr_tm) = itr end do ! plc tracers ! ! allocate ! ! loop over regions: do n = 1, nregions ! local grid sizes: imr = im(n) jmr = jm(n) lmr = lm(n) ! loop over tracers and proceses: do itr = 1, plc_ntr do ipr = 0, plc_npr ! no time stored yet: plc_dat(n,itr,ipr)%t = AnyDate() ! concentration arrays distributed over levels: allocate( plc_dat(n,itr,ipr)%rm_k(imr,jmr,lmloc) ) ! init to zero: plc_dat(n,itr,ipr)%rm_k = 0.0 ! temporary storage for some plc procs that are collected parallel over tracers ? ! only transported tracers could be parallel over tracers: if ( plc_itr_to_tm(itr) <= ntracet ) then ! only for some emissions and (wet) deposition: plc_dat(n,itr,ipr)%par_t = ( (ipr==plc_ipr_pemi) .and. (itr/=plc_itr_nox) ) .or. & (ipr==plc_ipr_lwdep) .or. & (ipr==plc_ipr_tcnvd) else plc_dat(n,itr,ipr)%par_t = .false. end if ! allocate storage ? if ( plc_dat(n,itr,ipr)%par_t ) then ! distributed to this processor ? if ( tracer_active(plc_itr_to_tm(itr)) ) then ! full 3d field: allocate( plc_dat(n,itr,ipr)%rm_t(imr,jmr,lmr) ) else ! allocate at least one value to avoid problems in debug mode: allocate( plc_dat(n,itr,ipr)%rm_t(1,1,1) ) end if else ! not used ... nullify( plc_dat(n,itr,ipr)%rm_t ) end if end do end do end do ! regions ! ! done ! ! ok status = 0 end subroutine PLC_Init ! *** subroutine PLC_Done( status ) use dims , only : nregions ! --- in/out ------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_Done' ! --- local -------------------------------- integer :: n integer :: itr, ipr ! --- begin ------------------------------- ! loop over regions: do n = 1, nregions ! loop over tracers and proceses: do itr = 1, plc_ntr do ipr = 0, plc_npr ! clear: deallocate( plc_dat(n,itr,ipr)%rm_k ) if ( associated(plc_dat(n,itr,plc_ipr_pemi)%rm_t) ) then deallocate( plc_dat(n,itr,plc_ipr_pemi)%rm_t ) end if end do end do end do ! regions ! ok status = 0 end subroutine PLC_Done ! *** ! ! call PLC_Set( 'action', itr, ipr, status [,ipr2=] [,fac=] ) ! ! action description ! ---------- --------------------------------------------------------- ! ! 'fill-tracer' copy current concentrations into specified plc process ! ! 'add-change' add difference between current tracer mass and ! ipr2 fields to specified process ! ! 'add' add ipr2 field to ipr ! ! 'add-t' distribute rm_t field over levels and add ! ! ! ipr ! ------------- ---------------------------------------------- ! 0 temporary storage for initial concentrations ! plc_ipr_conc instant concentrations ! plc_ipr_pchem chemical production ! plc_ipr_lchem chemical loss; now used for total tendency ! plc_ipr_pemi emission production ! plc_ipr_ldep deposition loss ! subroutine PLC_Set( action, itr, ipr, status, ipr2, fac, t ) use GO , only : TDate, rTotal, NewDate, AnyDate, IsAnyDate, operator(-), operator(>), wrtgol use dims , only : im, jm, lm use dims , only : idate use chem_param , only : iO3, iCO, iNO2, iNOx, iCH2O, iSO2 use chem_param , only : xmNOx, xmNO2 use partools , only : lmloc use partools , only : Par_Scatter_Over_Levels, tracer_id use tm5_tendency, only : plc_ntr, plc_trname use tm5_tendency, only : plc_ipr_conc, plc_ipr_pchem, plc_ipr_lchem use tm5_tendency, only : plc_npr, plc_prname ! --- in/out ------------------------------------ character(len=*), intent(in) :: action integer, intent(in) :: itr integer, intent(in) :: ipr integer, intent(out) :: status integer, intent(in), optional :: ipr2 real, intent(in), optional :: fac type(TDate), optional :: t ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_Set' ! --- local -------------------------------- integer :: n integer :: imr, jmr, lmr real, allocatable :: rm_k(:,:,:) character(len=32) :: rm_k_unit type(TDate) :: t_curr real :: afac ! --- begin ------------------------------- ! current time: t_curr = NewDate(time6=idate) if ( present(t) ) t_curr = t ! loop over regions: do n = 1, nregions ! local dims: imr = im(n) jmr = jm(n) lmr = lm(n) ! what to do ? select case ( action ) ! ! concentration ! case ( 'fill-tracer' ) ! temporary storage: allocate( rm_k(imr,jmr,lmloc) ) ! copy tracer fields in kg into rm_k: rm_k_unit = 'kg' call Tracer_Get_Slab( n, plc_itr_to_tm(itr), rm_k_unit, rm_k, status ) IF_NOTOK_RETURN(status=1) ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) ! fill slab: plc_dat(n,itr,ipr)%rm_k = rm_k ! store unit: plc_dat(n,itr,ipr)%unit = trim(rm_k_unit) ! store time: plc_dat(n,itr,ipr)%t = t_curr ! clear deallocate( rm_k ) ! ! add difference ! case ( 'add-change' ) ! requires second process ... if ( .not. present(ipr2) ) then write (gol,'("action `add-change` requires argument `ipr2` ...")'); call goErr TRACEBACK; status=1; return end if ! only if something has been stored ... if ( IsAnyDate(plc_dat(n,itr,ipr2)%t) ) then write (gol,'("action `add-change` requires valid time:")'); call goErr write (gol,'(" plc tracer name : ",a)') trim(plc_trname(itr)); call goErr write (gol,'(" plc proces name : ",a)') trim(plc_prname(ipr)); call goErr TRACEBACK; status=1; return end if ! add change if non-zero time step was taken since last storage: if ( t_curr > plc_dat(n,itr,ipr2)%t ) then ! temporary storage: allocate( rm_k(imr,jmr,lmloc) ) ! copy tracer fields in kg into rm_k: rm_k_unit = 'kg' call Tracer_Get_Slab( n, plc_itr_to_tm(itr), rm_k_unit, rm_k, status ) IF_NOTOK_RETURN(status=1) ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) ! add differnce between current and previous concentration: plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + & rm_k - plc_dat(n,itr,ipr2)%rm_k ! store unit: plc_dat(n,itr,ipr)%unit = trim(rm_k_unit) ! clear deallocate( rm_k ) else write (gol,'("WARNING - zero delay, thus no change in ",a," ",a," = rm - ",a)') & trim(plc_trname(itr)), trim(plc_prname(ipr)), trim(plc_prname(ipr2)); call goPr call wrtgol( 'WARNING - prev time : ', plc_dat(n,itr,ipr2)%t ); call goPr call wrtgol( 'WARNING - curr time : ', t_curr ); call goPr end if ! ! add two tendencies ! case ( 'add' ) ! requires second process ... if ( .not. present(ipr2) ) then write (gol,'("action `add` requires argument `ipr2` ...")'); call goErr TRACEBACK; status=1; return end if ! weight: afac = 1.0 if ( present(fac) ) afac = fac ! add fields: plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + & plc_dat(n,itr,ipr2)%rm_k * afac ! ! add budgets collected parallel over tracers: ! case ( 'add-t' ) ! only if rm_t is setup ... if ( .not. plc_dat(n,itr,ipr)%par_t ) cycle ! temporary slab of levels: allocate( rm_k(imr,jmr,lmloc) ) ! distribute over levels: call Par_Scatter_Over_Levels( plc_dat(n,itr,ipr)%rm_t, tracer_id(plc_itr_to_tm(itr)), rm_k, status ) IF_NOTOK_RETURN(status=1) ! add to correct tendency: plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + rm_k ! clear deallocate( rm_k ) ! ! unsupported ... ! case default write (gol,'("unsupported action : ",a)') trim(action); call goErr TRACEBACK; status=1; return end select end do ! regions ! ok status = 0 end subroutine PLC_Set ! *** subroutine PLC_Reset( itr, ipr, status ) use GO , only : TDate, NewDate use dims, only : idate ! --- in/out ------------------------------------ integer, intent(in) :: itr integer, intent(in) :: ipr integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_Reset' ! --- local -------------------------------- integer :: n integer :: itr type(TDate) :: t_curr ! --- begin ------------------------------- ! current time: t_curr = NewDate(time6=idate) ! loop over regions: do n = 1, nregions ! zero concentrations: plc_dat(n,itr,ipr)%rm_k = 0.0 if ( plc_dat(n,itr,ipr)%par_t ) plc_dat(n,itr,ipr)%rm_t = 0.0 ! dummy times: plc_dat(n,itr,ipr)%t = t_curr end do ! regions ! ok status = 0 end subroutine PLC_Reset ! *** ! ! Add value for single cell to rm_k array ! subroutine PLC_AddValue( region, ipr, i, j, k_loc, itrace, x, status ) ! --- in/out ----------------------------------- integer, intent(in) :: region integer, intent(in) :: ipr integer, intent(in) :: i, j, k_loc, itrace real, intent(in) :: x integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_AddValue' ! --- local -------------------------------- integer :: itr ! --- begin -------------------------------- ! plc tracer index: itr = plc_itr_from_tm(itrace) ! add value ? if ( itr > 0 ) then plc_dat(region,itr,ipr)%rm_k(i,j,k_loc) = plc_dat(region,itr,ipr)%rm_k(i,j,k_loc) + x end if ! ok status = 0 end subroutine PLC_AddValue ! *** ! ! Add value for single cell to rm_k array ! subroutine PLC_AddValue_t( region, ipr, i, j, k, itrace, x, status ) ! --- in/out ----------------------------------- integer, intent(in) :: region integer, intent(in) :: ipr integer, intent(in) :: i, j, k, itrace real, intent(in) :: x integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/PLC_AddValue_t' ! --- local -------------------------------- integer :: itr ! --- begin -------------------------------- ! plc tracer index: itr = plc_itr_from_tm(itrace) ! not a plc tracer ? if ( itr <= 0 ) then status=0; return end if ! add value if rm_t is setup: if ( plc_dat(region,itr,ipr)%par_t ) then plc_dat(region,itr,ipr)%rm_t(i,j,k) = plc_dat(region,itr,ipr)%rm_t(i,j,k) + x end if ! ok status = 0 end subroutine PLC_AddValue_t #endif ! ====================================================================== subroutine assign_massdat( reg, iaction ) use dims , only : nregions use dims , only : im, jm, lm use chem_param, only : ntracet use partools , only : ntracetloc, lmloc ! --- in/out ------------------------- integer, intent(in) :: reg integer, intent(in) :: iaction ! --- local ------------------------- integer :: imr, jmr, lmr, nt ! --- begin ------------------------- select case ( iaction ) !WP! swapping from tracer to k case ( 0 ) #ifdef MPI imr = im(reg) jmr = jm(reg) lmr = lmloc nt = ntracet allocate ( mass_dat(reg)%rm_k (-1:imr+2,-1:jmr+2,lmr,nt) ) #ifdef slopes allocate ( mass_dat(reg)%rxm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rym_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rzm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) #ifdef secmom allocate ( mass_dat(reg)%rxxm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rxym_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rxzm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%ryym_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%ryzm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rzzm_k(-1:imr+2,-1:jmr+2,lmr,nt) ) #endif #endif mass_dat(reg)%rm_k =0.0 #ifdef slopes mass_dat(reg)%rxm_k=0.0 mass_dat(reg)%rym_k=0.0 mass_dat(reg)%rzm_k=0.0 #ifdef secmom mass_dat(reg)%rxxm_k=0.0 mass_dat(reg)%rxym_k=0.0 mass_dat(reg)%rxzm_k=0.0 mass_dat(reg)%ryym_k=0.0 mass_dat(reg)%ryzm_k=0.0 mass_dat(reg)%rzzm_k=0.0 #endif #endif #endif ! swapping from k to tracer case ( 1 ) imr = im(reg) jmr = jm(reg) lmr = lm(reg) nt = ntracetloc allocate ( mass_dat(reg)%rm_t (-1:imr+2,-1:jmr+2,lmr,nt) ) #ifdef slopes allocate ( mass_dat(reg)%rxm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rym_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rzm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) #ifdef secmom allocate ( mass_dat(reg)%rxxm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rxym_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rxzm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%ryym_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%ryzm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) allocate ( mass_dat(reg)%rzzm_t(-1:imr+2,-1:jmr+2,lmr,nt) ) #endif #endif mass_dat(reg)%rm_t=0.0 #ifdef slopes mass_dat(reg)%rxm_t=0.0 mass_dat(reg)%rym_t=0.0 mass_dat(reg)%rzm_t=0.0 #ifdef secmom mass_dat(reg)%rxxm_t=0.0 mass_dat(reg)%rxym_t=0.0 mass_dat(reg)%rxzm_t=0.0 mass_dat(reg)%ryym_t=0.0 mass_dat(reg)%ryzm_t=0.0 mass_dat(reg)%rzzm_t=0.0 #endif #endif end select end subroutine assign_massdat ! *** subroutine free_massdat( reg, iaction ) use dims , only : nregions ! --- in/out ------------------------- integer, intent(in) :: reg integer, intent(in) :: iaction ! --- local ------------------------- integer :: imr, jmr, lmr, nt ! --- begin ------------------------- select case ( iaction ) ! swapping from tracer to k case ( 0 ) deallocate ( mass_dat(reg)%rm_t) #ifdef slopes deallocate ( mass_dat(reg)%rxm_t) deallocate ( mass_dat(reg)%rym_t) deallocate ( mass_dat(reg)%rzm_t) #ifdef secmom deallocate ( mass_dat(reg)%rxxm_t) deallocate ( mass_dat(reg)%rxym_t) deallocate ( mass_dat(reg)%rxzm_t) deallocate ( mass_dat(reg)%ryym_t) deallocate ( mass_dat(reg)%ryzm_t) deallocate ( mass_dat(reg)%rzzm_t) #endif #endif !WP! swapping from k to tracer case ( 1 ) #ifdef MPI deallocate ( mass_dat(reg)%rm_k) #ifdef slopes deallocate ( mass_dat(reg)%rxm_k) deallocate ( mass_dat(reg)%rym_k) deallocate ( mass_dat(reg)%rzm_k) #ifdef secmom deallocate ( mass_dat(reg)%rxxm_k) deallocate ( mass_dat(reg)%rxym_k) deallocate ( mass_dat(reg)%rxzm_k) deallocate ( mass_dat(reg)%ryym_k) deallocate ( mass_dat(reg)%ryzm_k) deallocate ( mass_dat(reg)%rzzm_k) #endif #endif #endif end select end subroutine free_massdat ! =================================================================== subroutine Par_Check_Domain( region1, others, to_check ) use dims , only : parent,children, okdebug use ParTools , only : myid, root use ParTools , only : previous_par use ParTools , only : allocate_mass #ifdef MPI use mpi_comm , only : barrier use mpi_comm , only : escape_tm_mpi #endif ! --- in/out ----------------------------------- integer,intent(in) :: region1 character(len=*),intent(in) :: others character(len=*),intent(in) :: to_check ! --- local ------------------------------------ #ifdef MPI integer :: reg,regstart,regstop,iaction integer :: nchild,regiop #endif ! --- begin ------------------------------------ #ifdef MPI select case(others) case('p') !WP! also check parents regiop=parent(region1) if ( regiop == 0 ) then regstart=region1 regstop=region1 else if ( regiop /= 0 ) then regstart=regiop regstop=region1 end if case('c') !WP! also check child nchild=children(region1,0) if ( nchild == 0 ) then regstart=region1 regstop=region1 else if ( nchild /= 0 ) then regstart=region1 !WP! take number of last child as regstop regstop=children(region1,nchild) end if case('n') !WP! check no parents or children regstart=region1 regstop=region1 case default call escape_tm_mpi('check_domain: illegal value of *others* ') end select select case(to_check) case('levels') iaction=0 case('tracer') iaction=1 case default call escape_tm_mpi('check_domain: illegal value of *to_check* ') end select do reg=regstart,regstop ! make sure data is in requested domain if ( previous_par(reg) /= to_check ) then if ( okdebug ) then write (gol,'("check_domain: data will be swapped")'); call goPr end if previous_par(reg)=to_check ! reset previous par ! assign new memory if ( allocate_mass ) call assign_massdat(reg,iaction) call barrier ! swap data for this region call swap_all_mass(reg,iaction) call barrier ! free previously allocated memory if ( allocate_mass ) call free_massdat(reg,iaction) else if ( okdebug ) then write (gol,'("check_domain: data is on proper communicator, nothing to do")'); call goPr end if end if end do #endif end subroutine Par_Check_Domain ! *** #ifdef MPI !WP! subroutine swaps the arrays m,rm,rxm,rym,and rzm !WP! from levels to tracers (iaction=1) !WP! or from levels to tracers (iaction=0) !WP! for requested regions (region1=>region2) !WP! after swapping, old data is set to zero to avoid mistakes !WP! with 'old' data !WP! 23 december 2002 subroutine swap_all_mass(region,iaction) use dims, only : okdebug, im,jm ,lm use chem_param, only : ntracet use mpi_const , only : myid, root use mpi_const , only : ntracetloc, ntracet_ar use mpi_comm , only : barrier use mpi_comm , only : scatter_after_read_k_3d, gather_tracer_k_3d use mpi_comm , only : tracer_to_k, k_to_tracer ! --- in/out ---------------------------- integer,intent(in) :: iaction integer,intent(in) :: region ! --- local ----------------------------- integer :: imr,jmr,lmr real,dimension(:,:,:,:),pointer :: rxmsource real,dimension(:,:,:,:),pointer :: rxmtarget integer :: nsend,xx, xx_end ! --- begin ----------------------------- imr=im(region) jmr=jm(region) lmr=lm(region) #ifdef slopes #ifndef secmom xx_end = 4 #else xx_end = 10 #endif #endif !WP! once for each x,y, and z do xx = 1, xx_end select case(xx) case( 1 ) rxmsource => mass_dat(region)%rm_t rxmtarget => mass_dat(region)%rm_k #ifdef slopes case( 2 ) rxmsource => mass_dat(region)%rxm_t rxmtarget => mass_dat(region)%rxm_k case( 3 ) rxmsource => mass_dat(region)%rym_t rxmtarget => mass_dat(region)%rym_k case( 4 ) rxmsource => mass_dat(region)%rzm_t rxmtarget => mass_dat(region)%rzm_k #ifdef secmom case( 5 ) rxmsource => mass_dat(region)%rxxm_t rxmtarget => mass_dat(region)%rxxm_k case( 6 ) rxmsource => mass_dat(region)%rxym_t rxmtarget => mass_dat(region)%rxym_k case( 7 ) rxmsource => mass_dat(region)%ryym_t rxmtarget => mass_dat(region)%ryym_k case( 8 ) rxmsource => mass_dat(region)%rxzm_t rxmtarget => mass_dat(region)%rxzm_k case( 9 ) rxmsource => mass_dat(region)%ryzm_t rxmtarget => mass_dat(region)%ryzm_k case( 10 ) rxmsource => mass_dat(region)%rzzm_t rxmtarget => mass_dat(region)%rzzm_k #endif #endif end select if ( iaction == 0 ) then call tracer_to_k( rxmtarget, rxmsource, & imr,jmr,lmr, 2,2,0, ntracet,ntracetloc, & ntracet_ar ) call barrier else if ( iaction == 1 ) then call k_to_tracer( rxmtarget, rxmsource, & imr,jmr,lmr, 2,2,0, ntracet, ntracetloc, & ntracet_ar) end if nullify(rxmsource) nullify(rxmtarget) if ( okdebug ) then write (gol,'("swap_all_mass: Scattered tracer mass and slopes ",i4)') xx; call goPr end if end do end subroutine swap_all_mass #endif ! ================================================================ ! PLS: Brought back to life to print min, max and total sum of air and tracer ! masses. Used for comparison with TM6. subroutine Par_Check_Mass( region, text, tr_id ) use dims , only : lm, im, jm #ifdef MPI use chem_param, only : ntracet use mpi_const , only : myid, root_t, root_k use mpi_const , only : ntracetloc, lmloc use mpi_const, only : com_trac use mpi_const, only : ierr,my_real,mpi_max,mpi_min, mpi_sum, ntracet_ar use mpi_comm, only : barrier_t use mpi_comm , only : dump_field3d, dump_field4d, gather_tracer_t ! use mpi_comm , only : summasr #endif use meteodata, only : m_dat use partools, only : which_par, previous_par use partools, only : tracer_active, tracer_loc ! --- in/out ----------------------------------- integer, intent(in) :: region character(len=*), intent(in) :: text integer, intent(in), optional :: tr_id ! default to 1 ! --- local ------------------------------------ real,dimension(:,:,:),pointer :: m real,dimension(:,:,:,:),pointer :: rm, rmt integer :: imr, jmr, lmr, trid, ijmax3(3), ijmin3(3), ijmax4(4), ijmin4(4) real :: min_one, max_one, min_all, max_all, tot_one, tot_all ! --- begin ------------------------------------ imr=im(region);jmr=jm(region);lmr=lm(region) trid=1 if(present(tr_id)) trid=tr_id #ifdef MPI which_par=previous_par(region) ! if(which_par == 'tracer'.and.ntracetloc == 0) return if(which_par == 'levels') then !.and.lmloc == 0) return print*, "check_mass not implemented for decomposition across levels" return end if #endif rm => mass_dat(region)%rm_t m => m_dat(region)%data max_one=maxval(m(1:imr,1:jmr,1:lmr)) min_one=minval(m(1:imr,1:jmr,1:lmr)) tot_one=sum(m(1:imr,1:jmr,1:lmr)) ijmax3=maxloc(m(1:imr,1:jmr,1:lmr) ) ijmin3=minloc(m(1:imr,1:jmr,1:lmr) ) #ifdef MPI ! call barrier_t ! call mpi_reduce(max_one, max_all, 1, my_real, mpi_max, root_t, com_trac, ierr) ! call mpi_reduce(min_one, min_all, 1, my_real, mpi_min, root_t, com_trac, ierr) ! call mpi_reduce(tot_one, tot_all, 1, my_real, mpi_max, root_t, com_trac, ierr) !same on all PEs, one is enough! if ( myid == root_t ) then print *, text//' [par_check_mass] AIR maxloc =',ijmax3 print *, text//' [par_check_mass] AIR minloc =',ijmin3 ! print *, text//' [par_check_mass] AIR min, max, sum: ', min_all, max_all, tot_all print *, text//' [par_check_mass] AIR min, max, sum: ', min_one, max_one, tot_one end if allocate(rmt(-1:imr+2,-1:jmr+2,1:lmr,ntracet)) call gather_tracer_t( rmt, im(region),jm(region),lm(region), 2,2,0, ntracet, rm, .false. ) if ( myid == root_t ) then max_all=maxval( rmt(1:imr,1:jmr,1:lmr,:) ) min_all=minval( rmt(1:imr,1:jmr,1:lmr,:) ) tot_all=sum( rmt(1:imr,1:jmr,1:lmr,:) ) ijmax4=maxloc( rmt(1:imr,1:jmr,1:lmr,:) ) ijmin4=minloc( rmt(1:imr,1:jmr,1:lmr,:) ) end if deallocate(rmt) if ( myid == root_t ) then print *, text//' [par_check_mass] RM maxloc =',ijmax4 print *, text//' [par_check_mass] RM minloc =',ijmin4 print *, text//' [par_check_mass] RM min, max, sum: ', min_all, max_all, tot_all end if #else print *,text//' [par_check_mass] AIR min, max, sum: ',min_one, max_one, tot_one max_one=maxval(rm(1:imr,1:jmr,1:lmr,:)) min_one=minval(rm(1:imr,1:jmr,1:lmr,:)) tot_one=sum(rm(1:imr,1:jmr,1:lmr,:)) print *,text//' [par_check_mass] RM min, max, sum: ',min_one, max_one, tot_one #endif ! for one tracer only if (tracer_active(trid)) then max_one=maxval( rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) ) min_one=minval( rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) ) tot_one=sum( rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) ) ! THIS IS OK ijmax3=maxloc(rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) ) ijmin3=minloc(rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) ) print *, text//' [par_check_mass] ID maxloc =',ijmax3 print *, text//' [par_check_mass] ID minloc =',ijmin3 print *, text//' [par_check_mass] ID, min, max, sum: ', trid, min_one, max_one, tot_one end if ! PREVIOUSLY, but was commented anyway: !#ifdef MPI ! ! call barrier ! if ( previous_par(region) == 'tracer' ) then ! if ( myid == root_t ) print*,'check_mass: checking mass in ',text ! call barrier ! m => mass_dat(region)%m_t ! rm => mass_dat(region)%rm_t ! lmr=lm(region) ! ! call dump_field3d(region,m,lmr,(/2,2,0/),'m','mtrac.hdf') ! call dump_field4d(region,rm,lmr,ntracetloc,(/2,2,0/),'rm','rmtrac.hdf') ! call summasr( region, m, rm ) ! else if ( previous_par(region) == 'levels' ) then ! if ( myid == root_k ) print*,'check_mass: checking mass in ',text ! call barrier ! m => mass_dat(region)%m_k ! rm => mass_dat(region)%rm_k ! lmr=lmloc ! ! call dump_field3d(region,m,lmr,(/2,2,0/),'m','mlev.hdf') ! call dump_field4d(region,rm,lmr,ntracet,(/2,2,0/),'rm','rmlev.hdf') ! call summasr( region, m, rm ) ! end if ! call barrier ! nullify(m) ! nullify(rm) ! !#endif ! end subroutine Par_Check_Mass ! ================================================================ ! ! Tracer slopes (kg/half cell) and second moments (kg/(half cell)^2) ! are adjusted with the same factor as applied to the total mass ! (background: if the concentration 'in the middle' is decreased ! by a factor 0.5, then the concentrations at the boundaries are ! also changed with a factor 0.5 : ! ! ^ | | | | ! tracer| 6 | o | | ! mass | | - | | | ! (kg) 4 | o | -> | | ! | - | 3 | _ o ! 2 o | 2 | _ o | ! | | 1 o | ! --------------- --------------- ! ! NOTE: not implemented for parallel yet ... ! ? should l be a local or a global level ? ! ? assume that tracer k is on this pe ? ! subroutine AdjustTracer( drm, region, i, j, l, k, status ) use GO , only : gol, goErr use ParTools, only : which_par, tracer_active, lmloc, offsetn ! --- in/out ----------------------------- real, intent(in) :: drm ! tracer mass change (kg) integer, intent(in) :: region integer, intent(in) :: i, j ! cell indices integer, intent(in) :: l ! local level integer, intent(in) :: k ! tracer index integer, intent(out) :: status ! --- const ------------------------------- character(len=*), parameter :: rname = mname//'/AdjustTracer' ! --- local ------------------------------- real, dimension(:,:,:,:), pointer :: rm #ifdef slopes real, dimension(:,:,:,:), pointer :: rxm, rym, rzm #ifdef secmom real, dimension(:,:,:,:), pointer :: rxxm, rxym, rxzm, ryym, ryzm, rzzm #endif #endif real :: rm_old real :: skf ! --- begin ----------------------------- ! set pointers select case ( which_par ) case ( 'levels' ) ! check local level index: if ( l > lmloc ) then write (gol,'("found strange local level:")'); call goErr write (gol,'(" l, lmloc : ", 2i4)') l, lmloc; call goErr TRACEBACK; status=1; return end if ! set pointers: rm => mass_dat(region)%rm_k #ifdef slopes rxm => mass_dat(region)%rxm_k rym => mass_dat(region)%rym_k rzm => mass_dat(region)%rzm_k #ifdef secmom rxxm => mass_dat(region)%rxxm_k rxym => mass_dat(region)%rxym_k rxzm => mass_dat(region)%rxzm_k ryym => mass_dat(region)%ryym_k ryzm => mass_dat(region)%ryzm_k rzzm => mass_dat(region)%rzzm_k #endif #endif case ( 'tracer' ) ! leave if this pe does not contain this tracer : if ( .not. tracer_active(k+offsetn) ) then status=0; return end if ! set pointers: rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #ifdef secmom rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t #endif #endif case default write (gol,'("unsupported which_par : ",a)') trim(which_par); call goErr TRACEBACK; status=1; return end select ! store old tracer mass: rm_old = rm(i,j,l,k) ! store new tracer mass: rm(i,j,l,k) = rm_old + drm #ifdef slopes ! adjust concentration gradients ! scale factor to be applied; trap devision by zero: if ( (rm(i,j,l,k) > tiny(1.0)) .and. (rm_old > tiny(0.0)) ) then skf = rm(i,j,l,k) / rm_old else skf = 0.0 end if ! adjust slopes: rxm(i,j,l,k) = rxm(i,j,l,k) * skf rym(i,j,l,k) = rym(i,j,l,k) * skf rzm(i,j,l,k) = rzm(i,j,l,k) * skf #ifdef secmom ! adjust second moments: rxxm(i,j,l,k) = rxxm(i,j,l,k) * skf rxym(i,j,l,k) = rxym(i,j,l,k) * skf rxzm(i,j,l,k) = rxzm(i,j,l,k) * skf ryym(i,j,l,k) = ryym(i,j,l,k) * skf ryzm(i,j,l,k) = ryzm(i,j,l,k) * skf rzzm(i,j,l,k) = rzzm(i,j,l,k) * skf #endif #endif ! ok status = 0 end subroutine AdjustTracer ! *** subroutine tracer_print( key, region, status ) use GO use dims, only : im, jm, lm use chem_param, only : ntrace, ntracet, names !use chem_param, only : io3, ico, ino use partools, only : myid, root, npes use partools, only : lmar, lmloc, offsetl use partools, only : tracer_active, tracer_loc use partools, only : previous_par, which_par use partools, only : par_barrier use meteodata, only : m_dat, phlb_dat ! --- in/out --------------------------------- character(len=*), intent(in) :: key integer, intent(in) :: region integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/tracer_print' ! --- local --------------------------------- integer :: i, j, ii, jj, ll, kk integer :: l, l_loc integer :: imr, jmr, lmr real, allocatable :: pat(:,:) #ifdef slopes real, allocatable :: patx(:,:) real, allocatable :: paty(:,:) real, allocatable :: patz(:,:) #endif real, allocatable :: patm(:,:) real, allocatable :: patp(:,:) integer :: itr logical :: filled character(len=20) :: pkey ! --- begin ---------------------------------- ! set those for the box/tracer used to debug ii=45 jj=41 ll=33 kk=18 !if (region/=3) return imr = im(region) jmr = jm(region) lmr = lm(region) allocate( pat(imr,jmr) ) #ifdef slopes allocate( patx(imr,jmr) ) allocate( paty(imr,jmr) ) allocate( patz(imr,jmr) ) #endif allocate( patm(imr,jmr) ) allocate( patp(imr,jmr) ) !if (tracer_active(1)) then !write (gol,'(a," region, previous_par : ",i1," ",a)') key, region, previous_par(region); call goBug !endif ! call par_barrier ! info ... !write (gol,*) key, region; call goBug do itr = 1, ntrace if (itr/=kk) cycle ! only requested tracer do l = 1, lmr if (l/=ll) cycle ! only requested level filled = .false. if ( itr <= ntracet ) then select case ( previous_par(region) ) case ( 'levels' ) l_loc = l - offsetl if ( (l_loc >= 1) .and. (l_loc <= lmloc) ) then write (pkey,'("[rm_k",2i3,"]")') l_loc, itr pat = mass_dat(region)%rm_k (1:imr,1:jmr,l_loc,itr) #ifdef slopes patx = mass_dat(region)%rxm_k(1:imr,1:jmr,l_loc,itr) paty = mass_dat(region)%rym_k(1:imr,1:jmr,l_loc,itr) patz = mass_dat(region)%rzm_k(1:imr,1:jmr,l_loc,itr) #endif filled = .true. end if case ( 'tracer' ) if ( tracer_active(itr) ) then write (pkey,'("[rm_t",2i3,"]")') l, tracer_loc(itr) pat = mass_dat(region)%rm_t (1:imr,1:jmr,l,tracer_loc(itr)) #ifdef slopes patx = mass_dat(region)%rxm_t(1:imr,1:jmr,l,tracer_loc(itr)) paty = mass_dat(region)%rym_t(1:imr,1:jmr,l,tracer_loc(itr)) patz = mass_dat(region)%rzm_t(1:imr,1:jmr,l,tracer_loc(itr)) #endif filled = .true. end if case default write (gol,'("unsupported which_par : ",a)') trim(which_par); call goErr TRACEBACK; status=1; return end select else if ( itr <= ntrace ) then l_loc = l - offsetl if ( (l_loc >= 1) .and. (l_loc <= lmloc) ) then write (pkey,'("[rmck",2i3,"]")') l_loc, itr pat = chem_dat(region)%rm_k(1:imr,1:jmr,l_loc,itr) ! ntracet+1:ntrace #ifdef slopes patx = 0.0 paty = 0.0 patz = 0.0 #endif filled = .true. end if end if patm = m_dat(region)%data(1:imr,1:jmr,l) patp = phlb_dat(region)%data(1:imr,1:jmr,l) if ( filled ) then !write (gol,'(a,i3," ",a,i3,3es14.4)') key, itr, names(itr), l, maxval(pat), sum(pat)/size(pat), pat(3,3); call goBug !#ifdef slopes ! write (gol,*) key, region, itr, names(itr), l, & ! sum(abs(patm)), sum(abs(patp)), & ! sum(abs(pat)), sum(abs(patx)), sum(abs(paty)), sum(abs(patz)); call goBug !#else ! write (gol,*) key, region, itr, names(itr), l, & ! sum(abs(patm)), sum(abs(patp)), & ! sum(abs(pat)); call goBug !#endif ! write (gol,*) key, trim(pkey), region, itr, names(itr), l, sum(abs(pat)); call goBug ! ! write (gol,*) key, region, itr, names(itr), l, pat(22+1,28+1), sum(abs(pat)); call goBug #ifdef slopes write (gol,*) key//" "//names(itr), pat(ii,jj), patx(ii,jj), paty(ii,jj), patz(ii,jj), patm(ii,jj) ; call goBug #else write (gol,*) key//" "//names(itr), pat(ii,jj), patm(ii,jj) ; call goBug #endif ! do j = 1, jmr ! do i = 1, imr ! if ((region==1).and.((i/=28).or.(j/=26))) cycle ! if ((region==2).and.((i/=25).or.(j/=8))) cycle ! if ((region==3).and.((i>12).or.(j>1))) cycle ! write (gol,*) key, region, itr, names(itr), i,j,l, & ! patm(i,j), patp(i,j), pat(i,j), patx(i,j), paty(i,j), patz(i,j); call goBug ! end do ! end do end if ! do j = 1, jmr ! do i = 1, imr ! if ((region==1).and.((i/=28).or.(j/=26))) cycle ! if ((region==2).and.((i/=25).or.(j/=8))) cycle ! !if ((region==3).and.((i>12).or.(j>1))) cycle ! if ((region==3).and.((i/=7).or.(j/=1))) cycle ! if (myid==root) then ! write (gol,*) key, region, i,j,l, patm(i,j), patp(i,j); call goBug ! endif ! call par_barrier ! if (myid==npes-1) then ! write (gol,*) key, region, i,j,l, patm(i,j), patp(i,j); call goBug ! endif ! call par_barrier ! end do ! end do ! call par_barrier end do ! levels ! call par_barrier end do ! tracers ! info ... ! write (gol,*) key, region; call goBug deallocate( pat ) #ifdef slopes deallocate( patx ) deallocate( paty ) deallocate( patz ) #endif deallocate( patm ) deallocate( patp ) ! ok status = 0 end subroutine tracer_print ! *** ! ! Return slab of tracer mass for levels on this processor. ! #ifdef with_tendencies subroutine Tracer_Get_Slab( region, itr, unit, rm_k, status ) use dims , only : im, jm, lm use chem_param , only : ntrace, ntracet, fscale use partools , only : myid, root use partools , only : previous_par, tracer_loc, tracer_id, tracer_active use partools , only : lmloc, lmar, offsetl use partools , only : Par_Scatter_Over_Levels use meteodata , only : m_dat ! --- in/out ------------------------------------- integer, intent(in) :: region integer, intent(in) :: itr character(len=*), intent(in) :: unit real, intent(out) :: rm_k(:,:,:) integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Tracer_Get_Slab' ! --- local --------------------------------- integer :: imr, jmr, lmr integer :: itr_loc ! --- begin ---------------------------------- ! local grid size imr = im(region) jmr = jm(region) lmr = lm(region) ! check input: if ( any(shape(rm_k)/=(/imr,jmr,lmloc/)) ) then write (gol,'("shape of output grid not ok:")'); call goErr write (gol,'(" shape(rm_k) : ",3i4)') shape(rm_k); call goErr write (gol,'(" imr,jmr,lmloc : ",3i4)') imr,jmr,lmloc; call goErr TRACEBACK; status=1; return end if ! transported or short lived tracer ? if ( (itr >= 1) .and. (itr <= ntracet) ) then ! how distributed ? select case ( previous_par(region) ) case ( 'tracer' ) ! local tracer index if requested tracer is on this pe, ! or just the first local index on other pe's to avoid array bound problems: if ( tracer_active(itr) ) then itr_loc = tracer_loc(itr) else itr_loc = 1 end if ! distribute full 3D field stored in mass_dat%rm_t on processor tracer_id(itr) ! over rm_k arrays available on each processor; ! distribution is over the last dimension (levels): call Par_Scatter_Over_Levels( & mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc), & tracer_id(itr), rm_k, status ) IF_NOTOK_RETURN(status=1) ! change unit ? select case ( unit ) case ( 'kg' ) ! rm already in kg ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) case ( 'mass-mixing-ratio' ) ! mass mixing ratio ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr) case default write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr TRACEBACK; status=1; return end select case ( 'levels' ) ! array already distributed over levels; just copy: rm_k = mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) ! change unit ? select case ( unit ) case ( 'kg' ) ! rm already in kg ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) case ( 'mass-mixing-ratio' ) ! mass mixing ratio ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr) case default write (gol,'("unsupported unit for par level : ",a)') unit; call goErr TRACEBACK; status=1; return end select case default write (gol,'("unsupported distribution:")'); call goPr write (gol,'(" region : ",i2)') region ;call goPr write (gol,'(" previous_par : ",a)') previous_par(region) ;call goPr TRACEBACK; status=1; return end select else if ( (itr > ntracet) .and. (itr <= ntrace) ) then ! extract short-lived tracer: rm_k = chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) ! change unit ? select case ( unit ) case ( 'kg' ) ! rm already in kg ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) case ( 'mass-mixing-ratio' ) ! mass mixing ratio ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale ! apply unit conversion from (kg tm tracer) to (kg plc tracer) rm_k = rm_k * plc_kg_from_tm(itr) rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr) case default write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr TRACEBACK; status=1; return end select else write (gol,'("unsupported tracer index : ",i4)') itr; call goErr write (gol,'(" ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr TRACEBACK; status=1; return end if ! transported or chem-only ! ok status = 0 end subroutine Tracer_Get_Slab #endif ! *** ! ! Fill slabs of tracer mass (each processor some levels) with given unit ! into tracer mass array (in kg); adust slopes too. ! #ifdef with_feedback subroutine Tracer_Fill_Slabs( region, itr, unit, rm_k, status ) use dims , only : im, jm, lm use chem_param , only : ntrace, ntracet, fscale use partools , only : myid, root use partools , only : previous_par, tracer_loc, tracer_id, tracer_active use partools , only : lmloc, lmar, offsetl use partools , only : Par_Gather_From_Levels use meteodata , only : m_dat ! --- in/out ------------------------------------- integer, intent(in) :: region integer, intent(in) :: itr character(len=*), intent(in) :: unit real, intent(in) :: rm_k(:,:,:) integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Tracer_Fill_Slabs' ! --- local --------------------------------- integer :: imr, jmr, lmr integer :: i, j, l integer :: itr_loc real, allocatable :: x_k(:,:,:) real, allocatable :: x_t(:,:,:) ! --- begin ---------------------------------- ! local grid size imr = im(region) jmr = jm(region) lmr = lm(region) ! check input: if ( any(shape(rm_k)/=(/imr,jmr,lmloc/)) ) then write (gol,'("shape of output grid not ok:")'); call goErr write (gol,'(" shape(rm_k) : ",3i4)') shape(rm_k); call goErr write (gol,'(" imr,jmr,lmloc : ",3i4)') imr,jmr,lmloc; call goErr TRACEBACK; status=1; return end if ! transported or short lived tracer ? if ( (itr >= 1) .and. (itr <= ntracet) ) then ! how distributed ? select case ( previous_par(region) ) ! rm is parallel over tracers: case ( 'tracer' ) ! local tracer index if requested tracer is on this pe, ! or just the first local index on other pe's to avoid array bound problems; ! allocate storage for 3d field on target pe: if ( tracer_active(itr) ) then itr_loc = tracer_loc(itr) allocate( x_t(imr,jmr,lmr) ) else itr_loc = 1 allocate( x_t(1,1,1) ) end if ! collect rm_k arrays available on each processor ! into x_t on processor tracer_id(itr) call Par_Gather_From_Levels( rm_k, x_t, tracer_id(itr), status ) IF_NOTOK_RETURN(status=1) ! now for processor with requested tracer only: if ( tracer_active(itr) ) then ! rm_k was in given unit; convert x_t to kg : select case ( unit ) case ( 'kg' ) ! x_t already in kg ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_t = x_t / plc_kg_from_tm(itr) case ( 'mass-mixing-ratio' ) ! mass mixing ratio ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_t = x_t / plc_kg_from_tm(itr) x_t = x_t * m_dat(region)%data(1:imr,1:jmr,1:lmr) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_t = x_t / plc_kg_from_tm(itr) x_t = x_t * m_dat(region)%data(1:imr,1:jmr,1:lmr) / fscale(itr) case default write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr TRACEBACK; status=1; return end select ! replace new mass x_t by difference with rm : x_t = x_t - mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc) ! kg ! add change in kg to tracer arrays, adjust slopes: do l = 1, lmr do j = 1, jmr do i = 1, imr ! adjust tracer arrays for local indices: call AdjustTracer( x_t(i,j,l), region, i, j, l, itr_loc, status ) IF_NOTOK_RETURN(status=1) end do end do end do end if ! tracer active ! clear deallocate( x_t ) ! rm is parallel over levels: case ( 'levels' ) ! storage for tracer change: allocate( x_k(imr,jmr,lmloc) ) ! copy input: x_k = rm_k ! rm_k was in given unit; convert x_k to kg : select case ( unit ) case ( 'kg' ) ! x_k already in kg ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_k = x_k / plc_kg_from_tm(itr) case ( 'mass-mixing-ratio' ) ! from mass mixing ratio to tracer mass: x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_k = x_k / plc_kg_from_tm(itr) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) / fscale(itr) ! apply unit conversion from (kg plc tracer) to (kg tm tracer) x_k = x_k / plc_kg_from_tm(itr) case default write (gol,'("unsupported unit for par level : ",a)') unit; call goErr TRACEBACK; status=1; return end select ! replace x_k by change between new and old tracer mass: x_k = x_k - mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) ! add change in kg to tracer arrays, adjust slopes: do l = 1, lmloc do j = 1, jmr do i = 1, imr ! adjust tracer arrays for local indices: call AdjustTracer( x_k(i,j,l), region, i, j, l, itr, status ) IF_NOTOK_RETURN(status=1) end do end do end do ! clear deallocate( x_k ) ! error ... case default write (gol,'("unsupported distribution:")'); call goPr write (gol,'(" region : ",i2)') region ;call goPr write (gol,'(" previous_par : ",a)') previous_par(region) ;call goPr TRACEBACK; status=1; return end select else if ( (itr > ntracet) .and. (itr <= ntrace) ) then ! storage for tracer change: allocate( x_k(imr,jmr,lmloc) ) ! copy input: x_k = rm_k ! rm_k was in given unit; convert x_k to kg : select case ( unit ) case ( 'kg' ) ! x_k already in kg case ( 'mass-mixing-ratio' ) ! mass mixing ratio to tracer mass: x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) / fscale(itr) case default write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr TRACEBACK; status=1; return end select ! replace short-lived tracer: chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) = x_k ! kg ! clear deallocate( x_k ) else write (gol,'("unsupported tracer index : ",i4)') itr; call goErr write (gol,'(" ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr TRACEBACK; status=1; return end if ! transported or chem-only ! ok status = 0 end subroutine Tracer_Fill_Slabs #endif ! *** ! ! Fill array on this pe with tracer (all levels) ! #ifdef with_mpi subroutine Tracer_Get_Here( region, itr, unit, rm_t, dest_id, status ) use mpi_const , only : MPI_DOUBLE_PRECISION, MPI_STATUS_SIZE, MPI_SUCCESS use dims , only : im, jm, lm use chem_param , only : ntrace, ntracet, fscale use partools , only : myid, root use partools , only : previous_par, tracer_loc, tracer_id, tracer_active use partools , only : localComm ! use partools , only : lmloc, lmar, offsetl ! use partools , only : Par_Scatter_Over_Levels use meteodata , only : m_dat ! --- in/out ------------------------------------- integer, intent(in) :: region integer, intent(in) :: itr character(len=*), intent(in) :: unit real(8), intent(out) :: rm_t(:,:,:) integer, intent(in) :: dest_id integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Tracer_Get_Here' ! --- local --------------------------------- integer :: imr, jmr, lmr integer :: itr_loc integer :: mpi_status(MPI_STATUS_SIZE) integer :: mpi_tag ! --- begin ---------------------------------- ! local grid size imr = im(region) jmr = jm(region) lmr = lm(region) ! check input: if ( any(shape(rm_t)/=(/imr,jmr,lmr/)) ) then write (gol,'("shape of output grid not ok:")'); call goErr write (gol,'(" region : ",i4)') region; call goErr write (gol,'(" imr,jmr,lmr : ",3i4)') imr,jmr,lmr; call goErr write (gol,'(" shape(rm_t) : ",3i4)') shape(rm_t); call goErr TRACEBACK; status=1; return end if ! transported or short lived tracer ? if ( (itr >= 1) .and. (itr <= ntracet) ) then ! how distributed ? select case ( previous_par(region) ) case ( 'tracer' ) ! local tracer index: itr_loc = tracer_loc(itr) ! already on destination pe ? if ( tracer_id(itr) == dest_id ) then ! is this the destination pe ? if ( myid == dest_id ) then ! copy into destination: rm_t = mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc) end if else ! transfer from source to destination pe ... ! set some tag: mpi_tag = 1234 ! is this the source or destination pe ? if ( myid == tracer_id(itr) ) then ! this is the source; send to destination: call MPI_Send( mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc), & imr*jmr*lmr, MPI_DOUBLE_PRECISION, & dest_id, mpi_tag, localComm, status ) IF_MPI_NOTOK_RETURN(status=1) else if ( myid == dest_id ) then ! this is the destination; receive from source: call MPI_Recv( rm_t, imr*jmr*lmr, MPI_DOUBLE_PRECISION, & tracer_id(itr), mpi_tag, localComm, mpi_status, status ) IF_MPI_NOTOK_RETURN(status=1) end if end if ! rm_t now filled on destination; post process if necessary: if ( myid == dest_id ) then ! change unit ? select case ( unit ) case ( 'kg' ) ! rm already in kg case ( 'mass-mixing-ratio' ) ! mass mixing ratio rm_t = rm_t / m_dat(region)%data(1:imr,1:jmr,1:lmr) case ( 'volume-mixing-ratio' ) ! volume mixing ratio = mass mixing ratio * fscale rm_t = rm_t / m_dat(region)%data(1:imr,1:jmr,1:lmr) * fscale(itr) case default write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr TRACEBACK; status=1; return end select end if ! case ( 'levels' ) ! ! ! array already distributed over levels; just copy: ! rm_k = mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) ! ! ! change unit ? ! select case ( unit ) ! case ( 'kg' ) ! ! rm already in kg ! case ( 'mass-mixing-ratio' ) ! ! mass mixing ratio ! rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) ! case ( 'volume-mixing-ratio' ) ! ! volume mixing ratio = mass mixing ratio * fscale ! rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr) ! case default ! write (gol,'("unsupported unit for par level : ",a)') unit; call goErr ! TRACEBACK; status=1; return ! end select case default write (gol,'("unsupported distribution:")'); call goPr write (gol,'(" region : ",i2)') region ;call goPr write (gol,'(" previous_par : ",a)') previous_par(region) ;call goPr TRACEBACK; status=1; return end select else if ( (itr > ntracet) .and. (itr <= ntrace) ) then write (gol,'("routine not implemented for not-transported tracers yet ...")'); call goErr TRACEBACK; status=1; return ! ! extract short-lived tracer: ! rm_k = chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) ! ! ! change unit ? ! select case ( unit ) ! case ( 'kg' ) ! ! rm already in kg ! case ( 'mass-mixing-ratio' ) ! ! mass mixing ratio ! rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) ! case ( 'volume-mixing-ratio' ) ! ! volume mixing ratio = mass mixing ratio * fscale ! rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr) ! case default ! write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr ! TRACEBACK; status=1; return ! end select else write (gol,'("unsupported tracer index : ",i4)') itr; call goErr write (gol,'(" ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr TRACEBACK; status=1; return end if ! transported or chem-only ! ok status = 0 end subroutine Tracer_Get_Here #endif end module tracer_data