ecearth3/sources/tm5mp/proj/cb05/chemistry_0d__o3tracer.F90

#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