View source for Moisture Updates in v11-01

= Problem Description =

GCST member [[User:Lizzie_Lundgren|Lizzie Lundgren]] (Harvard) discovered that tracer unit conversions between [v/v] and [kg] do not accurately handle moisture by using moist air mass with dry air molecular weight. Updates to fix this issue and other moisture-related discrepancies discovered during the update implementation were introduced in [[GEOS-Chem_v11-01#v11-01a|GEOS-Chem v11-01a]].

Validation of the v11-01a moisture corrections revealed an additional moisture handling issue in advection. Inert tracer dry mixing ratio maps contained patterns resembling the moisture signature in the atmosphere. This phenomenon is due to lack of moisture handling within transport. Tracer mass is distributed relative to moist rather than dry air pressure and therefore preferentially pools within grid boxes that contain more water vapor mass.

Meemong Lee and Richard Weidner (JPL) issued a [http://wiki.seas.harvard.edu/geos-chem/images/6/6f/GCA-M2O2-Final.pdf JPL Publication] detailing a transport moisture fix for the Adjoint model. Their adjoint fix consists of (1) deriving dry surface pressures from GMAO moist surface pressures and specific humidity, and (2) replacing moist pressures with dry pressures throughout the model. 

Lizzie Lundgren (GCST) adapted JPL’s Adjoint fix to the GEOS-Chem forward model in [[GEOS-Chem_v11-01#v11-01h|GEOS-Chem v11-01h]]. The GEOS-Chem implementation (1) corrects the transport mass distribution problem by applying JPL’s fix to transport and unit conversions, (2) preserves the moisture updates incorporated into v11-01a, and (3) preserves inert tracer mass conservation.

= Summary of Changes =

Prior to GEOS-Chem v11-01, tracer unit conversions between [v/v] and [kg] did not accurately take moisture into account by combining total air mass with dry air molecular weight. The air quantity updates remedy this problem and fix some additional moisture-related bugs discovered during implemention. The updates were primarily introduced in v11-01a with some additional changes in v11-01h. Related changes include updating tracer units, discussed separately at [[GEOS-Chem_species_units|our ''GEOS-Chem species units'' wiki page]]. 

Below is a high-level summary of the moisture updates added to v11-01. All of these updates were incorporated in GEOS-Chem v11-01a, except where noted otherwise:

*'''Dry air partial pressure replaces wet atmospheric pressure for:'''
#air number density unit conversions
#gas partial pressure calculations

*'''Units of tracer concentration (State_Chm%Tracers), previously in [v/v], are now:'''
#[kg/kg dry air] in convection and advection
#[v/v dry air] everywhere else

*'''New State_Met variables:'''
#PMID_DRY - dry air partial pressure at pressure mid-point [hPa]
#PEDGE_DRY - dry air partial pressure at level lower edge [hPa]
#MAIRDEN - moist air density (previously AIRDEN) [kg/m3]
#AIRNUMDEN - dry air density [molec/cm3]
#TV - virtual temperature [K]

*'''Pre-existing State_Met variables that have modified definitions and algorithms:'''
#AD - now defined as dry air mass (previously total air mass), and derived from surface dry pressure
#AIRDEN - now defined as dry air density (previously total air density, including H2O), and derived from surface dry pressure

*'''Pre-existing State_Met variables that have unchanged definitions but modified algorithms:'''
#RH - now uses Nordquist 1973 formulation for water vapor saturation pressure
#BXHEIGHT - now uses virtual temperature in the hypsometric equation
#AVGW - now defines specific humidity as kg H2O / kg total air in calculation of mol H2O / mol dry air

*'''Pre-existing State_Met variables that have modified values due to changes elsewhere:'''
#AIRVOL - impacted by change in box height calculation

*'''Moist air mass and density replace dry air mass and density when:'''
#MOISTQ (tendency in specific humidity [kg/kg/s]) is used with air quantities for unit conversion
#SPHU (specific humidity [g/kg]) is used with air quantities for unit conversion

NOTE: All State_Met variables listed above are calculated within routine AIRQNT in DAO_Mod. Relative humidity (RH) and specific humdiity (SPHU) are also calculated within routine SET_H2O_TRAC in UCX_Mod under certain conditions.

= Pressure and Humidity =

== Pressure Uses - Dry vs. Moist ==

GEOS-Chem now uses both wet atmospheric pressure and dry air partial pressure. Previously wet atmospheric pressure was used for all pressure applications. Below is a summary of how each pressure is used following initial air quantity fix updates.

'''Dry air partial pressure is now used for:'''
#Air number density conversions (ucx_mod, mercury_mod, gasconc)
#Partial pressure calculations (ucx_mod, sulfate_mod)
#Surface pressure diagnostic ND31 (per request of Alex Turner for CH4
#Where usage is unclear (ruralbox - this routine will soon be obsolete after Flexchem is implemented)

'''Wet atmospheric pressure is used for all other pressure uses:'''
#Pressure level diagnostic outputs
#Cloud top height calculation
#Potential temperature calculation
#Volume cloud fraction calculation
#Pressure level thresholding (i.e. height proxy)
#Barometric formula
#Hypsometric equation (now uses virtual temperature in stat_chem_mod and vdiff_mod)
#Mean free path calculations

There is also a dry surface pressure introduced in v11-01h to correct a [[GEOS-Chem_v11-01#Incorrect_incorporation_of_moisture_in_GEOS-Chem_air_quantities|transport issue]] identified following moisture updates in v11-01a. The dry surface pressure is a bug fix submitted by Meemong Lee and Richard Weidner of JPL. It is used in advection and to construct mass per level for use in unit conversions, calculation of dry air mass, and calculation of dry air density. See the [[Moisture_Updates_in_v11-01#Transport|Transport section]] of this wiki page for more information.

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 15:54, 18 May 2015 (UTC)

== Moist Air Pressure ==

'''<big>Legacy implementation</big>'''

State_Met%PMID is allocated and initialized to 0 but it is not updated or used. Instead, the mid-point pressure is retrieved via the GET_PCENTER subroutine in GeosUtil/pressure_mod.F.

State_Met%PEDGE is used only if EXTERNAL_GRID or EXTERNAL_FORCING are defined and is not updated internally.  Instead, the grid box lower edge pressure is retrieved via the GET_PEDGE subroutine implemented in GeosUtil/pressure_mod.F.

         ! Pressure at bottom edge of grid box [hPa]
         P1                        = GET_PEDGE( I, J, L )
 
         ! Pressure at top edge of grid box [hPa]
         P2                        = GET_PEDGE( I, J, L+1 )
 
         ! Pressure difference between top & bottom edges [hPa]
         State_Met%DELP(I,J,L)     = P1 - P2

'''<big>New implementation</big>'''

There are now several moist air pressures to choose from:
#PEDGE - wet air pressure at bottom edge of grid box 
#PMID - wet air pressure at grid box pressure mid-point, defined as arithmetic mean of edge pressures
#PMEAN - locally defined wet air pressure at grid box altitude-weighted mid-point, defined as vertically-integrated mean pressure

DELP will represent the difference between upper and lower grid box edge moist air pressures. This is the same as the difference between upper and lower grid box edge dry air partial pressures when assuming constant humidity across the grid box. Humidity is not spatially interpolated to grid box edges in GEOS-Chem.

DELP, PMID, PMID_DRY, PEDGE, and PEDGE_DRY will all be set in AIRQNT in dao_mod.F. GET_PCENTER and GET_PEDGE calls elsewhere in GEOS-Chem are replaced with State_Met%PMID_DRY, State_Met%PMID, State_Met%PEDGE, and State_Met%PEDGE_DRY where appropriate. 

Parameters in GeosCore/dao_mod.F:
    
 ! Conversion factors
 REAL(fp), PARAMETER   :: TOFFSET = -2.7315e+2_fp   ! K -&gt; deg C    
 REAL(fp), PARAMETER   :: PCONV   = 1.00e+2_fp      ! hPa -&gt; Pa      
 REAL(fp), PARAMETER   :: RHCONV  = 1.00e-2_fp      ! % -&gt; fraction     
 REAL(fp), PARAMETER   :: MCONV   = 1.00e-3_fp      ! g -&gt; kg

Implementation in GeosCore/dao_mod.F:      

         !=============================================================
         ! Set wet air pressures [hPa] 
         ! (lower edge, delta, centroid, and spatially-weighted mean)
         !=============================================================  
 
         ! Pressure at bottom edge of grid box [hPa]
         State_Met%PEDGE( I,J,L )  = GET_PEDGE( I,J,L )
 
         ! Pressure at top edge of grid box [hPa]
         PEdge_Top                 = GET_PEDGE( I,J,L+1 )
 
         ! Pressure at bottom edge of grid box [hPa] (level LLPAR+1 only)
         IF (L == LLPAR) THEN
            State_Met%PEDGE( I,J,L+1 ) = PEdge_Top
         ENDIF
 
         ! Pressure difference between top & bottom edges [hPa]
         State_Met%DELP( I,J,L )   = State_Met%PEDGE( I,J,L ) 
     &                                - PEdge_Top
 
         ! Arithmetic average pressure in grid box [hPa] defined as
         ! ( PEDGE(L) + PEDGE(L+1) ) / 2. This represents the grid box 
         ! mass centroid pressure. Use in the ideal gas law yields
         ! a local air density at the centroid.
         State_Met%PMID( I,J,L )   = GET_PCENTER( I,J,L )
 
         ! Mean altitude-weighted pressure in grid box [hPa] defined as 
         ! average P(z) over z. Use in the ideal gas law yields a
         ! mean density equivalent to total mass per volume in grid box.
         PMEAN = State_Met%DELP(I,J,L) / log( State_Met%PEDGE(I,J,L) /
     &                                   PEdge_Top ) 

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:43, 29 September 2016 (UTC)

== Water Vapor Saturation Pressure ==

'''<big>Legacy implementation</big>'''

The formulation used for calculating water vapor saturation pressure in legacy code is based on the NASA GTE PEM-Tropics handbook and is as follows:

     REAL(fp), PARAMETER   :: A =   23.5518e+0_fp
     REAL(fp), PARAMETER   :: B = 2937.4e+0_fp
     REAL(fp), PARAMETER   :: C =   -4.9283e+0_fp
 
        ! Temperature at grid box (I,J,L)
        TEMP = State_Met%T(I,J,L)
 
        ! Saturation water vapor pressure in mbar 
        ! (from NASA GTE PEM-Tropics handbook)
        ESAT = ( 10e+0_fp**( A - ( B / TEMP ) ) ) * ( TEMP**C )

'''<big>New implementation</big>'''

As part of the updates, the saturation pressure of water vapor is now calculated using the empirical formulation in Nordquist, 1973. Note that it does not span the full dynamic range and incorporating a formulation that spans the entire possible dynamic range of the atmosphere would be a future science update.

Parameters in GeosCore/dao_mod.F:

      ! Empirical parameters for water vapor saturation pressure
      ! (Source: Nordquist, 1973. "Numerical Approximiations of
      !  Selected Meteorological Parameters Related to Cloud Physics"
      !  Text quality clarifications from Stipanuk, 1973. "Algorithms
      !  for Generating a Skew-T, Log P Diagram and Computing Selected
      !  Meteorological Quantities") 
      REAL(fp), PARAMETER   :: ESATP1  = 2.3832241e+1_fp   
      REAL(fp), PARAMETER   :: ESATP2  = -5.02808e+0_fp       
      REAL(fp), PARAMETER   :: ESATP3  = 8.1328e-3_fp    
      REAL(fp), PARAMETER   :: ESATP4  = 3.49149e+0_fp     
      REAL(fp), PARAMETER   :: ESATP5  = -1.3028844e+3_fp
      REAL(fp), PARAMETER   :: ESATP6  = -1.3816e-7_fp
      REAL(fp), PARAMETER   :: ESATP7  = 1.1344e+1_fp
      REAL(fp), PARAMETER   :: ESATP8  = -3.03998e-2_fp
      REAL(fp), PARAMETER   :: ESATP9  = -2.949076e+3_fp 

Parameters in GeosCore/dao_mod.F:

         !=============================================================
         ! Calculate water vapor saturation pressure [hPa]
         ! from temperature
         !=============================================================   
 
         EsatA = ESATP1 + ESATP2 * log10( State_Met%T(I,J,L) ) 
         EsatB = ESATP3 * 10**( ESATP4 + ESATP5 / State_Met%T(I,J,L) )
         EsatC = ESATP6 * 10**( ESATP7 + ESATP8 * State_Met%T(I,J,L) )
         EsatD = ESATP9 / State_Met%T(I,J,L) 
 
         ! Saturation water vapor pressure [mb]
         Esat = 10**( EsatA + EsatB + EsatC + EsatD )  

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Water Vapor Partial Pressure and Relative Humidity ==

'''<big>Legacy implementation</big>'''

Relative humidity is externally input with MERRA, GEOS-5, and GEOS-fp. If using GEOS-4 or GCAP, relative humidity is calculated in subroutine MAKE_RH in GeosCore/dao_mod.F legacy code. In that subroutine, the partial pressure of water vapor is calculated based on an assumption that specific humidity is defined as mass water vapor per mass dry air. The calculation uses a hard-coded approximation of the ratio of dry air molecular wt in [kg/mol] to water molecular weight in [g/mol]. 

Implementation in GeosCore/dao_mod.F:

        ! Pressure at midpoint of box (I,J,L)
        PRES = State_Met%PMID(I,J,L)
           
        ! Specific humidity in mb
        SHMB = State_Met%SPHU(I,J,L) * 1.6072e-3_fp * PRES
           
        ! Relative humidity as a percentage
        State_Met%RH(I,J,L) = ( SHMB / ESAT ) * 100e+0_fp

'''<big>New implementation</big>'''

In the new implementation, relative humidity is calculated in AIRQNT rather than MAKE_RH. The calculate is only necessary if using GEOS-4 or GCAP met data. The water vapor partial pressure calculation is different from legacy code, using a specific humidity definition of mass water vapor per total mass air. I derived the equation I use by substituting partial pressure * molecular weight for mass in the definition of specific humidity and solving for the partial pressure of water vapor.

NOTE: The updated relative humidity formulation, like the legacy formulation, does not span the entire dynamic range of the atmosphere. Using formulations for the entire dynamic range will be a future science update.

Parameters in Headers/CMN_GTCM_mod.F:

      ! AIRMW : Molecular weight of air [g/mol]
      REAL(fp), PARAMETER :: AIRMW    = 28.97e+0_fp 
 
      ! H2OMW : Molecular weight of water [g/mol]
      REAL(fp), PARAMETER :: H2OMW   = 18.016e+0_fp

Implementation in GeosCore/dao_mod.F:

         ! Convert specific humidity from [g/kg] to [kg/kg]
         SPHU_kgkg = State_Met%SPHU(I,J,L) * MCONV     
 
 #if   !defined( GEOS_5 ) && !defined( MERRA ) && !defined( GEOS_FP )
 
         !=============================================================
         ! Calculate water vapor partial pressure from specific humidity
         ! if using GEOS-4 or GCAP
         !=============================================================  
 
         ! Mol water vapor per mol moist air from specific humidity
         ! (this formulation is derived from the ideal gas law
         !  and the definition of specific humidity)
         AVGW_moist = AIRMW * SPHU_kgkg / ( AIRMW * SPHU_kgkg 
     &                      + H2OMW * ( 1.0e+0_fp - SPHU_kgkg ) ) 
 
         ! Water vapor partial pressure at grid box centroid [hPa]
         Ev_mid = State_Met%PMID(I,J,L) * AVGW_moist 
 
         !=============================================================
         ! Set grid box relative humidity [%] if using GEOS-4 or GCAP
         !============================================================= 
         !
         ! RH is the ratio of water vapor partial pressure to the
         ! saturation pressure of water, expressed as a percentage.
         !
         ! If RH is not externally input, first calculate the
         ! water vapor partial pressure from specific humidity.
         ! Then use it with the saturation pressure of water
         ! calculated above to set RH.
         !    
         !          Grid box           mol H2O
         !  Ev   =  moist air * ---------------------- 
         ! [hPa]    pressure     mol H2O + mol dry air 
         !          [hPa]       
         !                             
         !          H2O vapor partial pressure
         !  RH   =  -------------------------- * 100
         !           H2O saturation pressure
         !
 
         ! Relative humidity [%]
         State_Met%RH(I,J,L) = ( Ev_mid / Esat ) * 100e+0_fp
 
 #elif defined( GEOS_5 ) || defined ( MERRA ) || defined( GEOS_FP )
 
         !=============================================================
         ! Calculate water vapor partial pressures [hPa] from relative
         ! humidity if using GEOS-5, MERRA, or GEOS-FP. 
         !============================================================= 
 
         Ev_mid = State_Met%RH( I,J,L ) * RHCONV * Esat
 
 #endif


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Dry Air Partial Pressure  ==

'''<big>Legacy implementation</big>'''

Dry air pressures are not calculated in GEOS-Chem legacy code. Only moist air pressure is used.

'''<big>New implementation</big>'''

There are now several 3D dry air pressures to choose from:
#PEDGE_DRY - dry air partial pressure at bottom edge of grid box 
#PMID_DRY - dry air partial pressure at grid box mid-point, defined as arithmetic mean of edge pressures 
#PMEAN_DRY - dry air partial pressure at grid box altitude-weighted mid-point, defined as vertically-integrated mean pressure

The dry air partial pressures are calculated as the moist air pressure minus the water vapor partial pressure. For MERRA and GCAP, specific humidity is used to calculate water vapor partial pressure since those data sets do not include relative humidity. For GEOS-4, GEOS-5, and GEOS-FP, relative humidity is used. 

Implementation in GeosCore/dao_mod.F:

         !=============================================================
         ! Set grid box dry air partial pressures [hPa]
         ! at grid box centroid (arithmetic mean pressure level)
         !=============================================================  
 
         ! Partial pressure of dry air at box centroid [hPa]
         State_Met%PMID_DRY( I,J,L ) = State_Met%PMID(I,J,L) - Ev_mid
 
         !=============================================================
         ! Set grid box dry air partial pressures at grid box edges [hPa]
         ! Assume constant humidity across grid box.
         !=============================================================  
 
         Ev_edge = State_Met%PEDGE(I,J,L) * Ev_mid
     &             / State_Met%PMID(I,J,L)
  
         ! Partial pressure of dry air at lower edge of grid box [hPa]
         State_Met%PEDGE_DRY( I,J,L ) = State_Met%PEDGE(I,J,L) - Ev_edge
 
         ! Set dry air partial pressure for level LLPAR+1 lower edge
         IF ( L == LLPAR ) THEN     
            State_Met%PEDGE_DRY( I,J,L+1 ) = Pedge_Top - Ev_edge
         ENDIF
 
         !=============================================================
         ! Set grid box dry air partial pressures at grid box 
         ! altitude-weighted mean pressure [hPa]
         ! Assume constant humidity across grid box.
         !=============================================================  
 
         Ev_mean = State_Met%PMEAN(I,J,L) * Ev_mid
     &             / State_Met%PMID(I,J,L)
 
         ! Altitude-weighted mean dry air partial pressure of grid box [hPa]
         State_Met%PMEAN_DRY( I,J,L ) = State_Met%PMEAN(I,J,L) - Ev_mean


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

= Other Air Quantities =

== Air Mass ==

The air mass is calculated using Newton's second law (F=ma). The legacy method calculates the total mass in the grid box, including any water vapor present. 

In v11-01a, we introduced a new implementation that calculated the mass of dry air (no H2O included) by applying a dry air mass fraction scaling factor calculated from the ratio of dry and moist air densities. The updated v11-01a grid box mass was therefore less than the legacy mass by an order of 1% since the dry air partial density is less than moist air density which includes water vapor.

In v11-01h, we updated the previous dry air mass calculation to use the delta dry pressures derived from surface dry pressure and GMAO vertical grid parameters. This change enabled implementation of the [[GEOS-Chem_v11-01#Incorrect_incorporation_of_moisture_in_GEOS-Chem_air_quantities|moisture fix in v11-01g]], an update to correct the misuse of moisture of advection that caused a moisture signature in the v/v dry air concentration maps. 

'''<big>Legacy implementation</big>'''

Parameters in <tt>Headers/CMN_GCTM_mod.F</tt>:

      ! g0    : Acceleration due to gravity at earth's surface [m/s^2]
      REAL(fp), PARAMETER :: g0       =   9.8e+0_fp   
 
      ! g0_100 : 100 / g0
      REAL(fp), PARAMETER :: g0_100   = 100e+0_fp / g0

Implementation in GeosCore/dao_mod.F:

 #if defined( ESMF_ )
            AREA_M2                   = State_Met%AREA_M2(I,J,1)
 #else
            AREA_M2                   = GET_AREA_M2( I,J,L )
 #endif
 
         !===========================================================
         ! AD = (dry) mass of air in grid box (I,J,L) in kg, 
         ! given by:        
         !
         !  Mass    Pressure       100         1        Surface area 
         !        = difference   *  ---      *  ---   *   of grid box 
         !          in grid box         1           g           AREA_M2
         !
         !   kg         mb          Pa      s^2           m^2
         !  ----  =    ----      * ----  * -----  *      -----
         !    1          1            mb       m             1
         !===========================================================
         State_Met%AD(I,J,L) = State_Met%DELP(I,J,L) 
     &                       * G0_100 
     &                       * AREA_M2

'''<big>New implementation</big>'''

Parameters in <tt>Headers/CMN_GCTM_mod.F</tt>:

      ! g0    : Acceleration due to gravity at earth's surface [m/s^2]
      REAL(fp), PARAMETER :: g0       =   9.8e+0_fp     
 
      ! g0_100 : 100 / g0
      REAL(fp), PARAMETER :: g0_100   = 100e+0_fp / g0

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         ! Update dry pressure difference as calculated from the
         ! dry surface pressure
         State_Met%DELP_DRY(I,J,L) = GET_DELP_DRY(I,J,L)

         !===============================================================
         ! Set mass of dry air in grid box [kg]
         !==============================================================
         
         ! AD = mass of dry air in grid box (I,J,L) 
         !
         ! given by:        
         ! 
         !        Dry air pressure   Grid box   100 [Pa]   1 [s2]
         ! Mass = difference       * surface  * -------- * -----
         !        across grid        area       1 [hPa]    g [m]
         !        box [hPa]          [m2]                        
         !
         ! NOTES: 
         ! Dry air pressure difference across grid box is calculate
         ! from the surface dry pressure and GMAO A and B parameters
         ! (see GeosUtil/pressure_mod.F). The dry surface pressure
         ! that is used is calculated from the GMAO wet surface pressure 
         ! (or time-interpolated value) and the A and B parameters 
         ! (see SET_DRY_SURFACE_PRESSURE). It is not derived from the
         ! wet edge pressures. This distinction is important for
         ! mass conservation.
         State_Met%AD(I,J,L) = (State_Met%DELP_DRY(I,J,L) * G0_100)
     &                         * AREA_M2

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Box Height ==

The box height is calculated using the hydrostatic equation with the ideal gas law. The legacy implementation uses the gas constant for dry air in the idea gas law with the moist air pressures without taking into account the presence of water vapor. The new implementation incorporates the presence of water vapor by using the virtual temperature in place of ambient temperature to enable continued use of the moist air pressures with the dry air gas constant. The updated box height is greater than the legacy box height by an order of 1% which is consistent with a greater effective gas constant in the presence of water vapor.

'''<big>Legacy implementation</big>'''

Parameters in <tt>Headers/CMN_GCTM_mod.F</tt>:
 
      ! Rd    : Gas Constant in Dry Air [J/K/kg] 
      REAL(fp), PARAMETER :: Rd       = 287.0e+0_fp                 
 
      ! g0    : Acceleration due to gravity at earth's surface [m/s^2]
      REAL(fp), PARAMETER :: g0       =   9.8e+0_fp   
 
      ! Rdg0   = Rd    / g0
      REAL(fp), PARAMETER :: Rdg0     = Rd / g0

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         ! Pressure at bottom edge of grid box [hPa]
         P1                        = GET_PEDGE( I, J, L )
 
         ! Pressure at top edge of grid box [hPa]
         P2                        = GET_PEDGE( I, J, L+1 )
 
         !===========================================================
         ! BXHEIGHT is the height (Delta-Z) of grid box (I,J,L)
         ! in meters.
         !
         ! The formula for BXHEIGHT is just the hydrostatic eqn. 
         ! Rd = 287 J/K/kg is the value for the ideal gas constant
         ! R for air (M.W = 0.02897 kg/mol),  or
         ! Rd = 8.31 J/(mol*K) / 0.02897 kg/mol.
         !===========================================================
         State_Met%BXHEIGHT(I,J,L) = Rdg0
     &                             * State_Met%T(I,J,L)
     &                             * LOG( P1 / P2 )

'''<big>New implementation</big>'''

The new implementation takes water vapor into account by incorporating a "virtual" temperature in the ideal gas law, as done in Chapter 2 of Wallace and Hobbes' "Atmospheric Science: An Introductory Survey." The virtual temperature is now stored in State_Met for use elsewhere in GEOS-Chem.

Parameters in <tt>Headers/CMN_GCTM_mod.F</tt>:

      ! AIRMW : Molecular weight of air [g/mol]
      REAL(fp), PARAMETER :: AIRMW    = 28.97e+0_fp
 
      ! H2OMW : Molecular weight of water [g/mol]
      REAL(fp), PARAMETER :: H2OMW   = 18.016e+0_fp
 
      ! Rd    : Gas Constant in Dry Air [J/K/kg] 
      REAL(fp), PARAMETER :: Rd       = 287.0e+0_fp                 
 
      ! g0    : Acceleration due to gravity at earth's surface [m/s^2]
      REAL(fp), PARAMETER :: g0       =   9.8e+0_fp   
 
      ! Rdg0   = Rd    / g0
      REAL(fp), PARAMETER :: Rdg0     = Rd / g0

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         !=============================================================
         ! Set grid box height [m] 
         !=============================================================
 
         ! BXHEIGHT is the height (Delta-Z) of grid box (I,J,L) 
         ! in meters. It is calculated using the hypsometric equation.
         ! A virtual temperature is calculated to enable use of 
         ! of moist air pressures with dry air molecular weight.
         !  
         !              Gas constant   virtual
         !              for dry air  * grid box
         !                [J/K/mol]      temp [K]       bottom edge P
         ! height [m] = ------------------------- * ln(---------------)  
         !                     g [m/s2]                  top edge P
         !                           
         !  where,
         !
         !                                  Grid box temperature [K]
         !    Virtual  = --------------------------------------------------------
         !    Temp [K]            H20 partial pressure                MW_H2O
         !                   1 - --------------------       *  ( 1 - ------------ )
         !                        moist air  pressure                 MW_dryair
         !
         ! Source: Wallace and Hobbes "Atmospheric Science: An 
         !         Introductory Survey" 
         !
 
         ! Grid box virtual temperature [K]
         State_Met%TV( I,J,L ) = State_Met%T(I,J,L) / (1 -  
     &                           Ev_mid / State_Met%PMID(I,J,L) 
     &                           * ( 1 - H2OMW / AIRMW ) )
 
         ! Grid box box height [m]
         State_Met%BXHEIGHT( I,J,L ) = Rdg0 * State_Met%TV(I,J,L)
     &                                 * LOG( State_Met%PEDGE(I,J,L) 
     &                                 / PEdge_Top )


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Air Volume ==

The air quantity updates result in a higher grid box volume due to changes to the box height calculation (see above).

'''<big>Legacy implementation</big>'''

Implementation in GeosCore/dao_mod.F:

         !===========================================================
         ! AIRVOL is the volume of grid box (I,J,L) in meters^3
         !
         ! AREA_M2 is the Delta-X * Delta-Y surface area of grid
         ! boxes (I,J,L=1), that is, at the earth's surface.
         !
         ! Since the thickness of the atmosphere is much smaller 
         ! than the radius of the earth, we can make the "thin 
         ! atmosphere" approximation, namely:
         !
         !               (Rearth + h) ~ Rearth
         !
         ! Therefore, the Delta-X * Delta-Y surface area of grid
         ! boxes that are above the earth's surface will be 
         ! approx. the same as AREA_M2.  Thus we are justified 
         ! in using AREA_M2 for grid boxes (I, J, L > 1)
         !===========================================================
         State_Met%AIRVOL(I,J,L) = State_Met%BXHEIGHT(I,J,L) * AREA_M2

'''<big>New implementation</big>'''

Same as legacy.

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Air Density ==

'''<big>Legacy implementation</big>'''

Legacy GEOS-Chem calculates air density as the total grid box moist air mass divided by volume (the product of grid box area and height). 

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         !===========================================================
         ! AIRDEN = density of air (AD / AIRVOL) in kg / m^3 
         !===========================================================
         State_Met%AIRDEN(I,J,L) = State_Met%AD(I,J,L) 
     &                           / State_Met%AIRVOL(I,J,L)

'''<big>New implementation</big>'''

In the new implementation there are two air densities stored in State_Met:

* <tt>MAIRDEN</tt> - wet air density (includes water vapor mass), in [kg/m3]
* <tt>AIRDEN</tt> - dry air density (does not include water vapor mass), in [kg/m3]
* <tt>AIRNUMDEN</tt> - dry air density (does not include water vapor mass), in [molec/cm3]

The moist air density is calculated using the ideal gas law with weighted partial pressures. The dry air density also used the ideal gas law in its formulation starting in v11-01a, but we changed this in v11-01g to use legacy mass divided by volume in order to conserve mass in the moisture fix. Updating the air density to use the ideal gas law in a way that conserves mass following the moisture fix will be addressed in v11-02.

The moist density values (MAIRDEN in v11-01, AIRDEN in v10-01) are consistent with the legacy air density definition of total mass per volume of grid box. Maintaining this definition while switching to use of the ideal gas law required using the altitude-weighted mean pressure (new variable PMEAN). 

Parameters in <tt>Headers/CMN_GTCM_mod.F</tt>:

      ! AIRMW : Molecular weight of air [g/mol]
      REAL(fp), PARAMETER :: AIRMW    = 28.97e+0_fp
 
      ! H2OMW : Molecular weight of water [g/mol]
      REAL(fp), PARAMETER :: H2OMW   = 18.016e+0_fp
 
      ! RSTARG : Universal gas constant [J/K/mol] 
      REAL(fp), PARAMETER :: RSTARG   = 8.31450e+0_fp

Parameters in <tt>GeosCore/dao_mod.F</tt>:

      ! Conversion factors
      REAL(fp), PARAMETER   :: PCONV   = 1.00e+2_fp         ! hPa -> Pa 
      REAL(fp), PARAMETER   :: MCONV   = 1.00e-3_fp          ! g -> kg

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         !===========================================================
         ! Set grid box wet air density [kg/m3]
         !==========================================================
 
         !  MAIRDEN = density of moist air [kg/m^3],
         ! given by:
         !          
         !            Partial          Molec     Partial              Molec
         !            pressure   * wt of   +  pressure of   * wt of
         !            of dry air      air           water vapor    water    
         ! Moist      [hPa]         [g/mol]   [hPa]              [g/mol] 
         ! Air     =  ------------------------------------------------       
         ! Density    Universal gas constant * Temp * 1000  * 0.01
         !                   [J/K/mol]                     [K]      [g/kg]   [hPa/Pa]
         !
         State_Met%MAIRDEN(I,J,L) = ( State_Met%PMEAN_DRY(I,J,L) * AIRMW 
     &                                + Ev_mean * H2OMW ) * PCONV * MCONV
     &                                / RSTARG / State_Met%T(I,J,L)
 
         ! Set grid box dry air density [kg/m3]
         !
         ! NOTE: Air density is derived from dry surface pressure
         ! following implementation of the moisture fix, and therefore
         ! its equation must be dry air mass divided by volume. 
         ! This is because the level pressures derived from the dry 
         ! surface pressure may be used for extracting mass per level, 
         ! but not a representative pressure for that level. The ideal 
         ! gas law requires a representative level pressure. Eventually 
         ! a representative pressure must be derived that, when used in
         ! the ideal gas law, will replicate the air density defined as
         ! mass/volume. Moving to use of the ideal gas law is necessary 
         ! for grid-independence (ewl, 9/16/16)
         State_Met%AIRDEN(I,J,L) = State_Met%AD(I,J,L)  
     &                             / State_Met%AIRVOL(I,J,L)

         ! Set grid box dry air number density [molec/cm3]
         State_Met%AIRNUMDEN(I,J,L) = State_Met%AIRDEN(I,J,L)
     &                                * 1e-3_fp *  AVO / AIRMW


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Water Vapor Mixing Ratio ==

AVGW is the water vapor volume mixing ratio in units of moles water vapor per mole dry air. This quantity existed prior to the addition of the H2O tracer. It is currently used in aerosols and in RRTMG.

'''<big>Legacy implementation</big>'''

The water vapor mixing ratio is set in subroutine MAKE_AVGW in GeosCore/dao_mod.F. It is calculated by multiplying specific humidity [g H2O / kg air] by the ratio of dry air to water molecular weights.  This approach assumes that specific humidity represents water vapor mass per dry air mass which is incorrect. In addition, the legacy implementation hard-codes the molecular weight of water using a low precision value.

Implementation in <tt>GeosCore/dao_mod.F</tt>:

      ! Conversion factor
      REAL(fp), PARAMETER   :: HCONV = 28.97e-3_fp / 18.0e+0_fp 
 
      State_Met%AVGW(I,J,L) = HCONV * State_Met%SPHU(I,J,L)

'''<big>New implementation</big>'''

AVGW is now computed in AIRQNT rather than MAKE_AVGW. The calculation uses existing global parameters for the molecular weights of air and water. It uses SPHU/(1-SPHU) rather than SPHU, effectively replacing [kg water vapor / kg total air] with [kg water vapor / kg dry air ].

Parameters in <tt>GeosCore/dao_mod.F</tt>:

      REAL(fp), PARAMETER   :: MCONV   = 1.00e-3_fp      ! g -> kg

Implementation in <tt>GeosCore/dao_mod.F</tt>:

         !=============================================================
         ! Set AVGW, the water vapor volume mixing ratio from 
         ! specific humidity [mol H2O / mol dry air]
         !=============================================================  
         !
         ! vol H2O       dry air    H2O         kg H2O      kg wet air
         ! ----------  = molec wt / molec wt * ---------- * ---------
         ! vol dry air   [g/mol]    [g/mol]    kg wet air   kg dry air
         !
         !      thus AVGW = AIRMW * SPHU_kgkg / ( H2OMW * (1-SPHU_kgkg) )
         !
         ! where, 
         !        SPHU_kgkg = specific humidity in [kg/kg]
         !
 
         ! Convert specific humidity from [g/kg] to [kg/kg]
         SPHU_kgkg = State_Met%SPHU(I,J,L) * MCONV     
         
         ! Water vapor volume mixing ratio [v/v dry air]
         State_Met%AVGW(I,J,L) = AIRMW * SPHU_kgkg 
     &                           / ( H2OMW * (1.0e+0_fp - SPHU_kgkg ) )


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

= Changes Beyond AirQnt =

== SET_H2O_TRAC in UCX_Mod ==

There are three distinct changes in routine SET_H2O_TRAC in UCX_mod. 

#The routine calculates the H2O tracer mass from specific humidity by multiplying by mass. Since State_Met%AD is now dry air mass, and tracer units are now dry mixing ratio, the updates include replacement of mass with specific humidity (a mass ratio) in the calculation. 
#The routine calculates relative humidity from SPHU if some conditions are met and uses the same formulation as the RH calculation in dao_mod. Since AIRQNT now uses an updated formulation for water vapor saturation pressure, the formulation SET_H2O_TRAC is now also updated to match it.
#Since moisture is now incorporated into several State_Met variables beyond SPHU and RH, AIRQNT is now called at the end of SET_H2O_TRAC to updated all air quantities if SPHU and RH were reset within this routine.

'''<big>Legacy implementation</big>'''
 
         IF ( READ_SPHU ) THEN
            STT(I,J,L,IDTH2O) = State_Met%SPHU(I,J,L) * 1.e-3_fp *
     &                          State_Met%AD(I,J,L)
         ELSE
            ! Set specific humidity to transported H2O
            State_Met%SPHU(I,J,L) = STT(I,J,L,IDTH2O) * 1.e+3_fp /
     &                              State_Met%AD(I,J,L)
 
             ...  (legacy RH calculation here) ...
 
            State_Met%RH(I,J,L) = ( SHMB / ESAT ) * 100e+0_fp
         ENDIF

'''<big>New implementation</big>'''

         IF ( READ_SPHU ) THEN
 
            ! Calculate specific humidity [g H2O/kg total air] as [kg/kg] 
            SPHU_kgkg = State_Met%SPHU(I,J,L) * 1.e-3_fp
             
            ! Set H2O species concentration [kg/kg dry] using SPHU [kg/kg]
            ! without using box mass (ewl, 7/18/16)
            Spc(I,J,L,id_H2O) = SPHU_kgkg / ( 1.0e+0_fp - SPHU_kgkg )
 
         ELSE
 
            ! Calculate specific humidity in [kg H2O / kg total air] 
            ! using transported H2O [kg/kg dry] (ewl, 7/18/16)
            SPHU_kgkg =  Spc(I,J,L,id_H2O) / 
     &                  ( 1.0e+0_fp + Spc(I,J,L,id_H2O) ) 
 
 
            ! Set State_Met specific humidity [g/kg]
            State_Met%SPHU(I,J,L) = SPHU_kgkg * 1.e+3_fp
   
            ... (new RH calculate here) ...
 
            State_Met%RH(I,J,L) = ( Ev_mid / Esat ) * 100e+0_fp 
 
         ENDIF
 
      ... 
 
      ! If humidity was updated, update all moist-dependent air quantities
      ! and species mixing ratio with the new moisture content (ewl, 4/29/15) 
      IF ( LActiveH2O .and. ( .not. SetStrat ) ) THEN
         CALL AIRQNT( am_I_Root, Input_opt, State_Met, State_Chm, RC, 
     &                update_mixing_ratio=.TRUE. )
      ENDIF

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 17:40, 29 September 2016 (UTC)

== MOIST_QQ in Wetscav_Mod ==

MOISTQ is the tendency in specific humidity (kg/kg/s) that comes with the GEOS-4, GEOS-5, and GCAP met fields. In wetscav_mod routine MAKE_QQ, MOISTQ is multipled by State_Met%AIRDEN to convert units. Since AIRDEN is now dry air density and MOISTQ is per moist air mass, we must change AIRDEN to MAIRDEN (moist air density) where it was used with MOISTQ. 

'''<big>Legacy implementation example</big>'''

            QQ(L,I,J) = FRAC * MOISTQ(I,J,L) * AIRDEN(I,J,L) / 1000.0

'''<big>New implementation example</big>'''

            QQ(L,I,J)    = FRAC                      *
     &                     State_Met%MOISTQ(I,J,L)   *
     &                     State_Met%MAIRDEN(I,J,L)   / 1e+3_fp

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Convection ==

In the legacy implementation, tracer concentration units in convection_mod are mixing ratio [v/v dry air] and TCVV ( dry air molecular wt / tracer molecular wt ) is used to convert units for diagnostic and mercury calculations. 

In the new Implementation:

#Tracer units are now [kg/kg dry air]
#Use delta dry pressure instead of total air mass (including H2O) for ND38 diagnostic (convective loss) unit conversions within DO_CONVECTION. Note that the delta dry pressure is derived from the dry surface pressure implemented as part of the [[GEOS-Chem_v11-01#v11-01h|v11-01h]] moisture fix.

Changes made in DO_CONVECTION are as follows:

'''<big>Legacy implementation</big>'''

 ! DO_CONVECT begins here!
                ....    
                TMP3D = TMP3D * State_Met%AD / Input_Opt%TCVV(IT) 
                ... 
     END SUBROUTINE DO_CONVECTION

'''<big>New implementation</big>'''

 ! DO_CONVECT begins here!
                ...
                TMP3D = TMP3D * State_Met%DELP_DRY(I,J,L) * 
     &                 G0_100 * AREA_M2
                ...
      END SUBROUTINE DO_CONVECTION

--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)

== Transport ==

In the legacy implementation, tracer concentration units in transport_mod and tpcore routines are mixing ratio [v/v dry air]. During advection, however, delta pressure is used as an air mass proxy even though it is represents total air rather than dry air mass. Without taking into account the mass to mole conversion for wet air, error is therefore introduced in the advection mass balance. In addition, tracer concentrations are adjusted within a residual mass correction following advection in order to preserve mass conservation. The correction uses the function GET_AIR_MASS to retrieve the before and after air mass for each grid box. However, that function returns the moist air mass rather than dry air mass which is inconsistent with the v/v dry air units. 

In v11-01a, we made the following updates to correct these discrepancies and simplify the code:

#Convert tracer units from v/v dry air to kg/kg total air just before advection and back to v/v dry air after advection. 
#Remove the residual mass correction since it is no longer necessary if tracer units are in kg/kg total air during advection. Removal of the residual mass correction also allows removal of function GET_AIR_MASS.
#Call AIRQNT to update all air quantities following the pressure update. Note that previously the post-advection AIRQNT call occurred in main.F. Moving it to within transport is not necessary but simplifies main.F.

Validation of the v11-01a moisture corrections revealed a new problem in advection. Inert tracer dry mixing ratio maps contained patterns resembling the moisture signature in the atmosphere. This was due to lack of moisture handling within advection, specifically that tracer mass was distributed relative to moist rather than dry air pressure and therefore preferentially pooled within grid boxes that contained more water vapor mass. While the dry vs. wet air mass discrepancy was resolved in v11-01a, the chosen unit, total mixing ratio, was at odds with the advection scheme which assumes zero moisture.

Meemong Lee and Richard Weidner (JPL) issued a [http://wiki.seas.harvard.edu/geos-chem/images/6/6f/GCA-M2O2-Final.pdf JPL Publication] recommending a fix for this problem in the Adjoint model. Their Adjoint fix consisted of (1) deriving dry surface pressures from GMAO moist surface pressures and specific humidity, and (2) replacing moist pressures with dry pressures throughout the model.

We adapted JPL’s Adjoint fix to the GEOS-Chem forward model in GEOS-Chem v11-01h. The GEOS-Chem implementation (1) corrects the transport mass distribution problem by applying JPL’s fix to transport and unit conversions, (2) preserves the moisture updates incorporated into v11-01a, and (3) preserves inert tracer mass conservation. 

For advection, the v11-01h update required the following: 

#Removal of the dry mixing ratio to total mixing ratio tracer unit conversions in advection routines such that dry mixing ratio is the concentration unit used.
#Set the surface pressure at the midpoint of the dynamic timestep (P_TP1 in advection) to the previous ("floating") dry surface pressure derived from the time-interpolated GMAO surface pressure with moisture in the column removed.
#Set the surface pressure at the end of the dynamic timestep (P_TP2 in advection) to the new dry surface pressure derived in the same way as the "floating" dry surface pressure.
#Pass the new pressures rather than the legacy total pressures (including moisture) to advection
#Update the dry "floating" pressure with the dry surface pressure post-advection.
#Update the dry mixing ratio post-advection with the dry air mass per level calculated using the new "floating" dry pressure.

It is very important to note that the '''new dry surface pressure should NOT be used for constructing a vertical profile of dry air pressure'''. The dry air pressure derived from dry surface pressure and GMAO vertical grid parameters A and B would not reflect an accurate height profile. For this reason, 3D dry pressures in State_Met%PEDGE_DRY and State_Met%PMID_DRY continue to be constructed using GMAO wet surface pressures, A and B parameters, and specific humidity.

It is also important to note that the dry air mass per level calculated using the dry surface pressure is close but not exactly the same as the dry air mass calculated using the legacy method (with State_Met variables PEDGE and SPHU). Both methods produce valid approximations of level mass, but '''using dry air mass derived from the new dry surface pressure is required for species unit conversions to ensure mass conservation'''. This is because species are advected in units of mass mixing ratio, with dry air mass in the denominator. Retaining the same global species mass before and after advection, therefore, requires unit conversions that use a dry air mass definition consistent with that used in advection. 
     
--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 18:43, 29 September 2016 (UTC)