#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" ! ! ! Chemistry routines per cell. ! ! Parameterized ozone scheme. ! module chemistry_0d implicit none ! --- in/out -------------------------------------- private public :: Cariolle_0d !public :: TropForce_0d public :: ColdTracer_0d ! --- const ----------------------------- character(len=*), parameter :: mname = 'chemistry_0d' contains ! ======================================================= ! :::::::::::::::::: CARIOLE CHEMISTRY SCHEME :::::::::::::::::: ! ! pseudo chemistry with Cariole parametrization (see also notes) ! ! dr/dt = c1 + c2*( r - c3 ) + c4*( T - c5 ) + c6*( Splus - c7 ) - k_h*r ! ! = c1 + c2*( r - c3 ) + c4*( T - c5 ) + c6*( Splus-f*r + f*r - c7 ) - k_h*r ! ! = [c1 - c2*c3 + c4*( T - c5 ) + c6*( Splus-f*r - c7 ) ] + [ c2 + c6*f - k_h ]*r ! ! = A + B*r = B ( r - -A/B ) ! = B ( r - C ) with C = -A/B ! ! with ! A = c1 - c2 c3 + c4 delT + c6(Splus-c7), ! B = c2 + term accounting for contribution of the current ! layer to the column above this layer ! - extra term k_h for hetero chem ! r = volume mixing ratio ! delT = temp difference TM3/Cariolle in K ! Splus = ozone above current layer (molec cm-2) ! f = fraction for which r contributes to Splus ! c1..c7 = Cariolle coefficients, t = time in seconds ! k_h = extra term, for example from cold tracer scheme ! ! local solutions: ! ! r(t+dt) = -A/B + ( r(t) - -A/B ) exp[ B dt ] , B /= 0 ! = C + ( r(t) - C ) exp[ B dt ] ! ! r(t+dt) = r(t) + A dt , B == 0 ! ! Henk Eskes, KNMI, May 1999 ! ! o Scheme implemented for use in TM5 ! Above top of the model prescribed ozone column from cariolle climatology ! Bram Bregman Wed Oct 30 11:41:09 MET 2002 ! ! o Generalized contribution of current layer to overhead column. ! Arjo Segers, september 2005. ! ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine Cariolle_0d( o3, o3o, f, temper, cc, k_h, & timesc_min, dt, do3 ) use GO , only : gol, goPr ! --- in/out ------------------------------ real, intent(in) :: o3 ! ozone volume mixing ratio (ppv) real, intent(in) :: o3o ! ozone overhead (mlc/m2) real, intent(in) :: f ! factor for which o3 contributes to o3o real, intent(in) :: temper ! temperature (K) real, intent(in) :: cc(7) ! Cariolle coeff real, intent(in) :: k_h ! extra reaction rate hetro chem real, intent(in) :: timesc_min ! minimum time scale sec : [0,..) real, intent(in) :: dt ! sec real, intent(out) :: do3 ! ozone increment ! --- local ------------------------------- real :: A, B real :: sgn, timesc real :: o3_new ! --- begin ------------------------------- ! calculate the A term, including extra contribution from dd : A = cc(1) - cc(2) * cc(3) + cc(4) * ( temper - cc(5) ) + cc(6) * ( o3o-f*o3 - cc(7) ) ! calculate the B term, including the ozone column fraction and extra term: B = cc(2) + cc(6)*f - k_h ! different solutions for B==0 and B/=0 : if ( abs(B) < tiny(1.0) ) then ! ! local solution: ! ! r(t+dt) = r(t) + A dt ! o3_new = o3 + A * dt else ! ! local solution: ( c2 < 0, and normally A > 0 and B < 0) ! ! r(t+dt) = ( r(t) + A/B ) exp[ B dt ] - A/B ! ! time scale of reaction: ! exp[ B dt ] = exp[ sgn |B| dt ] = exp[ sgn dt / (1/|B|) ] ! bound time scale to prevent too fast chemistry: sgn = sign( 1.0, B ) timesc = max( timesc_min, 1.0/abs(B) ) ! sec, [0,..) ! apply: o3_new = ( o3 + A/B )*exp(sgn*dt/timesc) - A/B end if ! trap negatives ... o3_new = max( 0.0, o3_new ) ! return increment: do3 = o3_new - o3 end subroutine Cariolle_0d ! ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ! ! ! This subroutine relaxes the tropospheric ozone to the ! ! climatology of Fortuin ! ! ! ! dr/dt = -1/tau ( r - r_clim ) ! ! ~= c2 ( r - c3 ) ! ! ! ! r(t+dt) = r_clim + ( r(t) - r_clim ) exp[ - dt / tau ] ! ! ! ! Relaxation time scale is altitude depended: ! ! ! ! tau(1) tau(2) ! ! (weeks) (months) ! ! | ! ! | * ------ plevs(1) (230 hPa) ! ! | o ! ! | ! ! | o ! ! | ! ! | o ! ! | * ------ plevs(2) (500 hPa) ! ! | o ! ! | ! ! | o ! ! | ! ! | o ! ! | ! ! | o ! ! --------------------------------> tau ! ! ! ! Henk Eskes, KNMI, May 2000 ! ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ! subroutine TropForce_0d( o3, o3c, p, plevs, taus, dt, do3, status ) ! ! ! --- in/out --------------------------------- ! ! real, intent(in) :: o3 ! ozone ! real, intent(in) :: o3c ! ozone climat ! real, intent(in) :: p ! pressure level (Pa) ! real, intent(in) :: dt ! sec ! real, intent(in) :: plevs(2) ! pressure levels (Pa) ! real, intent(in) :: taus(2) ! time scales (sec) ! real, intent(out) :: do3 ! ozone increment ! integer, intent(out) :: status ! ! ! --- const -------------------------------- ! ! character(len=*), parameter :: rname = mname//'/TropForce' ! ! ! --- local --------------------------------- ! ! real :: tau ! real :: alfa ! real :: o3_new ! ! ! --- begin ------------------------------------ ! ! ! only for troposphere: ! if ( p < plevs(1) ) then ! ! ! just copy: ! o3_new = o3 ! ! else ! ! ! troposphere; relax towards climat: ! ! ! set time scale: ! if ( p >= plevs(2) ) then ! ! lower troposphere ! tau = taus(1) ! else ! ! interpolate between plev(1) and plev(2) : ! alfa = ( p - plevs(1) ) / ( plevs(2) - plevs(1) ) ! tau = taus(1) * (1.0-alfa ) + taus(2) * alfa ! end if ! ! ! relax: ! o3_new = o3c + ( o3 - o3c ) * exp( - dt / tau ) ! ! end if ! ! ! return increment: ! do3 = o3_new - o3 ! ! ! ok ! status = 0 ! ! end subroutine TropForce_0d ! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ! Parameterization of PSC's ! ! A cold tracer X is used to mark and track air below a temperature threshold. ! ! Production and loss of X : ! ! dX/dt = k_p ( 1 - X ) - k_l X ! ! k_p = 1/tau_p : T < Tpsc and |lat| > 40 deg ! = 0 : otherwise ! ! k_l = 1/tau_l : sza > 0 (day) ! = 0 : otherwise (night) ! ! Rewrite: ! ! dX/dt = k_p - k_p X - k_l X = -(k_p+k_l) [ X - k_p/(k_p+k_l) ] ! ! with local solution: ! ! X(t+dt) = X(t) , if k_p = k_l = 0 ! ! = kp/(kp+kl) + [ X(t) - kp/(kp+kl) ] exp[ -(kp+kl) dt ] , otherwise ! ! Adding an additional term to the Cariolle Parameterization: ! ! dO3/dt = c1 + c2*(O3-c3) + c4*(T-c5) + c6*(S-c7) - kX O3 ! ^^^^^^^ ! where ! ! k = 1/tau : sza > 0 (day) ! = 0 : otherwise ! ! Peter Braesicke, Cambridge, Centre for Atmospheric Science ! ! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ! Routine returns: ! o new cold tracer value : X(t+dt) ! o factor to be used in Cariolle scheme : kX(t) ! ! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine ColdTracer_0d( X, p, act_plevs, T, csza, lat, dt, dX, kX, status ) use go, only : gol, goErr, goPr use go, only : TDate ! --- in/out ---------------------------- real, intent(in) :: X ! X(t) (0-1) real, intent(in) :: p ! pressure (Pa) real, intent(in) :: act_plevs(2) ! active pressure range real, intent(in) :: T ! temperature (K) real, intent(in) :: csza ! cos_sza [-1,1] real, intent(in) :: lat ! latitude (deg) real, intent(in) :: dt ! time step (sec) real, intent(out) :: dX ! X(t+dt)-X(t) real, intent(out) :: kX ! to be used in cariolle scheme integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/ColdTracer_0d' ! --- const ------------------------------ ! time scale for heterogeneous breakdown of ozone - 10 to 20 days real, parameter :: tau_het = 12.0 * 86400.0 ! 12 days ! time scale for creating activation real, parameter :: tau_prod = 4.0 * 3600.0 ! 4 hours ! time scale for loss of activation, ! Northern and Southern hemisphere - 5,10 d real, parameter :: tau_loss_nh = 5.0 * 86400.0 ! 5 days real, parameter :: tau_loss_sh = 10.0 * 86400.0 ! 10 days ! --- local -------------------------------- logical :: active_level real :: Tpsc real :: kp, kl real :: X_new ! --- begin ------------------------------------ ! check ... if ( (X < 0.0) .or. (X > 1.0+1e-3) ) then write (gol,'("cold tracer X should be in [0,1] ; found : ",es12.4)') X; call goErr write (gol,'("in ",a)') rname; call goErr; status=1; return end if ! active level ? active_level = (p >= minval(act_plevs)) .and. (p <= maxval(act_plevs)) ! ! hetro chem : ozone loss by psc ! ! hetero chem ? ! set factor in Cariolle scheme if day and cold tracer > 0.03 if ( (csza > 0.0) .and. (X > 0.03) .and. active_level ) then ! day kX = X / tau_het else ! night or not in strato: kX = 0.0 end if ! ! production and loss of activation field ! ! Compute Temp_NAT as function of pressure . ! Document 'Parameterization of PSCs' by Peter Braesicke ! gives table of pressures (hPa) and temperature (K) : ! ! pressure (hPa) temperature (K) ! -------------- --------------- ! 31.62 191.74 ! 38.31 192.77 ! 46.42 193.81 ! 56.23 194.85 ! 68.13 195.91 ! 82.54 196.98 ! 100.00 198.06 ! 120.23 199.12 ! 141.90 200.07 ! 164.40 200.93 ! 187.42 201.69 ! 210.43 202.38 ! 233.43 202.99 ! ! This points are on a line, thus probably created from a parameterization. ! Linear fit by idl through the points: ! temperature_K = 146.229 + 5.635 * LOG( Pressure_Pa ) ! Tpsc = 146.229 + 5.635*log(p) ! production rate: if ( (T < Tpsc) .and. (abs(lat) > 40.0) .and. active_level ) then kp = 1.0 / tau_prod else kp = 0.0 end if ! loss rate: if ( csza > 0.0 ) then ! daylight, thus destruction of psc if ( lat > 0.0 ) then ! northern hemisphere kl = 1.0 / tau_loss_nh else ! southern hemisphere kl = 1.0 / tau_loss_sh end if else ! night kl = 0.0 end if ! apply destruction and loss: if ( (kp < tiny(1.0)) .and. (kl < tiny(1.0)) ) then ! no production and loss, thus unchanged: X_new = X else ! local solution: X_new = kp/(kp+kl) + ( X - kp/(kp+kl) ) * exp( -(kp+kl)*dt ) end if ! truncate: X_new = min( max( 0.0, X_new ), 1.0 ) ! return increment: dX = X_new - X ! ! done ! ! ok status = 0 end subroutine ColdTracer_0d end module chemistry_0d