Skip to content
Snippets Groups Projects
Normal view Permalink
ch_exqssa.f90 8.28 KiB
Newer Older
  • Learn to ignore specific revisions
    • ESCOBAR MUNOZ Juan's avatar
    Juan 9/01/2014: typo on MNH_LIC : 2014 & LICENSE  
    ESCOBAR MUNOZ Juan committed
    1 
    !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
    ESCOBAR MUNOZ Juan's avatar
    Juan 9/01/2014: add header MNH_LIC to all MNH/*.f* files  
    ESCOBAR MUNOZ Juan committed
    2 
    !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
    ESCOBAR MUNOZ Juan's avatar
    Juan 9/01/2014: typo on MNH_LIC : 2014 & LICENSE  
    ESCOBAR MUNOZ Juan committed
    3 
    !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
    ESCOBAR MUNOZ Juan's avatar
    Juan 9/01/2014: add header MNH_LIC to all MNH/*.f* files  
    ESCOBAR MUNOZ Juan committed
    4 
    !MNH_LIC for details. version 1.
    WAUTELET Philippe's avatar
    Beginning of open source history
    WAUTELET Philippe committed
    5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 
    !-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
! MASDEV4_7 chimie 2006/05/18 13:07:25
!-----------------------------------------------------------------
!!    ##################### 
    MODULE MODI_CH_EXQSSA
!!    ##################### 
!!
!
INTERFACE
!!
SUBROUTINE CH_EXQSSA(PTSIMUL, PDTACT, PCONC, PNEWCONC, KEQ, KVECNPT, KMI, &
                       PSLOW, PFAST, PDTMAX)
IMPLICIT NONE
REAL,    INTENT(IN) :: PTSIMUL  ! time of simulation
REAL,    INTENT(IN) :: PDTACT   ! actual time-step
INTEGER, INTENT(IN) :: KEQ      ! dimension of the problem to solve
INTEGER, INTENT(IN) :: KVECNPT
REAL, INTENT(IN),  DIMENSION(KVECNPT,KEQ) :: PCONC    
                              ! concentration vector at PTSIMUL
REAL, INTENT(OUT), DIMENSION(KVECNPT,KEQ) :: PNEWCONC 
                              ! solution at PTSIMUL + PDTACT
INTEGER, INTENT(IN) :: KMI    ! model number
                              ! reac. rates, auxiliary variables
REAL,    INTENT(IN) :: PSLOW, PFAST, PDTMAX
END SUBROUTINE CH_EXQSSA
!!
END INTERFACE
!!
END MODULE MODI_CH_EXQSSA
!!    ######################################################################### 
SUBROUTINE CH_EXQSSA(PTSIMUL, PDTACT, PCONC, PNEWCONC, KEQ, KVECNPT, KMI, &
                       PSLOW, PFAST, PDTMAX)
!!    #########################################################################
!!
!!*** *CH_EXQSSA*
!!
!!    PURPOSE
!!    -------
!    solve one time-step of the chemical differential equation d/dt C = P - L*C
!!
!!**  METHOD
!!    ------
!!    Quasy steady state approximation (EXQSSA):
!!      C_n+1 = Q(C_n,h) = P_n / L_n + ( C_n - P_n / L_n ) exp( - L_n Dt)
!!
!!    Extrapolation
!!    C_n+1(h) = Q(C_n,h)
!!    C_n+0.5(h/2)=Q(C_n,h/2)
!!    C_n+1(h/2)=Q(C_n+0.5,h/2)
!!    C_n+1 = 2 C_n+1(h/2) - C_n+1(h)
!!
!!    REFERENCE
!!    ---------
!!    O. Hertel et al ,
!!    Test of two Numerical schemes for use in atmospheric transport 
!!    chemistry models,
!!    Atmospheric environment, 16, 2591--2611, 1993.
!!
!!    A. Sandu et al
!!    Benchmarking Stiff Ode Solvers For Atmospheric Chemistry Problems I:
!!    Implicit versus Explicit
!!    submitted to Atmospheric Environment 26/01/1996
!!
!!    AUTHOR
!!    ------
!!    A. FASSI FIHRI    *Laboratoire d'Aerologie*
!!
!!    MODIFICATIONS
!!    -------------
!!    Original 13/02/96
!!    10/04/96 (K. Suhre) take out performance control statements
!!    31/07/96 (K. Suhre) restructured
!!    05/03/98 (V. Crassier) vectorized with fixed timestep
!!    31/03/99 (K. Suhre) change minor details (i.e. numerical pbs with
!!                        LOSS=0 species under Absoft and HP f90 compilers
!!                        and control of correct time stepping if user input
!!                        is uncorrect) and restructure
!!    01/08/01 (C. Mari)   add PSLOW, PFAST, PDTMAX arguments
!!    01/12/03 (D. Gazen)  change Chemical scheme interface
!!
!!    EXTERNAL
!!    --------
USE MODI_CH_PRODLOSS
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!------------------------------------------------------------------------------
!
!*       0.   DECLARATIONS
!        -----------------
IMPLICIT NONE
!
!*       0.1  declaration of arguments
!
REAL,    INTENT(IN) :: PTSIMUL  ! time of simulation
REAL,    INTENT(IN) :: PDTACT   ! actual time-step
INTEGER, INTENT(IN) :: KEQ      ! dimension of the problem to solve
INTEGER, INTENT(IN) :: KVECNPT
REAL, INTENT(IN),  DIMENSION(KVECNPT,KEQ) :: PCONC    
                              ! concentration vector at PTSIMUL
REAL, INTENT(OUT), DIMENSION(KVECNPT,KEQ) :: PNEWCONC 
                              ! solution at PTSIMUL + PDTACT
INTEGER, INTENT(IN) :: KMI    ! model number
                              ! reac. rates, auxiliary variables
REAL,    INTENT(IN) :: PSLOW, PFAST, PDTMAX
!
!*       0.2  declaration of local variables
!
REAL, DIMENSION(KVECNPT,KEQ)     :: ZCONC, ZPROD, ZLOSS
REAL, DIMENSION(KVECNPT,KEQ)     :: ZCONC1, ZCONC2
REAL, DIMENSION(KVECNPT,KEQ)     :: ZCONC1A, ZCONC2A
REAL, DIMENSION(KVECNPT,KEQ)     :: ZCONC1B, ZCONC2B
REAL, DIMENSION(KVECNPT,KEQ)     :: ZCONC1C, ZCONC2C
REAL, DIMENSION(KVECNPT,KEQ)     :: ZTMP, ZRATIO
REAL                             :: ZTIME, ZSTEP
INTEGER                          :: JI
REAL, DIMENSION(KVECNPT,KEQ)     :: ZA, ZB, ZC
REAL, PARAMETER                  :: ZEPS = 1E-14 ! to avoid division by zero
INTEGER                          :: IITER  ! number of internal timesteps
!
!------------------------------------------------------------------------------
!
!*       1.   INITIALIZE SOME VARIABLES
!        -------------------------------
!
ZTIME  = PTSIMUL 
ZCONC  = PCONC
IITER  = MAX(1,INT(1.E-3+(PDTACT/PDTMAX)))
ZSTEP  = PDTACT / IITER
!
!*       2.   INTERATION LOOP
!        --------------------
!
DO JI=1, IITER
!
!*       2.1  calculate production and loss terms
!
  CALL CH_PRODLOSS(ZTIME,ZCONC,ZPROD,ZLOSS,KMI,KVECNPT,KEQ)
!
! elimiate negative loss rates that may be inherited from the transport scheme
!
  ZLOSS(:,:) = MAX(0.,ZLOSS(:,:))
  ZA(:,:) = 0.
  ZB(:,:) = 0.
  WHERE ( ZLOSS .LE. 1./(PSLOW*ZSTEP) ) ZA(:,:) = 1.
  WHERE ( ZLOSS .GT. 1./(PFAST*ZSTEP) ) ZB(:,:) = 1.
  ZC(:,:) = 1. - ZA(:,:) - ZB(:,:)
!        
!*       2.2  first guess
!
! especes lentes ZLOSS<1./(PSLOW*ZSTEP)
! calcul des constantes multiplicatives
!***************************************

  ZTMP(:,:)    = ZSTEP*(ZPROD(:,:) - ZLOSS(:,:)*ZCONC(:,:)) 
  ZCONC1A(:,:) = ZCONC(:,:) + ZTMP(:,:) 
  ZCONC2A(:,:) = ZCONC(:,:) + 0.5 * ZTMP(:,:)
  
! especes rapides ZLOSS>1./(PFAST*ZSTEP)
! calcul des constantes multiplicatives
!***************************************
  
  ZCONC1B(:,:) = ZPROD(:,:) / MAX(ZEPS, ZLOSS(:,:))
  ZCONC2B(:,:) = ZCONC1B(:,:)
  
! especes intermediaires
! calcul des constantes multiplicatives
!***************************************

  ZRATIO(:,:)  = ZPROD(:,:) / MAX(ZEPS, ZLOSS(:,:))
  ZTMP(:,:)    = EXP(-ZLOSS(:,:)*ZSTEP*0.5)
  ZCONC1C(:,:) = ZTMP(:,:) * ZTMP(:,:) * ( ZCONC(:,:) - ZRATIO(:,:) ) &
                 + ZRATIO(:,:)
  ZCONC2C(:,:) = ZTMP(:,:) * ( ZCONC(:,:) - ZRATIO(:,:) ) + ZRATIO(:,:)

! calcul des concentrations intermediaires
!*****************************************

  ZCONC1(:,:)= ZA(:,:)*ZCONC1A(:,:) &
             + ZB(:,:)*ZCONC1B(:,:) &
             + ZC(:,:)*ZCONC1C(:,:)
  ZCONC2(:,:)= ZA(:,:)*ZCONC2A(:,:) &
             + ZB(:,:)*ZCONC2B(:,:) &
             + ZC(:,:)*ZCONC2C(:,:)

!
!*      2.3  calculate new production and loss terms
!
  CALL CH_PRODLOSS(ZTIME+0.5*ZSTEP,ZCONC2,ZPROD,ZLOSS,KMI,KVECNPT,KEQ)
!
! elimiate negative loss rates that may be inherited from the transport scheme
!
  ZLOSS(:,:) = MAX(0.,ZLOSS(:,:))
  ZA(:,:) = 0.
  ZB(:,:) = 0.
  WHERE ( ZLOSS .LE. 1./(PSLOW*ZSTEP) ) ZA(:,:) = 1.
  WHERE ( ZLOSS .GT. 1./(PFAST*ZSTEP) ) ZB(:,:) = 1.
  ZC(:,:) = 1. - ZA(:,:) - ZB(:,:)
!
! especes lentes ZLOSS<1./(PSLOW*ZSTEP)
! calcul des constantes multiplicatives
!***************************************
  
  ZTMP(:,:) = ZSTEP*(ZPROD(:,:) - ZLOSS(:,:)*ZCONC(:,:)) 
  ZCONC2A(:,:) = ZCONC2(:,:) + 0.5 * ZTMP(:,:)
!
! especes rapides ZLOSS>1./(PFAST*ZSTEP)
! calcul des constantes multiplicatives
!***************************************

  ZCONC2B(:,:) = ZPROD(:,:) / MAX(ZEPS, ZLOSS(:,:))
  
!
! especes intermediaires
! calcul des constantes multiplicatives
!***************************************
  
  ZRATIO(:,:)  = ZPROD(:,:) /  MAX(ZEPS, ZLOSS(:,:))
  ZTMP(:,:)    = EXP(-ZLOSS(:,:)*ZSTEP*0.5)
  ZCONC2C(:,:) = ZTMP(:,:) * ( ZCONC2(:,:) - ZRATIO(:,:) ) + ZRATIO(:,:) 
  
! calcul des concentrations intermediaires
!*****************************************

  ZCONC2(:,:) = ZA(:,:)*ZCONC2A(:,:) &
              + ZB(:,:)*ZCONC2B(:,:) &
              + ZC(:,:)*ZCONC2C(:,:)
  
! IF ZSTEP*ZLOSS < 1. --> A=1 else B=1
!*************************************

  ZA(:,:) = 0.
  WHERE ( (ZSTEP*ZLOSS) .LT. 1. ) ZA(:,:) = 1.
  ZB(:,:) = 1. - ZA(:,:)

  ZCONC(:,:) = ZA(:,:)*MAX(0.,2.*ZCONC2(:,:)-ZCONC1(:,:)) &
             + ZB(:,:)*MAX(0.,ZCONC2(:,:))

  ZTIME=ZTIME+ZSTEP

ENDDO
!
PNEWCONC = ZCONC
!     
END SUBROUTINE CH_EXQSSA

    Documentation | Mentions légales | CGU | Accessibilité : non conforme