! #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: TOOLBOX ! ! !DESCRIPTION: General tools for chemistry/emission treatment. !\\ !\\ ! !INTERFACE: ! MODULE TOOLBOX ! ! !USES: ! use GO, only : gol, goErr, goPr, goLabel, goBug use dims, only : okdebug IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: zfarr ! calculation of Arrhenius expression for rate constants public :: zf3bod ! calculation of the reaction rate of 3 body reaction public :: ltropo ! find tropopause level from T gradient public :: ltropo_ifs ! find tropopause level from T gradient, IFS way public :: lvlpress ! return index of highest level below a pressure value public :: dumpfield ! write 4D field (real) to file (debug tool) public :: dumpfieldi ! write 4D field (real) to file (debug tool) public :: coarsen_emission ! convert GLOBAL field to each region public :: escape_tm ! error handler public :: distribute_emis2D ! distribute vertically 2D field [work on MPI tiles] public :: distribute1x1 ! distribute vertically GLOBAL 1x1 2D field public :: distribute1x1b ! alternate distribute1x1 public :: tropospheric_columns ! integrate tropospheric ozone ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'toolbox' ! ! !INTERFACE: ! interface coarsen_emission module procedure coarsen_emission_1d module procedure coarsen_emission_2d module procedure coarsen_emission_3d module procedure coarsen_emission_d23 end interface ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - optional input to distribute_emis2D; protex ! 24 Jan 2012 - P. Le Sager - fixed coarsen_emission_2d and 3d for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------- !BOP ! ! !FUNCTION: ltropo ! ! !DESCRIPTION: determine tropopause level from temperature gradient ! reference: WMO /Bram Bregman !\\ !\\ ! !INTERFACE: ! integer function ltropo(region,tx,gphx,lmx) ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region integer,intent(in) :: lmx real,dimension(lmx),intent(in) :: tx real,dimension(lmx+1),intent(in) :: gphx ! ! !REVISION HISTORY: ! programmed by fd 01-11-1996 ! changed to function fd 12-07-2002 ! 16 Jul 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC integer :: ltmin,ltmax,l real :: dz,dt ! ltropo is highest tropospheric level ! above is defined as stratosphere ! min tropopause level 450 hPa (ca. 8 km) ltmin=lvlpress(region,45000.,98400.) ! max tropopause level 70 hPa (ca. 20 km) ltmax=lvlpress(region,7000.,98400.) ltropo=ltmin do l=ltmin,ltmax dz=(gphx(l)-gphx(l-2))/2. dt=(tx(l)-tx(l-1))/dz if ( dt < -0.002) ltropo=l !wmo 85 criterium for stratosphere ! increase upper tropospheric level end do !l end function ltropo !EOC !----------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------- !BOP ! ! !FUNCTION: ltropo_ifs ! ! !DESCRIPTION: determine tropopause level from temperature gradient ! as it is done in IFS cmip6_strato_aero_process.F90 (EC-Earth 3.2.3) !\\ !\\ ! !INTERFACE: ! integer function ltropo_ifs(region,i,j,tx,lmx) use binas, only : rgas_air, grav use dims, only : at, bt, lm use meteodata, only : phlb_dat, m_dat ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region, i, j integer,intent(in) :: lmx real,dimension(lmx),intent(in) :: tx ! ! !REVISION HISTORY: ! ! !REMARKS: ! returns the highest level *in* the troposhere, like !EOP !------------------------------------------------------------------------ !BOC integer :: ltmin, ltmax, l, ltropo_ifs_default, jlev real :: dz, dt, lapse_rate, lapse_rate_mean real, dimension(lmx) :: papf real, dimension(lmx+1) :: zpth logical :: ltest real, dimension(:,:,:), pointer :: phlb phlb => phlb_dat(region)%data ltest=.false. ! temperature at half levels do l=2,lmx zpth(l) = (tx(l) + tx(l-1))/2. enddo zpth(1)=tx(1) zpth(lmx+1)=tx(lmx) ! pressure at full-levels do l=1,lmx papf(l)=( phlb(i,j,l) + phlb(i,j,l+1) ) / 2. enddo ltropo_ifs = lmx ltropo_ifs_default = 1 DO l=1,lmx ! register highest troposheric level if (papf(l) .gt. 10000.) ltropo_ifs_default = max(ltropo_ifs_default,l) if (papf(l)>7500. .and. papf(l)<55000.) then ! the lapse rate is computed for each layer (here we use the gas ! equation and the hydrostatic approximation): lapse_rate = (zpth(l)-zpth(l+1)) / (phlb(i,j,l)-phlb(i,j,l+1)) * papf(l)/tx(l)*grav/rgas_air if (lapse_rate <= 0.002) then dz=0 jlev=l do while (dz <= 2000) ! compute the delta(p) corresponding to 2km dz=dz+(papf(jlev)-papf(jlev+1))/phlb(i,j,jlev+1)*zpth(jlev+1)*rgas_air/grav jlev = jlev+1 enddo lapse_rate_mean = (tx(l)-tx(jlev))/dz ! lapse rate over 2 km if (lapse_rate_mean <= 0.002) then ltropo_ifs =min(l,ltropo_ifs) ltest=.true. endif endif endif enddo if (.not. ltest) then ltropo_ifs = ltropo_ifs_default end if end function ltropo_ifs !--------------------------------------------------------------------------- ! TM5 ! !--------------------------------------------------------------------------- !BOP ! ! !FUNCTION: lvlpress ! ! !DESCRIPTION: lvlpress determines the index of the level with a pressure ! greater than press i.e. in height below press ! determines level independent of vertical resolution lm ! uses hybrid level coefficients to determine lvlpress !\\ !\\ ! !INTERFACE: ! integer function lvlpress(region,press,press0) ! ! !USES: ! use dims, only: at, bt, lm ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region real, intent(in) :: press real, intent(in) :: press0 ! ! !REVISION HISTORY: ! programmed by Peter van Velthoven, 6-november-2000 ! 16 Jul 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC real :: zpress0, zpress integer :: l,lvl if ( press0 == 0.0 ) then zpress0 = 98400. else zpress0 = press0 endif lvl = 1 ! l increases from bottom (l=1) to top (l=lm) do l = 1, lm(region) zpress = (at(l)+at(l+1) + (bt(l)+bt(l+1))*zpress0)*0.5 if ( press < zpress ) then lvl = l endif end do lvlpress = lvl end function lvlpress !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: zfarr ! ! !DESCRIPTION: ZFARR calculation of Arrhenius expression for rate constants !\\ !\\ ! !INTERFACE: ! real function zfarr(rx1,er,ztrec) ! ! !INPUT PARAMETERS: ! real,intent(in) :: rx1,er,ztrec ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - ! !EOP !------------------------------------------------------------------------ !BOC zfarr=rx1*exp(er*ztrec) end function zfarr !EOC real function zf3bod(rx1,rx2,rx3) !---------------------------------------------------------------------- ! ! function to calculate the reaction rate of 3 body reaction ! !------------------------------------------------------------------ ! real,intent(in) :: rx1,rx2,rx3 ! zf3bod=rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2)) ! end function zf3bod !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DumpField ! ! !DESCRIPTION: Write 4D field (type real) to HDF file !\\ !\\ ! !INTERFACE: ! subroutine dumpfield(flag,fname,im,jm,lm,nt,field,namfield) ! ! !USES: ! #ifdef with_hdf4 use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: im, jm, lm, nt integer, intent(in) :: flag real, dimension(im,jm,lm,nt), intent(in) :: field ! 4D field character(len=*), intent(in) :: fname ! file name character(len=*), intent(in) :: namfield ! field name ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - protex ! !EOP !------------------------------------------------------------------------ !BOC #ifdef with_hdf4 integer :: sfstart, io, sfend if ( flag==0 ) then io = sfstart(fname,DFACC_CREATE) else io = sfstart(fname,DFACC_WRITE) if ( io == -1 ) io = sfstart(fname,DFACC_CREATE) end if call io_write(io,im,'X',jm,'Y',lm,'Z',nt,'N',field,namfield) io = sfend(io) #endif end subroutine dumpfield !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: dumpfieldi ! ! !DESCRIPTION: Write 4D field (type integer) to HDF file !\\ !\\ ! !INTERFACE: ! subroutine dumpfieldi(flag,fname,im,jm,lm,nt,field,namfield) ! ! !USES: ! #ifdef with_hdf4 use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write #endif ! ! !INPUT PARAMETERS: ! integer,intent(in) :: im, jm, lm, nt integer,intent(in) :: flag integer,dimension(im,jm,lm,nt),intent(in) :: field ! 4D field character(len=*),intent(in) :: fname ! file name character(len=*),intent(in) :: namfield ! field name ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - protex ! !EOP !------------------------------------------------------------------------ !BOC #ifdef with_hdf4 integer :: sfstart,io,sfend if ( flag == 0 ) then io = sfstart(fname,DFACC_CREATE) else io = sfstart(fname,DFACC_WRITE) if (io==-1) io = sfstart(fname,DFACC_CREATE) end if call io_write(io,im,'X',jm,'Y',lm,'Z',nt,'N',field,namfield) io = sfend(io) #endif end subroutine dumpfieldi !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: coarsen_emission_1d ! ! !DESCRIPTION: Transform the 1D field available on e.g. 1 degree resolution ! to the desired zoom geometry !\\ !\\ ! !INTERFACE: ! subroutine coarsen_emission_1d( name_field, jm_emis, field_in, field, avg, status) ! ! name_field : name, only for printing reasons ! jm_emis : the dimension of the input field ! field_in : the input field ! field : type d2_data (defined in chem_param) ! avg : flags the mode: avg = 1 means that a surface area ! weighted area is calulated (e.g. land fraction) ! avg = 0 means that the sum of the underlying field ! is taken (e.g. emissions in kg/month). ! ! !USES: ! use dims, only: nregions, dy, yref, ybeg, jm, dxy11, idate use global_types, only: d2_data implicit none ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: name_field integer, intent(in) :: avg ! 0=no 1=yes integer, intent(in) :: jm_emis real, dimension(jm_emis), intent(in) :: field_in ! ! !INPUT/OUTPUT PARAMETERS: ! type(d2_data), dimension(nregions), target :: field ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Jun 2012 - P. Le Sager - got rid of escape_tm ! ! !REMARKS: ! (1) FIXME : should work as is for lon-lat MPI domain decomposition, but ! not tested. ! (2) This is limited to a 1 degree zonal global input [-90,+90] ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/coarsen_emission_1d' real, dimension(:), pointer :: coarse_field integer :: region real :: scale,w,wtot integer :: ystart integer :: yres integer :: j,j_index,jj integer :: jm_start ! start do region=1,nregions ! main loop over regions yres=dy/yref(region) if ( yres < 1 ) then write (gol,'("1 degree minimum resolution in emissions")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if if ( jm_emis /= 180 ) then write (gol,'("input resolution should be 1 degree")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if !WP! find index of southmost gridpoint based on 1x1 degree jm_start=ybeg(region)+91 if ( jm_start <= -1 ) then write (gol,'("requested emission fields outside valid region")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if coarse_field => field(region)%d2 do j=1,jm(region) !cycle through 1x1 grid with steps for this region j_index=jm_start+(j-1)*yres coarse_field(j) = 0.0 wtot = 0.0 do jj=0,yres-1 if(avg==1) then w = dxy11(j_index+jj) wtot = wtot + w coarse_field(j)=coarse_field(j) + field_in(j_index+jj)*w else coarse_field(j)=coarse_field(j) + field_in(j_index+jj) end if end do if ( avg == 1 ) coarse_field(j) = coarse_field(j)/wtot end do !if ( avg == 0 ) then ! write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') & ! ' coarsen_emission_1d: region ',region,name_field, & ! ' Field coarsened month :',idate(2),'total:',& ! sum(coarse_field) ; call goPr !else ! write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') & ! ' coarsen_emission_1d: region:',region,name_field, & ! ' Field averaged month :',idate(2) ; call goPr !end if nullify(coarse_field) end do ! loop over regions.... end subroutine coarsen_emission_1d !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COARSEN_EMISSION_2D ! ! !DESCRIPTION: Transform a 2D *global* emissions, at a given resolution, ! to the geometry of all model regions. !\\ !\\ ! !INTERFACE: ! SUBROUTINE COARSEN_EMISSION_2D(name_field, im_emis, jm_emis, field_in, field, avg, status) ! ! !USES: ! use GO, only : gol, goPr, goErr use Grid, only : TllGridInfo, Init, Done, FillGrid use dims, only : nregions use global_types, only : emis_data use meteodata, only : global_lli ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: name_field ! name, only for printing reasons integer, intent(in) :: im_emis, jm_emis ! grid size of global input field real, intent(in) :: field_in(im_emis,jm_emis) ! input field integer, intent(in) :: avg ! avg = 1 means that a surface area ! weighted area is calulated (e.g. land fraction) ! avg = 0 means that the sum of the underlying field ! is taken (e.g. emissions in kg/month). ! ! !INPUT/OUTPUT PARAMETERS: ! type(emis_data), target :: field(nregions) ! per region, %surf(im,jm) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/coarsen_emission_2d' integer :: region type(TllGridInfo) :: lli_in character(len=32) :: combkey ! --- begin --------------------------------------- ! define input grid call Init( lli_in, -180.0+0.5*360.0/im_emis, 360.0/im_emis, im_emis, & -90.0+0.5*180.0/jm_emis, 180.0/jm_emis, jm_emis, status ) IF_NOTOK_RETURN(status=1) ! sum or average status=0 select case ( avg ) case ( 0 ) ; combkey = 'sum' case ( 1 ) ; combkey = 'area-aver' case default write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr status=1 end select IF_NOTOK_RETURN(status=1) ! main loop over regions do region = 1, nregions ! convert grid: call FillGrid( global_lli(region), 'n', field(region)%surf, & lli_in, 'n', field_in, trim(combkey), status) IF_NOTOK_RETURN(status=1) end do ! done call Done( lli_in, status ) IF_NOTOK_RETURN(status=1) END SUBROUTINE COARSEN_EMISSION_2D !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COARSEN_EMISSION_3D ! ! !DESCRIPTION: Transform a 3D *global* emissions, at a given resolution, ! to the geometry of all model regions. !\\ !\\ ! !INTERFACE: ! SUBROUTINE COARSEN_EMISSION_3D( name_field, im_emis, jm_emis, lm_emis, & field_in, field, avg, status ) ! ! !USES: ! use GO, only : gol, goPr, goErr use Grid, only : TllGridInfo, Init, Done, FillGrid use dims, only : nregions use global_types, only : d3_data use meteodata, only : global_lli ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: name_field integer, intent(in) :: avg integer, intent(in) :: im_emis, jm_emis, lm_emis real, intent(in) :: field_in(im_emis,jm_emis,lm_emis) ! ! !INPUT/OUTPUT PARAMETERS: ! type(d3_data), target :: field(nregions) ! contains, per region, 3d field %d3(im,jm,lm) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/coarsen_emission_3d' integer :: region integer :: l type(TllGridInfo) :: lli_in character(len=32) :: combkey ! --- begin --------------------------------------- ! define input grid call Init( lli_in, -180.0+0.5*360.0/im_emis, 360.0/im_emis, im_emis, & -90.0+0.5*180.0/jm_emis, 180.0/jm_emis, jm_emis, status ) IF_NOTOK_RETURN(status=1) ! sum or average status=0 select case ( avg ) case ( 0 ) ; combkey = 'sum' case ( 1 ) ; combkey = 'area-aver' case default write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr status=1 end select IF_NOTOK_RETURN(status=1) ! main loop over regions do region = 1, nregions ! convert grid: do l = 1, lm_emis call FillGrid( global_lli(region), 'n', field(region)%d3(:,:,l), & lli_in, 'n', field_in(:,:,l), trim(combkey), status) IF_NOTOK_RETURN(status=1) end do !! info ... !select case ( combkey ) ! case ( 'aver', 'area-aver' ) ! write (gol,'("coarsen_emission : region ",i1," ",a," field averaged")') & ! region, trim(name_field); call goPr ! case ( 'sum' ) ! write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distr.; total ",es12.3)') & ! region, trim(name_field), sum(field(region)%d3); call goPr ! case default ! write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distributed")') & ! region, trim(name_field); call goPr !end select end do ! regions ! done call Done( lli_in, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE COARSEN_EMISSION_3D !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COARSEN_EMISSION_D23 ! ! !DESCRIPTION: Transform the 2D emissions available on lat-pressure resolution ! to the desired zoom geometry !\\ !\\ ! !INTERFACE: ! subroutine coarsen_emission_d23( name_field, jm_emis, lm_emis, & field_in, field, avg, status ) ! ! name_field : name, only for printing reasons ! jm_emis : the y-dimension of the input field ! lm_emis : the z-dimension of the input field ! field_in : the input field ! field : type d23_data (defined in chem_param) ! avg : flags the mode: avg = 1 means that a surface area ! weighted area is calulated (e.g. land fraction) ! avg = 0 means that the sum of the underlying field ! is taken (e.g. emissions in kg/month). ! ! !USES: ! use dims, only : nregions, dy, yref, ybeg, jm, lm, dxy11, idate use global_types, only : d23_data ! ! !INPUT PARAMETERS: ! character(len=*),intent(in) :: name_field integer,intent(in) :: avg !0=no 1=yes integer,intent(in) :: jm_emis,lm_emis real,dimension(jm_emis,lm_emis),intent(in) :: field_in ! ! !INPUT/OUTPUT PARAMETERS: ! type(d23_data),dimension(nregions),target :: field ! contains, per region, field %d23(jm,lm) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Jun 2012 - P. Le Sager - got rid of escape_tm ! ! !REMARKS: ! (1) FIXME : should work as is for lon-lat MPI domain decomposition, but ! not tested. ! (2) First dimension of input is limited to a 1 degree global [-90,+90] ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/coarsen_emission_d23' real, dimension(:,:), pointer :: coarse_field integer :: region real :: scale,w,wtot integer :: ystart integer :: yres integer :: j,j_index,jj,l integer :: jm_start ! start do region=1,nregions ! main loop over regions yres=dy/yref(region) if ( yres < 1 ) then write (gol,'("1 degree minimum resolution in emissions")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if if ( jm_emis /= 180 ) then write (gol,'("input resolution should be 1 degree")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if !WP! find index of southmost gridpoint based on 1x1 degree jm_start=ybeg(region)+91 if ( jm_start <= -1 ) then write (gol,'("requested emission fields outside valid region")'); call goErr status=1 IF_NOTOK_RETURN(status=1) end if coarse_field => field(region)%d23 do l=1,lm(region) do j=1,jm(region) !cycle through 1x1 grid with steps for this region j_index=jm_start+(j-1)*yres coarse_field(j,l) = 0.0 wtot = 0.0 do jj=0,yres-1 if(avg==1) then w = dxy11(j_index+jj) wtot = wtot + w coarse_field(j,l)=coarse_field(j,l) + & field_in(j_index+jj,l)*w else coarse_field(j,l)=coarse_field(j,l) + & field_in(j_index+jj,l) end if end do if ( avg == 1 ) coarse_field(j,l) = coarse_field(j,l)/wtot end do !j end do !l !if ( avg == 0 ) then ! ! write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') & ! ' coarsen_emission_d23: region ',region,name_field// & ! ' Field coarsened month ',idate(2),' total ',& ! sum(coarse_field); call goPr !else ! write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') & ! ' coarsen_emission_d23: region ',region,name_field// & ! ' Field averaged month ',idate(2); call goPr !end if nullify(coarse_field) end do ! loop over regions.... end subroutine coarsen_emission_d23 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: ESCAPE_TM ! ! !DESCRIPTION: abnormal termination of a model run !\\ !\\ ! !INTERFACE: ! subroutine escape_tm( msg ) ! ! !USES: ! use GO, only : GO_Print_Done use dims, only : okdebug, kmain, itau use datetime, only : wrtgol_tstamp use partools, only : Par_StopMPI #ifdef with_prism #ifdef oasis3 use mod_oasis #endif use TM5_Prism, only : comp_id #endif ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: msg ! character string to be printed on unit gol ! ! !REVISION HISTORY: ! 19 Jun 2012 - P. Le Sager - use par_stopmpi instead of stop; protex ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/escape_tm' integer :: status ! --- begin --------------------------- write (gol,'(1x,39("--"))'); call goErr write (gol,'(1x,39("--"))'); call goErr call wrtgol_tstamp(itau,'ERROR - ESCAPE_TM'); call goErr write (gol,'(1x,a)') trim(msg); call goErr write (gol,'(1x,39("--"))'); call goErr write (gol,'(1x,39("--"))'); call goErr ! try to display at least the messages in a proper way ... call GO_Print_Done( status ) if (status/=0) write (*,'("ERROR status from GO_Print_Done : ",i6)') status #ifdef with_prism #ifdef oasis3 call oasis_abort( comp_id, rname, 'called prism_abort' ) #endif #endif call Par_StopMPI end subroutine escape_tm !EOC !----------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: DISTRIBUTE_EMIS2D ! ! !DESCRIPTION: subroutine to distribute the emissions given as a TM5 2D field ! into a TM5 3D emission field. Hlow and Hhigh are the bounds of ! the 2D emission fields. They give the height RELATIVE to oro ! From that the distribution is calculated ! employing the geopotential height (relative to oro) !\\ !\\ ! !INTERFACE: ! SUBROUTINE DISTRIBUTE_EMIS2D( region, emis2D, emis3D, hlow, hhigh, status, xfrac, lat_start, lat_end) ! ! !USES: ! use tm5_distgrid, only : Get_DistGrid, dgrid use dims, only : lm, im, jm, dy, yref, ybeg use meteodata, only : gph_dat use global_types, only : d3_data, emis_data ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region type(emis_data), intent(in), target :: emis2D ! 2D emission field (kg/gridbox/month) type(d3_data), intent(inout), target :: emis3D ! 3D emission field (kg/gridbox/month) real, intent(in) :: hlow ! lowest level (m) real, intent(in) :: hhigh ! highest level (m) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status real, intent(in), optional :: xfrac ! fraction of emissions to put b/w HLOW and HHIGH (default=1) real, intent(in), optional :: lat_start ! lower limit (default=-90) of application domain real, intent(in), optional :: lat_end ! upper limit (default=+90) of application domain ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - Added Lat_start and lat_end optional inputs. ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !---------------------------------------------------------------------- !BOC character(len=*), parameter :: rname = mname//'/distribute_emis2D' real, dimension(:,:,:),pointer :: gph ! geopotential height (m) real, dimension(:,:,:),pointer :: e3d ! 3D emissions real, dimension(:,:),pointer :: e2d ! 2D emissions integer :: i,j,l,j_st,j_end real, dimension(lm(1)+1) :: height real, dimension(lm(1)) :: dz real :: dy1,dze real :: totw, f, tot2d, tot3db, tot3da, fraction integer, dimension(3) :: ubound_e3d integer :: lmmax real :: hhighb real :: lat_low,lat_high integer :: i1, i2, j1, j2 real :: l_status, g_status ! --- begin --------------------------------------------- status = 0 if (present(xfrac)) then fraction = xfrac else fraction = 1.0 endif if (fraction < spacing(fraction)) return ! degenerated cases if (present(lat_start)) then lat_low = lat_start else lat_low = -90. endif if (present(lat_end)) then lat_high = lat_end else lat_high = 90. endif CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! get indices of application-region dy1=dy/yref(region) j_st = floor((lat_low - ybeg(region))/dy1) + 1 j_end = floor((lat_high - ybeg(region))/dy1) + 1 ! is tile in appl. region? if ((j_end .lt. j1 ) .or. (j_st .gt. j2) ) then return end if ! limit range if (j_end .gt. j2 ) j_end = j2 if (j_st .lt. j1 ) j_st = j1 ! check inputs dze = hhigh - hlow if ( (hhigh <= 0.0) .or. (hlow < 0.0) .or. (dze < 0.0) ) then write (gol,'("found strange emission heights:")'); call goErr write (gol,'(" hhigh (> 0.0 ?) : ",es12.4)') hhigh; call goErr write (gol,'(" hlow (>= 0.0 ?) : ",es12.4)') hlow; call goErr write (gol,'(" hhigh-hlow (> 0.0 ?) : ",es12.4)') dze; call goErr TRACEBACK; status=1; return end if ! init nullify(gph) nullify(e2d) nullify(e3d) gph => gph_dat(region)%data e2d => emis2d%surf e3d => emis3d%d3 ! get highest possible level ubound_e3d = ubound(e3d) lmmax = ubound_e3d(3) ! total of to-be-redistributed before redistribution tot2d = sum(e2d(:,j_st:j_end)*fraction) ! total of target array before adding redistributed data tot3db = sum(e3d(:,j_st:j_end,:)) ! do work jlo: do j=j_st,j_end do i=i1, i2 ! zero based height: height(1) = 0.0 do l=1, lm(region) dz(l) = gph(i,j,l+1) - gph(i,j,l) height(l+1) = height(l) + dz(l) enddo totw = 0.0 if( hhigh > height(lmmax+1) ) then write (gol,'("try to put emissions higher than array allows:")') ; call goErr write (gol,'(" hhigh : ",es12.4)') hhigh ; call goErr write (gol,'(" height(lmmax+1) : ",es12.4)') height(lmmax+1) ; call goErr write (gol,'(" at i/j=",i3,1x,i3)') i,j ; call goErr do l = 1, size(height) write (gol,*) 'height : ', l, height(l); call goErr end do do l = 1, size(gph,3) write (gol,*) 'gph : ', l, gph(i,j,l); call goErr end do TRACEBACK; status=1 exit jlo endif zz: do l=1, lmmax if (hhigh > height(l+1)) then if ( hlow < height(l) ) then f = dz(l)/dze totw = totw + f e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j) else if( hlow < height(l+1)) then f = (height(l+1)-hlow)/dze totw = totw + f e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j) endif else if ( hlow < height(l)) then f = (hhigh - height(l))/dze totw = totw + f e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j) else totw = totw + 1.0 e3d(i,j,l) = e3d(i,j,l) + fraction*e2d(i,j) endif exit zz endif enddo zz if ( abs(totw-1.0) > 1e-14 ) then write (gol,'("sum weight /= 1 : ",es12.4)') totw-1.0; call goErr TRACEBACK; status=1 exit jlo end if enddo !i enddo jlo !j IF_NOTOK_RETURN(status=1) ! Total of target array after adding redistributed data, and check ! DO NOT FIXME ?? : the following test will fail if tot2d << tot3da (strictly: if ! tot2d is less than the last significant digit of tot3da). However this ! is good, since we actually expect that tot3da to be 0 (as used in ! chemistry): so we can catch some bad initialization. tot3da = sum(e3d(:,j_st:j_end,:)) if (abs((tot3da-tot3db)-tot2d) > 1e-8*abs(tot2d) ) then write (gol,'("emissions have not been distributed mass-conserving")'); call goErr write(gol,*)"totals to add : ", tot2d, tot3db ; call goErr write(gol,*)"total after : ", tot3da ; call goErr write(gol,*)"with bounds : ", j1, j2, j_st, j_end; call goErr write(gol,*)"shapes e3d",shape(e3d) ; call goErr write(gol,*)"lbound e3d",lbound(e3d); call goErr write(gol,*)"ubound e3d",ubound(e3d); call goErr write(gol,*)"minmax e3d",minval(e3d), maxval(e3d); call goErr write(gol,*)"shapes e2d",shape(e2d) ; call goErr write(gol,*)"lbound e2d",lbound(e2d); call goErr write(gol,*)"ubound e2d",ubound(e2d); call goErr write(gol,*)"minmax e2d",minval(e2d), maxval(e2d); call goErr TRACEBACK; status=1; return end if ! done nullify(gph) nullify(e2d) nullify(e3d) END SUBROUTINE DISTRIBUTE_EMIS2D !EOC !--------------------------------------------------------------------------- ! TM5 ! !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: distribute1x1 ! ! !DESCRIPTION: subroutine to distribute 1x1 emissions linearly between ! hlow and hhigh. The vertical level is determined by ! the orography which is read from the surface file... ! A simple scale height vertical structure is assumed. ! ! To be called by one processor, with a global 2D oro input !\\ !\\ ! !INTERFACE: ! SUBROUTINE DISTRIBUTE1X1( emi1x1, hlow, hhigh, emis3d, oro, status, xfrac) ! ! !USES: ! use Binas, only : grav use dims, only : nlon360, nlat180, lm, itau, staggered, at, bt ! ! !INPUT PARAMETERS: ! real, dimension(nlon360,nlat180), intent(in) :: emi1x1 ! (kg/1x1 gridbox) 2D field of emissions real, dimension(nlon360,nlat180), intent(in) :: hlow ! (m) lower bound of emission real, dimension(nlon360,nlat180), intent(in) :: hhigh ! (m) higher bound of emission real, dimension(nlon360,nlat180), intent(in) :: oro ! (m m/s2) real, intent(in), optional :: xfrac ! fraction of emissions to put ! ! !INPUT/OUTPUT PARAMETERS: ! real, dimension(nlon360,nlat180,lm(1)), intent(inout) :: emis3d ! (kg/box) distributed in height integer, intent(out) :: status ! ! !REVISION HISTORY: ! 8 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/distribute1x1' real, parameter :: scalh = 8000.0 real :: p0, pt real, dimension(lm(1)) :: height, dz integer :: i,j,l real :: hh,hl,dze,totw,f,e3da,e3db, fract ! --- begin ----------------------------------- if (present(xfrac)) then fract = xfrac else fract = 1.0 endif e3db = sum(emis3d) do j=1,nlat180 do i=1,nlon360 if (fract*emi1x1(i,j) > 1e-14) then height(1) = oro(i,j)/grav p0 = 1e5*exp(-height(1)/scalh) do l=1,lm(1)-1 ! bug reported by FD: alog(0) crashes! pt = at(l+1) + bt(l+1)*p0 height(l+1) = height(1)-scalh*alog(pt/p0) dz(l) = height(l+1)-height(l) enddo hl = max(height(1),hlow(i,j)) hh = max(hhigh(i,j),height(1)) dze = hh-hl if(dze < 0.0) then status=1 IF_NOTOK_RETURN(status=1) else if ( dze == 0.0) then ! this somehow happens! hh = height(1)+1.0 hl = height(1) endif totw = 0.0 zz: do l=1, lm(1) if (hh > height(l+1)) then if ( hl < height(l) ) then f = dz(l)/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) else if( hl < height(l+1)) then f = (height(l+1)-hl)/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) endif else if ( hl < height(l)) then f = (hh - height(l))/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) else totw = totw + 1.0 emis3d(i,j,l) = emis3d(i,j,l) + fract*emi1x1(i,j) endif exit zz endif enddo zz if ( abs(totw-1.0) > 1e-14 ) then status=1 IF_NOTOK_RETURN(status=1) end if endif enddo enddo e3da = sum(emis3d) if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) then status=1 IF_NOTOK_RETURN(status=1) end if status=0 END SUBROUTINE DISTRIBUTE1X1 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DISTRIBUTE1X1B ! ! !DESCRIPTION: same as DISTRIBUTE1X1, but with scalar for HLOW and HHIGH ! ! Is it still used ???? ! ! subroutine to distribute 1x1 emissions linearly between ! hlow and hhigh. The vertical level is determined by ! the orography which is read from the surface file... ! A simple scale height vertical structure is assumed. ! same as distribute1x1 but hlow, hhigh now scalar ! ALSO: the height is now defined relative to the orography! !\\ !\\ ! !INTERFACE: ! SUBROUTINE DISTRIBUTE1X1B( emi1x1, hlow, hhigh, emis3d, oro, status, xfrac) ! ! !USES: ! use Binas, only : grav use dims, only : nlon360, nlat180, lm, itau, staggered, at, bt ! ! !INPUT PARAMETERS: ! real, dimension(nlon360,nlat180), intent(in) :: emi1x1 ! (kg/1x1 gridbox) 2D field of emissions real, intent(in) :: hlow ! (m) lower bound of emission real, intent(in) :: hhigh ! (m) higher bound of emission real, dimension(nlon360,nlat180), intent(in) :: oro ! (m m/s2) real, intent(in), optional :: xfrac ! fraction of emissions to put ! ! !INPUT/OUTPUT PARAMETERS: ! real, dimension(nlon360,nlat180,lm(1)), intent(inout) :: emis3d ! (kg/box) distributed in height integer, intent(out) :: status ! ! !REVISION HISTORY: ! 8 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/distribute1x1b' real, parameter :: scalh = 8000.0 real :: p0, pt real, dimension(lm(1)) :: height, dz integer :: i,j,l real :: hh,hl,dze,totw,f,e3da,e3db, fract, hlow_oro, hhigh_oro ! --- begin ----------------------------------- if (present(xfrac)) then fract = xfrac else fract = 1.0 endif e3db = sum(emis3d) do j=1,nlat180 do i=1,nlon360 if (fract*emi1x1(i,j) > 1e-14) then height(1) = oro(i,j)/grav hlow_oro = hlow + height(1) hhigh_oro = hhigh + height(1) p0 = 1e5*exp(-height(1)/scalh) do l=1,lm(1)-1 ! bug reported by FD: alog(0) crashes! pt = at(l+1) + bt(l+1)*p0 height(l+1) = height(1)-scalh*alog(pt/p0) dz(l) = height(l+1)-height(l) enddo hl = max(height(1),hlow_oro) hh = max(hhigh_oro,height(1)) dze = hh-hl if(dze < 0.0) then status=1 IF_NOTOK_RETURN(status=1) else if ( dze == 0.0) then ! this somehow happens! hh = height(1)+1.0 hl = height(1) endif totw = 0.0 zz: do l=1, lm(1) if (hh > height(l+1)) then if ( hl < height(l) ) then f = dz(l)/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) else if( hl < height(l+1)) then f = (height(l+1)-hl)/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) endif else if ( hl < height(l)) then f = (hh - height(l))/dze totw = totw + f emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j) else totw = totw + 1.0 emis3d(i,j,l) = emis3d(i,j,l) + fract*emi1x1(i,j) endif exit zz endif enddo zz if ( abs(totw-1.0) > 1e-14 ) then status=1 IF_NOTOK_RETURN(status=1) end if endif enddo enddo e3da = sum(emis3d) if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) then status=1 IF_NOTOK_RETURN(status=1) end if status=0 END SUBROUTINE DISTRIBUTE1X1B !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: TROPOSPHERIC_COLUMNS ! ! !DESCRIPTION: Routine to integrate tropospheric (ozone) ! note: routine is now written for Ozone, but may be changed ! to be more general. The definition of tropopause is critical for ozone. ! here, the slope is used in the interpolation. ! multiple tropopause values are ignored, but are known to ! occur. The lowest crossing of thres is used. ! In the lower atmosphere > 600 hPa, values > thres are ignored. !\\ !\\ ! !INTERFACE: ! subroutine tropospheric_columns(region, field, slope, column, thres, xmo3, status) ! ! !USES: ! use binas, only : Dobs, xmair, Avog use global_data, only : region_dat use meteodata, only : phlb_dat, m_dat use Dims, only : lm, isr, ier, jsr, jer, CheckShape, im, jm ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region in the TM5 zoom definition real, dimension(:,:,:), intent(in) :: field ! 3D field of tracer (O3) real, dimension(:,:,:), intent(in) :: slope ! 3D field of tracer z -slope (O3) ! ! !OUTPUT PARAMETERS: ! real, dimension(:,:),intent(out) :: column ! output: tropospheric column in DU real, intent(in) :: thres ! ppb threshold for tropospheric column real, intent(in) :: xmo3 ! mol mass tracer integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Jun 2012 - P. Le Sager - fix calls to CheckShape ! ! !REMARKS: ! (1) FIXME : must be adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/tropospheric_columns' #ifdef with_zoom integer,dimension(:,:),pointer :: zoomed #endif real,dimension(:,:,:),pointer :: m real,dimension(:,:,:),pointer :: phlb real,dimension(:),pointer :: dxyp real, dimension(2*lm(1)+1) :: o3a real :: o3mm, o3mmu, o3mml, o3trop, o3ml, o3mu, frac integer :: i,j,l,ip ! FIXME write (gol,'("ERROR - routine not converted to TM5-MP")') ; call goErr status=1 IF_NOTOK_RETURN(status=1) ! start call CheckShape( (/im(region),jm(region)/), shape(column), status ) call CheckShape( (/im(region),jm(region), lm(region)/), shape(field), status ) call CheckShape( (/im(region),jm(region), lm(region)/), shape(slope), status ) dxyp=> region_dat(region)%dxyp #ifdef with_zoom zoomed => region_dat(region)%zoomed #endif phlb => phlb_dat(region)%data m => m_dat(region)%data ! collect ozone on mid levels and at boundaries ! average the estimates from upper/lower gridboxes ! except for bottom and top. ! column(:,:) = 0.0 do j=jsr(region), jer(region) do i=isr(region), ier(region) #ifdef with_zoom if(zoomed(i,j)/=region) cycle #endif do l=1,lm (region) ip = 2*l-1 o3mm = xmair/xmo3*1e9*field(i,j,l)/m(i,j,l) ! level o3mmu = xmair/xmo3*1e9*(field(i,j,l)+slope(i,j,l))/m(i,j,l) ! top o3mml = xmair/xmo3*1e9*(field(i,j,l)-slope(i,j,l))/m(i,j,l) ! bottom if(l == 1) then o3a(ip) = o3mml else o3a(ip) = o3a(ip) + 0.5*o3mml endif o3a(ip+1) = o3mm if(l == lm(region)) then o3a(ip+2) = o3mmu else o3a(ip+2) = 0.5*o3mmu endif enddo o3trop = 0.0 height: do l=1,lm (region) ip = 2*l-1 ! split gridbox in upper and lower part ! for more accurate interpolation o3ml = 0.5*field(i,j,l) - 0.25*slope(i,j,l) o3mu = 0.5*field(i,j,l) + 0.25*slope(i,j,l) if (phlb(i,j,l) > 60000.0) then ! avoid surface o3>150ppb o3trop = o3trop + field(i,j,l) cycle height endif if (phlb(i,j,l+1) < 7000.0) exit height ! now about time 70 hPa at top if(o3a(ip) < thres) then ! bottom value less than thres if(o3a(ip+1) >= thres) then ! but central value is not frac = (thres - o3a(ip))/(o3a(ip+1)-o3a(ip)) o3trop = o3trop + frac*o3ml exit height else if (o3a(ip+2) >= thres) then ! but upper is not frac = (thres - o3a(ip+1))/(o3a(ip+2)-o3a(ip+1)) o3trop = o3trop + o3ml + frac*o3mu exit height else ! entire layer is not o3trop = o3trop + field(i,j,l) endif else exit height endif enddo height column(i,j) = o3trop*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs ! kg ->dobs enddo ! i enddo !j call dumpfield(0,'dump.hdf', im(region), jm(region), lm(region), 1, field, 'o3') call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region), 1, slope, 'o3slope') call dumpfield(1,'dump.hdf', im(region), jm(region), 1, 1, column, 'o3column') call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region)+1, 1, phlb, 'phlb') call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region), 1, m, 'm') call dumpfield(1,'dump.hdf', jm(region),1, 1, 1, dxyp, 'dxyp') nullify(dxyp) #ifdef with_zoom nullify(zoomed) #endif nullify(phlb) nullify(m) end subroutine tropospheric_columns !EOC END MODULE TOOLBOX