!### 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" !------------------------------------------------------------------------- ! TM5 ! !------------------------------------------------------------------------- !BOP ! ! !MODULE: toolbox ! ! !DESCRIPTION: General tools for chemistry/emission treatment. !\\ !\\ ! !INTERFACE: ! module toolbox ! ! !USES: ! use GO, only : gol, goErr, goPr, goLabel use dims, only : okdebug implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: zfarr public :: ltropo, lvlpress, print_totalmass public :: dumpfield, dumpfieldi public :: coarsen_emission, coarsen_emission_2d public :: escape_tm public :: distribute_emis2D public :: distribute1x1 public :: distribute1x1b public :: tropospheric_columns ! ! !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 ! ! !REMARKS: ! !EOP !----------------------------------------------------------------------- contains ! ============================================================================ !----------------------------------------------------------------------- ! !**** tropo determine tropopause level ! ! programmed by fd 01-11-1996 ! changed to function fd 12-07-2002 ! ! purpose ! ------- ! ! interface ! --------- ! i=ltropo(region,tx,gphx,lmx) ! ! method ! ------ ! determine tropopause level ! from temperature gradient ! ! externals ! --------- ! function lvlpress ! ! reference ! --------- ! WMO /Bram Bregman !---------------------------------------------------------------------- integer function ltropo(region,tx,gphx,lmx) ! in/out integer,intent(in) :: region integer,intent(in) :: lmx real,dimension(lmx),intent(in) :: tx real,dimension(lmx+1),intent(in) :: gphx ! local 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 !----------------------------------------------------------------------- ! !**** lvlpress determines the index of the level with a pressure ! greater than press i.e. in height below press ! ! programmed by Peter van Velthoven, 6-november-2000 ! ! purpose ! ------- ! determines level independent of vertical resolution lm ! ! interface ! --------- ! i = lvlpress(press,press0) ! ! method ! ------ ! use hybrid level coefficients to determine lvlpress ! ! input ! --------- ! press : pressure in Pa ! press0 : surface pressure in Pa ! !---------------------------------------------------------------------- integer function lvlpress(region,press,press0) use dims, only: at, bt, lm ! in/out integer, intent(in) :: region real, intent(in) :: press real, intent(in) :: press0 ! local 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 ! subroutine calc_region_coordinates(region,xlat,xlon) ! ! ! ! ! ! ! use dims, only: im, jm, dx, dy, xref, yref, xbeg, ybeg ! implicit none ! ! ! in/out ! integer, intent(in) :: region ! real, dimension(im(region)+1),intent(out) :: xlon ! real, dimension(jm(region)+1),intent(out) :: xlat ! ! ! local ! integer :: i,j,l,n ! real :: dx1,dy1 ! gridsize in degrees longitude and latitude ! ! dx1=dx/xref(region) !longitude ! dy1=dy/yref(region) !latitude ! ! ! do j=1,jm(region)+1 ! ! southern border,last value is North Border ! xlat(j)=ybeg(region)+dy1*(j-1) ! enddo ! ! do i=1,im(region)+1 ! ! western border, last value is East Border ! xlon(i)=xbeg(region)+dx1*(i-1) ! enddo ! ! end subroutine calc_region_coordinates real function zfarr(rx1,er,ztrec) !------------------------------------------------------------------ ! !**** ZFARR calculation of Arrhenius expression for rate constants ! !------------------------------------------------------------------ ! implicit none real,intent(in) :: rx1,er,ztrec ! zfarr=rx1*exp(er*ztrec) ! end function zfarr subroutine print_totalmass(region) use dims, only : im,jm,lm use global_data, only : mass_dat implicit none ! in/out integer,intent(in) :: region ! local integer :: i real :: value logical,save :: first=.true. value = & sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),1)) if ( first ) then first = .false. write (gol,*) 'print_totalmass: Initial Total mass tracer 1 ',value; call goPr else write (gol,*) 'print_totalmass: Total mass tracer 1 ',value; call goPr end if end subroutine print_totalmass subroutine dumpfield(flag,fname,im,jm,lm,nt,field,namfield) ! ! Write 4D field (type real) to HDF file ! use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write implicit none ! in/out 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 ! local 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) end subroutine dumpfield subroutine dumpfieldi(flag,fname,im,jm,lm,nt,field,namfield) ! ! Write 4D field (type integer) to HDF file ! use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write implicit none ! in/out 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 ! local 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) end subroutine dumpfieldi subroutine coarsen_emission_1d(name_field,jm_emis,field_in,field,avg) ! ! Purpose: Transform the 1D field available on e.g. 1 degree resolution ! to the desired zoom geometry ! 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). ! use dims, only: nregions, dy, yref, ybeg, jm, dxy11, idate use global_types, only: d2_data implicit none ! in/out 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 type(d2_data),dimension(nregions),target :: field ! contains per region field d2(jm) ! local 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 call escape_tm('coarsen_emission_1d:'//& ' 1 degree minimum resolution in emissions') end if if ( jm_emis /= 180 ) then call escape_tm( & 'coarsen_emission_1d: input resolution should be 1 degree') end if !WP! find index of southmost gridpoint based on 1x1 degree jm_start=ybeg(region)+91 if ( jm_start <= -1 ) then call escape_tm('coarsen_emission_1d:'//& ' requested emission fields outside valid region') 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 ! ! Purpose: Transform the 2D emissions available on ! e.g. 1x1 degree resolution ! to the desired zoom geometry ! name_field : name, only for printing reasons ! im_emis : the x-dimension of the input field ! jm_emis : the y-dimension of the input field ! field_in : the input field ! field : type emis_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). ! subroutine coarsen_emission_2d(name_field,im_emis,jm_emis,field_in,field,avg) 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 : lli ! --- in/out ------------------------------------ character(len=*),intent(in) :: name_field integer,intent(in) :: avg ! 0=no 1=yes integer,intent(in) :: im_emis,jm_emis real, intent(in) :: field_in(im_emis,jm_emis) type(emis_data), target :: field(nregions) ! contains per region field surf(im,jm) ! --- local ------------------------------------- integer :: region type(TllGridInfo) :: lli_in character(len=32) :: combkey integer :: status ! --- 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 (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from init lli') ! sum or average ? select case ( avg ) case ( 0 ) ; combkey = 'sum' case ( 1 ) ; combkey = 'area-aver' case default write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr call escape_tm('ERROR - in toolbox/coarsen_emission_2d') end select ! main loop over regions do region = 1, nregions ! convert grid: call FillGrid( lli(region), 'n', field(region)%surf, & lli_in , 'n', field_in , trim(combkey), status ) if (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from fillgrid') !! 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)%surf); 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 (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from done lli') ! ok status = 0 end subroutine coarsen_emission_2d ! ! Purpose: Transform the 3D emissions available on ! e.g. 1x1 degree resolution ! to the desired zoom geometry ! name_field : name, only for printing reasons ! im_emis : the x-dimension of the input field ! 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 d3_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). ! subroutine coarsen_emission_3d( name_field, im_emis, jm_emis, lm_emis, & field_in, field, avg ) 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 : lli ! --- in/out ------------------------------------ character(len=*), intent(in) :: name_field integer, intent(in) :: avg !0=no 1=yes integer, intent(in) :: im_emis, jm_emis, lm_emis real, intent(in) :: field_in(im_emis,jm_emis,lm_emis) type(d3_data), target :: field(nregions) ! contains per region 3d field d3(im,jm,lm) ! --- local ----------------------------- integer :: region integer :: l type(TllGridInfo) :: lli_in character(len=32) :: combkey integer :: status ! --- 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 (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from init lli') ! sum or average ? select case ( avg ) case ( 0 ) ; combkey = 'sum' case ( 1 ) ; combkey = 'area-aver' case default write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr call escape_tm('ERROR - in toolbox/coarsen_emission_2d') end select ! main loop over regions do region = 1, nregions ! convert grid: do l = 1, lm_emis call FillGrid( lli(region), 'n', field(region)%d3(:,:,l), & lli_in , 'n', field_in(:,:,l) , trim(combkey), status ) if (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from fillgrid') 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 (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from done lli') ! ok status = 0 end subroutine coarsen_emission_3d subroutine coarsen_emission_d23( name_field, jm_emis, lm_emis, & field_in, field, avg ) ! ! Purpose: Transform the 2D emissions available on lat-pressure resolution ! to the desired zoom geometry ! 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). ! use dims, only: nregions, dy, yref, ybeg, jm, lm, dxy11, idate use global_types, only: d23_data implicit none ! in/out 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 type(d23_data),dimension(nregions),target :: field ! contains per region field d23(jm,lm) ! local 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 call escape_tm('coarsen_emission_d23:'//& ' 1 degree minimum resolution in emissions') end if if ( jm_emis /= 180 ) then call escape_tm( 'coarsen_emission_d23:'//& 'input resolution should be 1 degree') end if !WP! find index of southmost gridpoint based on 1x1 degree jm_start=ybeg(region)+91 if(jm_start<=(-1)) then call escape_tm( 'coarsen_emission_d23:'//& 'requested emission fields outside valid region') 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 !-------------------------------------------------------------- ! ! abnormal termination of a model run ! ! msg character string to be printed on unit kmain ! !------------------------------------------------------------- subroutine escape_tm( msg ) use GO , only : GO_Print_Done use dims , only : okdebug, kmain, itau use datetime , only : wrtgol_tstamp #ifdef MPI use mpi_const, only : myid, root use mpi_comm , only : abortmpi #endif #ifdef with_prism #ifdef oasis3 use mod_prism_proto #endif #ifdef oasis4 use PRISM , only : prism_abort #endif use TM5_Prism, only : comp_id #endif ! --- in/out ----------------------------- character(len=*), intent(in) :: msg ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/escape_tm' ! --- local ------------------------------ 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 PRISM_Abort_Proto( comp_id, rname, 'called prism_abort' ) #endif #ifdef oasis4 call prism_abort( comp_id, rname, 'called prism_abort' ) #endif #endif #ifdef MPI call abortmpi #endif stop 'Fortran stop in escape_tm' end subroutine escape_tm ! !----------------------------------------------------------------------- ! 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 dims, only : lm, im, jm, dy, yref, ybeg use meteodata, only : gph_dat use global_types, only : d3_data, emis_data use OMP_partools, only : my_omp_domain, my_omp_max ! ! !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 real,intent(in), optional :: lat_start real,intent(in), optional :: lat_end ! latitude range for which the evaluation should be done. ! ! !REVISION HISTORY: ! 16 Jul 2010 - Achim Strunk - Added Lat_start and lat_end optional inputs. ! ! !REMARKS: ! !EOP !---------------------------------------------------------------------- !BOC ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/distribute_emis2D' ! --- local --------------------------------------------- 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 ! Openmp paramters integer :: js, je real :: l_status, g_status ! --- begin --------------------------------------------- if (present(xfrac)) then fraction = xfrac else fraction = 1.0 endif 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 dy1=dy/yref(region) j_st = floor ((lat_low - ybeg(region))/dy1) + 1 j_end = floor((lat_high - ybeg(region))/dy1) + 1 if ((j_end .lt. 1 ) .or. (j_st .gt. jm(region)) )return ! outside requested range if ((j_end .gt. jm(region)) .and. (j_st .lt. jm(region))) j_end = jm(region) ! limit region if ((j_st .lt. 1 ) .and. (j_end .gt. 1 )) j_st = 1 ! limit region 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 write (gol,'("in ",a)') rname; call goErr; status=1; return end if nullify(gph) nullify(e2d) nullify(e3d) gph => gph_dat(region)%data e2d => emis2d%surf e3d => emis3d%d3 ubound_e3d = ubound(e3d) lmmax = ubound_e3d(3) tot2d = sum(e2d(:,j_st:j_end)*fraction) tot3db = sum(e3d(:,j_st:j_end,:)) !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (region, g_status, gph, hlow, hhigh, lmmax, & !$OMP dze, e2d, e3d, fraction, j_st, j_end) & !$OMP private (i, j, js, je, l_status, height, dz, totw, gol, f) call my_omp_domain (j_st,j_end,js,je) l_status = 0.0 do j=js,je do i=1,im(region) ! 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 do l = 1, size(height) write (gol,*) 'height : ', l, height(l); call goErr end do write (gol,'("in ",a)') rname; call goErr; l_status=1.0 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 write (gol,'("in ",a)') rname; call goErr; l_status=1.0 end if enddo !i enddo !j call my_omp_max(l_status,g_status) !$OMP END PARALLEL if (g_status > 0.5) then status = 1 return end if tot3da = sum(e3d(:,j_st:j_end,:)) if (abs((tot3da-tot3db)-tot2d) > 1e-8*tot2d ) then write (gol,'("emissions have not been distributed mass-conserving")'); call goErr write (gol,'("in ",a)') rname; call goErr; status=1; return end if nullify(gph) nullify(e2d) nullify(e3d) ! ok status = 0 end subroutine distribute_emis2D !EOC ! ! !--------------------------------------------------------------------------- ! TM5 ! !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: distribute1x1 ! ! !DESCRIPTION: subroutine to distribute 1x1 emissions liniairly 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. !\\ !\\ ! !INTERFACE: ! subroutine distribute1x1(emi1x1, hlow, hhigh, emis3d, xfrac) use Binas , only : grav use dims , only : nlon360, nlat180, lm, itau, staggered, at, bt use dims , only : iglbsfc use MeteoData, only : oro_dat ! --- in/out ---------------------------------- 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,lm(1)), intent(inout) :: emis3d ! (kg/box) distributed in height real,intent(in), optional :: xfrac ! fraction of emissions to put ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/distribute1x1' !--- local ------------------------------------ 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 ----------------------------------- ! check ... if ( .not. oro_dat(iglbsfc)%used ) then write (gol,'("oro_dat(iglbsfc) not in use ...")'); call goErr TRACEBACK; call escape_tm('escape') end if 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_dat(iglbsfc)%data(i,j,1)/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 call escape_tm('In distribute1x1: dze <= 0') 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 ) & call escape_tm(' Routine distribute1x1 : sum weight /= 1') endif enddo enddo e3da = sum(emis3d) if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) & call escape_tm(' Routine distribute1x1 : distributed amount differs!') end subroutine distribute1x1 ! ! subroutine to distribute 1x1 emissions liniairly 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!!! ! subroutine distribute1x1b(emi1x1, hlow, hhigh, emis3d, xfrac) use Binas , only : grav use dims , only : nlon360, nlat180, lm, itau, staggered, at, bt use dims , only : iglbsfc use MeteoData, only : oro_dat ! --- in/out ---------------------------------- 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,lm(1)), intent(inout) :: emis3d ! (kg/box) distributed in height real,intent(in), optional :: xfrac ! fraction of emissions to put ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/distribute1x1b' !--- local ------------------------------------ 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 ----------------------------------- ! check ... if ( .not. oro_dat(iglbsfc)%used ) then write (gol,'("oro_dat(iglbsfc) not in use ...")'); call goErr TRACEBACK; call escape_tm('escape') end if 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_dat(iglbsfc)%data(i,j,1)/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 call escape_tm('In distribute1x1b: dze <= 0') 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 ) & call escape_tm(' Routine distribute1x1 : sum weight /= 1') endif enddo enddo e3da = sum(emis3d) if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) & call escape_tm(' Routine distribute1x1 : distributed amount differs!') end subroutine distribute1x1b ! 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. subroutine tropospheric_columns(region,field,slope, column,thres, xmo3) 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 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) 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 ! locals and pointers #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 ! start call CheckShape( (/im(region),jm(region)/), shape(column) ) call CheckShape( (/im(region),jm(region), lm(region)/), shape(field) ) call CheckShape( (/im(region),jm(region), lm(region)/), shape(slope) ) 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 end module toolbox