#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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" ! !################################################################# ! ! Parameterized ozone chemistry: ! o Cariolle, Cold Tracer ! o relaxation to climatology ! ! Relaxation time scale is altitude depended: ! ! | ------ 0 ! | ------ ^ ! | ^ | ! | | ColdTracer | ! | | |Cariolle ! | v | ! | ------ | ! | v ! | * ------ cariolle_plev (230 hPa) ! | o ! | * ------ tropo_plev (500 hPa) ! | o ! | o ! ----------+--------------+------> tau ! tropo_tau_lower tropo_tau_upper ! (weeks) (months) ! ! ! 22 Jun 2011 - P. Le Sager - Bug fix in O3 overhead. Plus added two ! checks on negative mixing ratio from Vincent. ! !### macro's ##################################################### ! module Chemistry_Cariolle use GO , only : gol, goPr, goErr, goLabel use dims , only : nregions use global_types, only : chem_data use Climat , only : TClimat #ifdef with_budgets use budget_global,only : nbudg, nbud_vg use chem_param, only : ncar #endif implicit none ! --- in/out ---------------------------- private public :: Cariolle_Init, Cariolle_Done, Cariolle_Chemie public :: o3clim_vmr, l_cariolle, with_trop_force ! --- const ----------------------------- character(len=*), parameter :: mname = 'cariolle' ! ! overhead fields ! ! 3d fields with total column ABOVE each level ! ! ! rm_ovh(:,:,lm ,iovh(k)) = 0.0 ! rm_ovh(:,:,lm-1,iovh(k)) = rm(:,:,lm,k) ! : ! rm_ovh(:,:,l ,iovh(k)) = sum( rm(:,:,l+1:lm,k) , 3 ) ! ! Total column is thus : ! ! rm(:,:,1,k) + rm_ovh(:,:,1,iovh(k)) ! ! number of overhead fields: integer, parameter :: novh = 1 ! indices: integer, parameter :: iovh_o3 = 1 ! --- var ------------------------------ ! overhead fields (4d arrays, parallel over layers) type(chem_data) :: ovh_dat(nregions) ! ozone climatology character(len=256) :: climat_fname type(TClimat),save :: o3clim_vmr(nregions) ! cariolle or linoz coefficients character(len=256) :: cariolle_fname type(TClimat),save :: o3coef(nregions) type(TClimat),save :: o3coefp(nregions) integer, parameter :: nc = 7 ! with cold tracer logical :: with_cold_tracer ! with forcing troposphere to climatology logical :: with_trop_force ! compute temper climat from model temperatures ? logical :: apply_taver ! minimum time scale for slower chemistry: real :: cariolle_tau_min ! pressure range over which Cariolle is applied: real :: cariolle_plev ! lowest model level for Cariolle chemistry: integer :: l_cariolle(nregions) ! Pressure range over which cold tracer is active: ! o production of cold tracer ! o hetrochem ! Loss of cold tracer is applied everywhere ... real :: coldtracer_plevs(2) ! 35.0e2,228.0e2 Pa ! tropospheric forcing to climatology: real :: tropo_plev ! Pa real :: tropo_tau ! sec #ifdef with_budgets real,dimension(nbudg, nbud_vg, ncar ),public :: budcarg #endif contains subroutine Cariolle_Init( status ) use GO , only : TrcFile, Init, Done, ReadRc use dims , only : im, jm, lm use global_data , only : rcfile use Climat , only : Init use file_cariolle, only : Cariolle_Info use toolbox , only : lvlpress ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Cariolle_Init' ! --- local ------------------------------ type(TrcFile) :: rcF integer :: region integer :: imr, jmr, lmr real :: tau_day, plev_hPa integer :: nlevs ! --- begin -------------------------------- ! ! read settings ! call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'chem.o3tracer.climat', climat_fname, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'chem.o3tracer.cariolle.coeff', cariolle_fname, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'chem.o3tracer.cariolle.taver', apply_taver, status ) IF_NOTOK_RETURN(status=1) ! pressure level below which Cariolle is applied: call ReadRc( rcF, 'chem.o3tracer.cariolle.strato.plev', plev_hPa, status ) IF_NOTOK_RETURN(status=1) cariolle_plev = plev_hPa * 1e2 ! Pa call ReadRc( rcF, 'chem.o3tracer.cariolle.strato.tau_min', tau_day, status ) IF_NOTOK_RETURN(status=1) cariolle_tau_min = tau_day * 24.0 * 3600.0 ! sec ! with cold tracer? call ReadRc( rcF, 'chem.o3tracer.coldtracer', with_cold_tracer, status ) IF_NOTOK_RETURN(status=1) ! pressure levels between which cold tracer scheme is applied: call ReadRc( rcF, 'chem.o3tracer.coldtracer.plevtop', coldtracer_plevs(1), status ) ! hPa IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'chem.o3tracer.coldtracer.plevbottom', coldtracer_plevs(2), status ) ! hPa IF_NOTOK_RETURN(status=1) coldtracer_plevs = coldtracer_plevs * 1.0e2 ! Pa call ReadRc( rcF, 'chem.o3tracer.tropforce', with_trop_force, status ) IF_NOTOK_RETURN(status=1) ! tropospheric forcing to climatology: ! on and below plev0 (hPa), use relax timescale tau0 (days) ! between plev0 and plev1 (hPa), linear interpol between tau0 and tau1 ! call ReadRc( rcF, 'chem.o3tracer.cariolle.tropo.plev', plev_hPa, status ) IF_NOTOK_RETURN(status=1) tropo_plev = plev_hPa * 1e2 ! Pa call ReadRc( rcF, 'chem.o3tracer.cariolle.tropo.tau', tau_day, status ) IF_NOTOK_RETURN(status=1) tropo_tau = tau_day * 24.0 * 3600.0 ! sec call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! ! other ! ! loop over regions: do region = 1, nregions ! regional grid sizes: imr = im(region) jmr = jm(region) lmr = lm(region) ! setup climatology; linear time interpolation call Init( o3clim_vmr(region), 'O3', 'ppv', 'linear', 1, jmr, lmr, status ) IF_NOTOK_RETURN(status=1) ! setup cariole coeff and pressure; constant within month ! grid sizes: call Cariolle_Info( status, nlevs=nlevs ) IF_NOTOK_RETURN(status=1) call Init( o3coef(region), 'coeff', 'unit', 'constant', nc, jmr, nlevs, status ) IF_NOTOK_RETURN(status=1) call Init( o3coefp(region), 'pressure', 'Pa', 'constant', 1, jmr, nlevs, status ) IF_NOTOK_RETURN(status=1) ! overhead fields (concentration arrays, thus 2 halo cells): allocate( ovh_dat(region)%rm_k(-1:imr+2,-1:jmr+2,lmr,novh) ) ! lowest model level for Cariolle chemistry ! the use of 1013.25 hPa allows reading full level pressures from ECMWF tables l_cariolle=lvlpress(region,cariolle_plev, 101325.) end do ! regions ! ok status = 0 end subroutine Cariolle_Init ! *** subroutine Cariolle_Done( status ) use Climat, only : Done ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Cariolle_Done' ! --- var ---------------------------------- integer :: region ! --- begin -------------------------------- ! loop over regions: do region = 1, nregions ! clear climatology: call Done( o3clim_vmr(region), status ) IF_NOTOK_RETURN(status=1) ! clear cariolle coeff: call Done( o3coef(region), status ) IF_NOTOK_RETURN(status=1) call Done( o3coefp(region), status ) IF_NOTOK_RETURN(status=1) ! overhead fields: deallocate( ovh_dat(region)%rm_k ) end do ! regions ! ok status = 0 end subroutine Cariolle_Done ! ================================================================ subroutine Cariolle_Chemie( region, tr, status ) use binas , only : Avog, xm_o3,xmair use GO , only : TDate, Get, DayNumber, rTotal, wrtgol use GO , only : operator(+), operator(-), operator(/) use toolbox , only : dumpfield,escape_tm use Num , only : Interp_Lin use Grid , only : AreaOper use Phys , only : cos_sza use dims , only : im, jm, lm use dims , only : isr, ier, jsr, jer use chem_param , only : fscale, io3, ipsc,xmo3 use global_data , only : mass_dat, region_dat use tracer_data , only : AdjustTracer use meteodata , only : phlb_dat, m_dat use meteo , only : lli, levi use meteo , only : temper_dat use chemistry_0d, only : Cariolle_0d, ColdTracer_0d use ParTools , only : tracer_active, lmloc, offsetl use ParTools , only : which_par, previous_par, myid use photolysis_data, only : phot_dat use mpi_comm , only : dump_field3d use boundary , only : use_o3du #ifdef with_budgets use budget_global,only : budg_dat,nzon_vg use chem_param , only : ra #endif implicit none ! --- in/out ----------------------------- integer, intent(in) :: region type(TDate), intent(in) :: tr(2) integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Chemie' ! real, parameter :: conv = 0.1*Avog/xm_o3 ! from kg/m2 --> #/cm2 real, parameter :: conv = 0.1*Avog/xmo3 ! from kg/m2 --> #/cm2 real, parameter :: todu = 3.767e-17 ! from #/cm2 --> du ! --- local ------------------------------- type(TDate) :: tmid integer :: imr, jmr, lmr integer :: i, j, l, ll, lll integer :: level, lglob #ifdef with_budgets integer :: nzone,nzone_v #endif real, pointer :: rm(:,:,:,:) real, pointer :: vo3(:,:,:) real, pointer :: ozone_klimtop(:) real, pointer :: phlb(:,:,:) real, pointer :: m(:,:,:) real, pointer :: temper(:,:,:) #ifdef with_zoom integer, pointer :: zoomed(:,:) #endif real, pointer :: area(:) ! cell area (m2) for each lat band real, allocatable :: Tclim(:,:) ! real,allocatable :: o3_overhead(:,:,:) ! in mlc/cm2 real,allocatable :: o3_overhead(:,:) ! in mlc/cm2 real :: lon, lat integer :: daynr, hour, minu, sec real :: csza real :: pf ! full level pressure (Pa) real :: dt_sec real :: o3_vmr ! ozone mixing ratio (ppv) real :: o3_ohc ! overhead column (mlc/cm2) real :: o3_ohc_fac real :: o3c_vmr ! ozone mixing ratio (ppv) real :: o3c_ohc ! overhead column (mlc/cm2) real :: kg_TO_mlc_m2 real :: Xpsc, dXpsc, kXpsc real :: do3_vmr, do3 real :: cc(nc) integer :: ilast real :: k_h real :: rm_new,A,B,sgn,timesc !type(TTimeSeriesHDF) :: F ! --- begin --------------------------------- if ( lmloc == 0 ) then status=0; return end if #ifdef MPI which_par = previous_par(region) if (which_par == 'tracers') & call escape_tm('Model should be parallel over levels in Cariolle chemistry') #endif call goLabel( rname ) if (myid == 0) then call wrtgol( 'ozone parchem : ', tr(1), ' - ', tr(2) ); call goPr if ( with_cold_tracer ) then write (gol,*) 'ozone parchem with cold tracer' ; call goPr else write (gol,*) 'ozone parchem without cold tracer' ; call goPr endif if ( with_trop_force ) then write (gol,*) 'ozone parchem with nudged tropospheric climatology' ; call goPr else write (gol,*) 'ozone parchem without nudged tropospheric climatology' ; call goPr endif end if ! regional grid sizes: imr = im(region) jmr = jm(region) lmr = lm(region) ! set pointers #ifdef MPI rm => mass_dat(region)%rm_k ! kg tracer, parallel over levels #else rm => mass_dat(region)%rm_t ! kg tracer, parallel over tracers #endif m => m_dat(region)%data ! kg air phlb => phlb_dat(region)%data ! air pressure at boundaries (1:lm+1) vo3 => phot_dat(region)%vo3 ! overhead O3, #/cm2, parallel over levels ozone_klimtop => phot_dat(region)%o3klim_top ! top layers overhead O3, #/cm2 temper => temper_dat(region)%data ! K area => region_dat(region)%dxyp ! m2 #ifdef with_zoom zoomed => region_dat(region)%zoomed #endif allocate(o3_overhead(jmr,lmr)) ; o3_overhead = 0.0 ! mid of interval: tmid = tr(1) + (tr(2)-tr(1))/2 ! ! ** setup climatology ! ! interpolate to this time, eventually read new fields: call Setup_O3Clim( lli(region), levi, O3clim_vmr(region), tmid, status ) IF_NOTOK_RETURN(status=1) if (use_o3du) then ! overhead ozone column for model top is calculated (in photolysis routine) ! based on Fortuin-Kelder climatology scaled by observations ! same should be done here for climatological overhead column do j=1,jmr o3_overhead(j,lmr) = ozone_klimtop(j) enddo else ! convert climatology from ppv to mlc/cm2 using standard airmass/m2 (p=1e5 Pa) !PRIOR ! o3_overhead(:,:,lmr) = 0.5 * O3clim_vmr(region)%data(:,:,lmr)*levi%m0(lmr) * conv / fscale(io3) !NOW o3_overhead(:,lmr) = 0.5 * O3clim_vmr(region)%data(1,:,lmr)*levi%m0(lmr) * conv / fscale(io3) endif do l = lmr-1,1,-1 ! convert climatology from ppv to mlc/cm2 using standard airmass/m2 (p=1e5 Pa) !PRIOR ! o3_overhead(:,:,l) = o3_overhead(:,:,l+1) + & ! 0.5 * ( O3clim_vmr(region)%data(:,:,l)*levi%m0(l) + & ! O3clim_vmr(region)%data(:,:,l+1)*levi%m0(l+1) ) * conv / fscale(io3) !NOW o3_overhead(:,l) = o3_overhead(:,l+1) + & 0.5 * ( O3clim_vmr(region)%data(1,:,l)*levi%m0(l) + & O3clim_vmr(region)%data(1,:,l+1)*levi%m0(l+1) ) * conv / fscale(io3) end do ! ! ** setup coefficients ! ! tables: ! bsf15k:/usr/people/eskes/fdscia/input/cariolle/OZONE_64x26.ASCII ! linoz/linoz_coeff.dat ! input routines: ! /usr/people/segers/work/tm/TM5/proj/chem/ozone.bak/src/chem_ozonepar.F90 ! check if current coeff are valid for this month, ! eventually (re)load: call Setup_Cariolle( lli(region), o3coef(region), o3coefp(region), tmid, status ) IF_NOTOK_RETURN(status=1) ! ! zonal average temperature ! ! compute smoothed temperature climatology ? if ( apply_taver ) then ! setup zonal field: allocate( Tclim(jmr,lmr) ) ! fill smoothed temperature: call Setup_Tclim( lli(region), temper, Tclim, status ) IF_NOTOK_RETURN(status=1) end if ! ! apply chemistry ! ! time step dt_sec = rTotal( tr(2) - tr(1), 'sec' ) ! loop over all grid cells: do j = jsr(region), jer(region) do i = isr(region), ier(region) ! skip if this cell is calculated in an overlying zoom region: #ifdef with_zoom if ( region_dat(region)%zoomed(i,j) /= region ) cycle #endif ! grid cell location: lon = lli(region)%lon_deg(i) ! deg lat = lli(region)%lat_deg(j) ! deg ! compute cos(solar_zenith_angle) daynr = DayNumber( tmid ) call Get( tmid, hour=hour, min=minu, sec=sec ) csza = cos_sza( daynr, hour, minu, sec, lon, lat ) ! conversion from (kg o3) to (mlc o3/m2) : ! ( kg O3 / m2 ) / (kg o3 / mol) * (mlc/mol) * (1e-4 m2/cm2) = mlc/cm2 kg_TO_mlc_m2 = 1.0 / area(j) / xm_o3 * Avog * 1.0e-4 ! loop over levels ilast = 1 do level = 1,lmloc lglob = offsetl + level ! Cariolle chemistry not applied below l_cariolle ! except when tropospheric nudging is active if ( (.not.with_trop_force) .and. (lglob .lt. l_cariolle(region)) ) exit ! full level pressure: pf = ( phlb(i,j,lglob) + phlb(i,j,lglob+1) )/2.0 ! Pa ! ozone volume mixing ratios: o3_vmr = fscale(io3) * rm(i,j,level,io3) / m(i,j,lglob) ! ppv o3_vmr = max(o3_vmr,0.0) ! can become negative... o3c_vmr = o3clim_vmr(region)%data(1,j,lglob) ! overhead column (including contribution of upper-half of current layer): o3_ohc = vo3(i,j,level) ! mlc/cm2 o3c_ohc = o3_overhead(j,lglob) ! mlc/cm2 o3_ohc_fac = 0.5*kg_TO_mlc_m2 / fscale(io3) * m(i,j,lglob) o3_ohc_fac = max(o3_ohc_fac,0.0) !if (i.eq.10.and.j.eq.10)then ! write(*,'a15,3i3,5e10.2')'DEBUG CARIOLLE',i,j,lglob,o3_ohc*todu,o3c_ohc*todu,o3_vmr*o3_ohc_fac*todu,o3_vmr,o3_ohc_fac !endif ! ** cold tracer scheme ! standard put values to zero, overwritten if with_cold_tracer = T Xpsc = 0. ! convert to mass again: dXpsc = 0. ! hetero chem (might be zero): k_h = 0. !kXpsc if ( with_cold_tracer ) then ! cold tracer (marked air mass) Xpsc = rm(i,j,level,ipsc) / m(i,j,lglob) ! ~ 0-1 Xpsc = min( max( 0.0, Xpsc ), 1.0 ) ! 0-1 ! hetrogeneous chemistry (cold tracer); ! fills kXpsc to be used as hetrochem reaction rate in Cariolle scheme: call ColdTracer_0d( Xpsc, pf, coldtracer_plevs, & temper(i,j,lglob), csza, lat, dt_sec, & dXpsc, kXpsc, status ) IF_NOTOK_RETURN(status=1) ! convert to mass again: dXpsc = dXpsc * m(i,j,lglob) ! kg ! add tracer mass changes: if (dXpsc .gt. 0.) then call AdjustTracer( dXpsc, region, i, j, level, ipsc, status ) IF_NOTOK_RETURN(status=1) end if #ifdef with_budgets nzone=budg_dat(region)%nzong(i,j) ! global budget nzone_v=nzon_vg(lglob) budcarg(nzone,nzone_v,2)=budcarg(nzone,nzone_v,2)+ dXpsc*1000./ra(ipsc) ! mol psc per cell #endif ! hetero chem (might be zero): k_h = kXpsc end if ! ** end applying cold tracer ! ** Cariolle scheme ! interpolate to pressure level: call Interp_Lin( o3coefp(region)%data(1,j,:), & o3coef(region)%data(:,j,:), pf, cc, ilast, status ) IF_NOTOK_RETURN(status=1) ! replace some coeffs by climatological data ! (have no impact if corresponding reaction rates are zero) cc(3) = o3c_vmr ! zonal aver ozone mixing ratio cc(7) = o3c_ohc ! zonal aver ozone overhead column ! use zonal aver temperature ? ! (has no impact if corresponding reaction rate is zero) if ( apply_taver ) cc(5) = Tclim(j,lglob) ! apply parameterized ozone chemistry for current cell: call Cariolle_0d( o3_vmr, o3_ohc, o3_ohc_fac, temper(i,j,lglob), & cc, k_h, cariolle_tau_min, dt_sec, do3_vmr ) IF_NOTOK_RETURN(status=1) ! if (i.eq.5.and. (j.eq.1 .or. j.eq.11 .or. j.eq. 80 .or. j.eq.90))then ! write(*,'a18,i3,i3,8e11.3')'DEBUG CARIOLLE cc',lglob,j,cc,dt_sec ! write(*,'a18,i3,i3,7e11.3')'DEBUG CARIOLLE o3',lglob,j,o3_vmr,o3_ohc,o3_ohc_fac,temper(i,j,lglob),pf,k_h,do3_vmr ! ! write(*,'a18,i3,i3,7e11.3')'DEBUG CAR timesc', lglob, j, 1./abs(cc(2) + cc(6)*o3_ohc_fac - k_h) ! endif ! ozone mass change: do3 = do3_vmr * m(i,j,lglob) / fscale(io3) ! kg ! add tracer mass changes: call AdjustTracer( do3, region, i, j, level, io3, status ) IF_NOTOK_RETURN(status=1) #ifdef with_budgets nzone=budg_dat(region)%nzong(i,j) ! global budget nzone_v=nzon_vg(lglob) budcarg(nzone,nzone_v,1)=budcarg(nzone,nzone_v,1)+do3*1000./ra(io3) ! mol o3 per cell #endif end do ! levels end do ! i end do ! j ! ! done ! nullify( rm ) nullify( m ) nullify( phlb ) nullify( vo3 ) nullify( ozone_klimtop ) nullify( temper ) #ifdef with_zoom nullify( zoomed ) #endif nullify( area ) #ifdef MPI rm => mass_dat(region)%rm_k ! kg tracer, parallel over levels #else rm => mass_dat(region)%rm_t ! kg tracer, parallel over tracers #endif m => m_dat(region)%data ! kg air nullify( rm ) nullify( m ) if ( apply_taver ) deallocate( Tclim ) deallocate (o3_overhead) ! ok call goLabel(); status=0 end subroutine Cariolle_Chemie ! ================================================================ subroutine Setup_Cariolle( lli, carcoef, carcoefp, t, status ) use Binas , only : p0 ! standard pressure (1e5 Pa) use GO , only : TDate, Get, NewDate use Num , only : Interp_Lin use Grid , only : TllGridInfo use Climat , only : TClimat, Set, Setup use file_cariolle, only : Cariolle_Info, CARIOLLE_READ ! --- in/out -------------------------------- type(TllGridInfo), intent(in) :: lli type(TClimat), intent(inout) :: carcoef ! 1-7 type(TClimat), intent(inout) :: carcoefp ! Pa type(TDate), intent(in) :: t integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Setup_Cariolle' ! --- local ------------------------------ integer :: year, month integer :: nlats, nlevs real, allocatable :: lats(:) real, allocatable :: pp(:,:) real, allocatable :: ppX(:,:,:) real, allocatable :: cc(:,:,:) real, allocatable :: ccX(:,:,:) type(TDate) :: tr_in(2) integer :: j, l, ic ! --- begin ------------------------------ ! try to setup for requested time; ! return status=-1 if not filled yet or filled for wrong time: call Setup( carcoef, t, status ) if (status==0) return ! data needs to be (re)setup ... ! grid sizes: call Cariolle_Info( status, nlats=nlats, nlevs=nlevs ) IF_NOTOK_RETURN(status=1) ! setup storage: allocate( lats(nlats) ) allocate( pp(nlats,nlevs) ) allocate( ppX(1,lli%nlat,nlevs) ) allocate( cc(nlats,nlevs,nc) ) allocate( ccX(nc,lli%nlat,nlevs) ) ! extract month: call Get( t, year=year, month=month ) ! load coeff for this month (ppmv) call Cariolle_Read( trim(cariolle_fname), month, lats, pp, cc, status ) IF_NOTOK_RETURN(status=1) ! interpolate pressure field to latitudes do l = 1, nlevs call Interp_Lin( lats, pp(:,l), lli%lat_deg, ppX(1,:,l), status ) IF_NOTOK_RETURN(status=1) end do ! spatial interpolation; ! loop over coeffs: do ic = 1, nc ! interpolate to latitudes do l = 1, nlevs call Interp_Lin( lats, cc(:,l,ic), lli%lat_deg, ccX(ic,:,l), status ) IF_NOTOK_RETURN(status=1) end do !! interpolate to standard pressures !do j = 1, lli%nlat ! ! Cariolle coefs from top->surf, thus no negation required: ! call Interp_Lin( ppX(j,:), ccX(j,:), levi%fp0, cc_in(ic,j,:) ) !end do end do ! coeff loop ! set time range for which this data is valid: ! o begin of current month: tr_in(1) = NewDate( year, month, 01, 00, 00 ) ! o begin of next month: month = month + 1 if ( month > 12 ) then year = year + 1 month = 1 end if tr_in(2) = NewDate( year, month, 01, 00, 00 ) ! store: call Set( carcoef, status, data=ccX, tr=tr_in ) IF_NOTOK_RETURN(status=1) call Set( carcoefp, status, data=ppX, tr=tr_in ) IF_NOTOK_RETURN(status=1) ! set coeff such that in troposphere only relaxation to climatology remains: if ( with_trop_force ) then call Setup_TropForce( carcoef, carcoefp, status ) IF_NOTOK_RETURN(status=1) end if ! clear storage: deallocate( lats ) deallocate( pp ) deallocate( ppX ) deallocate( cc ) deallocate( ccX ) ! ok status = 0 end subroutine Setup_Cariolle ! ================================================================ subroutine Setup_TropForce( carcoef, carcoefp, status ) ! --- in/out -------------------------------- type(TClimat), intent(inout) :: carcoef ! c1-c7 type(TClimat), intent(inout) :: carcoefp ! Pa integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Setup_TropForce' ! --- local ------------------------------ integer :: nlat, nlev integer :: j, l integer :: ilev0 real :: plev0, tau0 real :: plev, tau real :: cc(7) real :: alfa ! --- begin ------------------------------ ! grid size nlat = size(carcoefp%data,2) nlev = size(carcoefp%data,3) ! loop over latitudes do j = 1, nlat ! first level above troposphere: if ( carcoefp%data(1,j,1) > carcoefp%data(1,j,nlev) ) then ! upwards ilev0 = nlev ! model top do l = 1, nlev if ( carcoefp%data(1,j,l) <= cariolle_plev ) then ilev0 = l exit end if end do else ! downwards ilev0 = 1 ! model top do l = nlev, 1, -1 if ( carcoefp%data(1,j,l) <= cariolle_plev ) then ilev0 = l exit end if end do end if ! pressure in lowest cariolle level: plev0 = carcoefp%data(1,j,ilev0) ! Pa tau0 = - 1.0/carcoef%data(2,j,ilev0) ! sec ! loop over levels do l = 1, nlev ! current pressure and coeffs: cc = carcoef%data(:,j,l) plev = carcoefp%data(1,j,l) ! different pressure ranges: if ( plev >= tropo_plev ) then ! ~~ relaxation to climatology, slow ! relaxation towards climatology via Cariolle coeff: cc = 0.0 cc(2) = -1.0/tropo_tau else if ( plev >= cariolle_plev ) then ! ~~ relaxation to climatology, slower downwards ! interpolate tau between lower constant value and cariolle value: alfa = ( plev - plev0 ) / ( tropo_plev - plev0 ) tau = tau0 * (1.0-alfa ) + tropo_tau * alfa ! relaxation towards climatology via Cariolle coeff: cc = 0.0 cc(2) = -1.0/tau end if ! restore: carcoef%data(:,j,l) = cc end do ! levs end do ! lats ! ok status = 0 end subroutine Setup_TropForce ! ================================================================ subroutine Setup_O3clim( lli, levi, o3c, t, status ) use GO , only : TDate, Get, NewDate use Num , only : Interp_Lin use Grid , only : TllGridInfo, TLevelInfo use Climat , only : TClimat, Set, Setup use file_fortkeld, only : FortuinKelder_Info, FortuinKelder_Read ! --- in/out -------------------------------- type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi type(TClimat), intent(inout) :: o3c type(TDate), intent(in) :: t integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Setup_O3clim' ! --- local ------------------------------ integer :: year, month, day type(TDate) :: t_in(2) integer :: it integer :: nlats, nlevs real, allocatable :: lats_deg(:) real, allocatable :: pres_Pa(:) real, allocatable :: o3_ppv(:,:) real, allocatable :: o3X(:,:) real, allocatable :: sp(:) real, allocatable :: o3_in(:,:,:) integer :: j, l ! --- begin ------------------------------ ! try to interpolate in time to t; ! return status=-1 if not filled yet or filled for wrong time: call Setup( o3c, t, status ) if (status==0) return ! data needs to be (re)setup ... ! extract grid sizes: call FortuinKelder_Info( status, nlats=nlats, nlevs=nlevs ) IF_NOTOK_RETURN(status=1) ! setup arrays: allocate( lats_deg(nlats) ) allocate( pres_Pa(nlevs) ) allocate( o3_ppv(nlats,nlevs) ) allocate( o3X(lli%nlat,nlevs) ) allocate( sp(lli%nlat) ) allocate( o3_in(1,lli%nlat,levi%nlev) ) ! loop over interpolation times: do it = 1, 2 ! set times to 16th 00:00 call Get( t, year=year, month=month, day=day ) if ( it == 1 ) then ! previous mid of month if ( day <= 15 ) then ! mid of previous month month = month - 1 if ( month < 1 ) then year = year - 1 month = 12 end if end if else ! next mid of month if ( day >= 16 ) then ! mid of next month month = month + 1 if ( month > 12 ) then year = year + 1 month = 1 end if end if end if ! load climatology for this month (ppmv) call FortuinKelder_Read( trim(climat_fname), month, & lats_deg, pres_Pa, o3_ppv, status ) IF_NOTOK_RETURN(status=1) ! convert from ppmv to ppv: o3_ppv = o3_ppv * 1.0e-6 ! interpolate to latitudes do l = 1, nlevs call Interp_Lin( lats_deg, o3_ppv(:,l), lli%lat_deg, o3X(:,l), status ) IF_NOTOK_RETURN(status=1) end do ! interpolate to standard pressures do j = 1, lli%nlat ! negate pressure ax; should be increasing ... call Interp_Lin( -pres_Pa, o3X(j,:), -levi%fp0, o3_in(1,j,:), status ) IF_NOTOK_RETURN(status=1) end do ! time is mid of month: t_in(it) = NewDate( year, month, 16, 00, 00 ) ! store: call Set( o3c, status, t_in=t_in(it), data_in=o3_in, it_in=it ) IF_NOTOK_RETURN(status=1) end do ! times ! clear: deallocate( lats_deg ) deallocate( pres_Pa ) deallocate( o3_ppv ) deallocate( o3X ) deallocate( sp ) deallocate( o3_in ) ! set time range for which data could be interpolated: call Set( o3c, status, tr=t_in ) IF_NOTOK_RETURN(status=1) ! apply interpolation: should succeed now: call Setup( o3c, t, status ) if ( status < 0 ) then write (gol,'("time interpolation failed, while data has just been read.")'); call goErr write (gol,'("in ",a)') rname; call goErr; status=1; return else if ( status > 0 ) then write (gol,'("in ",a)') rname; call goErr; status=1; return end if ! ok status = 0 end subroutine Setup_O3clim ! ================================================================ subroutine Setup_Tclim( lli, T, Taver, status ) use Grid, only : TllGridInfo ! --- in/out -------------------------------- type(TllGridInfo), intent(in) :: lli real, intent(in) :: T(:,:,:) real, intent(out) :: Taver(:,:) integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Setup_Tclim' ! latitude range for model temperature average real, parameter :: latrange = 6.2 ! --- local ------------------------------ integer :: im, jm, lm integer :: i, j, l, k integer :: kmin, kmax real :: wk real :: taccu, waccu ! --- begin ------------------------------ im = size(T,1) jm = size(T,2) lm = size(T,3) do l = 1, lm do j = 1, jm taccu = 0.0 waccu = 0.0 ! cells within lat+[-latrange,+latrange] kmin = max(j-int(latrange/lli%dlat_deg),1) kmax = min(j+int(latrange/lli%dlat_deg),jm) ! loop over smoothing cells: do k = kmin, kmax ! weight decreasing with distance: wk = 1.0 - abs(real(k-j))*lli%dlat_deg/latrange if ( wk <= 0.0 ) then write (gol,'("internal error : wk = ",es12.4)') wk write (gol,'("in ",a)') rname; call goErr; status=1; return end if ! accumulate do i = 1, im taccu = taccu + wk * T(i,k,l) waccu = waccu + wk end do ! lon end do ! smoothing lats ! average: Taver(j,l) = taccu / waccu end do ! lats end do ! levs ! ok status = 0 end subroutine Setup_Tclim end module Chemistry_Cariolle