Air Quantity Updates for v11-01: Difference between revisions

From Geos-chem
Jump to navigation Jump to search
 
(7 intermediate revisions by one other user not shown)
Line 1: Line 1:
= 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 =
= 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 introduced in v11-01a remedy this problem and fix some additional moisture-related bugs discovered during implemention. In addition, where tracer units were previously [v/v], they are now consistently either [kg/kg total air] (advection and convection) or [v/v dry air] (everywhere else).  
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 species units, discussed separately at the [[GEOS-Chem_species_units|GEOS-Chem species units wiki page]].  


Below is a high-level summary of the updates. All of these updates were incorporated in GEOS-Chem v11-01a:
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:'''
*'''Dry air partial pressure replaces wet atmospheric pressure for:'''
Line 9: Line 19:
#gas partial pressure calculations
#gas partial pressure calculations


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


*'''New State_Met variables:'''
*'''New State_Met variables:'''
#PMID_DRY - dry air partial pressure at pressure mid-point [hPa]
#PMID_DRY - dry air partial pressure at pressure mid-point [hPa]
#PEDGE_DRY - dry air partial pressure at level lower edge [hPa]
#PEDGE_DRY - dry air partial pressure at level lower edge [hPa]
#PMEAN - altitude-weighted mean pressure [hPa]
#PMEAN_DRY- dry air partial pressure at altitude-weight mean pressure [hPa]
#MOISTMW - molecular wt of moist air in grid box [g/mol]
#ADMOIST - moist air mass (previously AD) [kg/m3]
#MAIRDEN - moist air density (previously AIRDEN) [kg/m3]
#MAIRDEN - moist air density (previously AIRDEN) [kg/m3]
#AIRNUMDEN - dry air density [molec/cm3]
#TV - virtual temperature [K]
#TV - virtual temperature [K]


*'''Pre-existing State_Met variables that have modified definitions and algorithms:'''
*'''Pre-existing State_Met variables that have modified definitions and algorithms:'''
#AD - now defined as dry air mass (previously total air mass)
#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)
#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:'''
*'''Pre-existing State_Met variables that have unchanged definitions but modified algorithms:'''
Line 39: Line 46:
#SPHU (specific humidity [g/kg]) is used with air quantities for unit conversion
#SPHU (specific humidity [g/kg]) is used with air quantities for unit conversion


NOTE: All of the above changes are explained in detail within this wiki page. 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 some circumstances.
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 and Humidity =
Line 87: Line 94:
#PEDGE - wet air pressure at bottom edge of grid box  
#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
#PMID - wet air pressure at grid box pressure mid-point, defined as arithmetic mean of edge pressures
#PMEAN - wet air pressure at grid box altitude-weighted mid-point, defined as vertically-integrated mean pressure
#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 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.
Line 132: Line 139:
         ! average P(z) over z. Use in the ideal gas law yields a
         ! average P(z) over z. Use in the ideal gas law yields a
         ! mean density equivalent to total mass per volume in grid box.
         ! mean density equivalent to total mass per volume in grid box.
         ! NOTE: Use of PMEAN is a temporary fix to preserve current
         PMEAN = State_Met%DELP(I,J,L) / log( State_Met%PEDGE(I,J,L) /
        ! usage of AIRDEN as total mass per volume.
        State_Met%PMEAN( I,J,L ) = State_Met%DELP(I,J,L) /
    &                              log( State_Met%PEDGE(I,J,L) /
     &                                  PEdge_Top )  
     &                                  PEdge_Top )  


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


== Water Vapor Saturation Pressure ==
== Water Vapor Saturation Pressure ==
Line 345: Line 349:
= Other Air Quantities =
= Other Air Quantities =


== Moist Air Molecular Weight ==
== Air Mass ==
 
'''<big>Legacy implementation</big>'''
 
Moist air molecular weight is not calculated in GEOS-Chem legacy code. Dry air molecular weight is used with moist air quantities.
 
'''<big>New implementation</big>'''
 
Moist air molecular weight is calculated by weighting the dry air and H2O molecular weights by the respective partial pressures.
 
        !=============================================================
        ! Set grid box moist molecular weight [g/mol]
        ! This quantity is useful if using TCVV with moist mixing ratio.
        ! Since we assume constant moisture and T across the grid box,
        ! MOISTMW is constant within the box.
        !============================================================= 
        ! Grid box moist air molecular weight [g/mol]
        State_Met%MOISTMW(I,J,L) = ( State_Met%PMID_DRY(I,J,L)
    &                              * AIRMW + Ev_mid * H2OMW )
    &                              / State_Met%PMID(I,J,L)


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)
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.


== Air Mass ==
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.


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. The new implementation calculates the mass of dry air only by applying a dry air mass fraction scaling factor calculated from the ratio of dry and moist air densities. The updated grid box mass is 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.  
 
NOTE: Grid box air mass is sometimes calculated using air density and volume rather than delta pressure and area. As we change units from grid box mass to grid box mass mixing ratio, these instances will be changed to calculate mass/area using delta pressure rather than air density and air volume or box height. This will enable isolating air density use to local condition representation only (not used to extrapolate air mass for the entire grid box) and thus use of PMID in the ideal gas law calculation of air density will be universally appropriate. For now, PMEAN is used to calculate air density to preserve the total mass/volume definition.  


'''<big>Legacy implementation</big>'''
'''<big>Legacy implementation</big>'''
Line 418: Line 400:
       ! g0_100 : 100 / g0
       ! g0_100 : 100 / g0
       REAL(fp), PARAMETER :: g0_100  = 100e+0_fp / g0
       REAL(fp), PARAMETER :: g0_100  = 100e+0_fp / g0
      ! RSTARG : Universal gas constant [J/K/mol]
      REAL(fp), PARAMETER :: RSTARG  = 8.31450e+0_fp


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


    ! Convert specific humidity from [g/kg] to [kg/kg]
        ! Update dry pressure difference as calculated from the
    SPHU_kgkg = State_Met%SPHU(I,J,L) * MCONV   
        ! 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]
         ! Set mass of dry air in grid box [kg]
         !==============================================================
         !==============================================================
       
        ! AD = mass of dry air in grid box (I,J,L)
         !
         !
         ! AD = mass of dry air in grid box (I,J,L)
         ! given by:       
         !   = mass of wet air * mass ratio of dry air to wet air
        !
        !        Dry air pressure  Grid box   100 [Pa]  1 [s2]
         ! Mass = difference      * surface  * -------- * -----
        !        across grid        area      1 [hPa]    g [m]
        !        box [hPa]          [m2]                       
         !
         !
         ! given by:      
         ! NOTES:  
         !
         ! Dry air pressure difference across grid box is calculate
        !        Air pressure   Grid box  Mass      100 [Pa]  1 [s2]
         ! from the surface dry pressure and GMAO A and B parameters
        ! Mass = difference   * surface  * ratio of * -------- * -----
         ! (see GeosUtil/pressure_mod.F). The dry surface pressure
        !        across grid   area      dry to    1 [hPa]    g [m]
         ! that is used is calculated from the GMAO wet surface pressure
        !        box [hPa]      [m2]      wet air               
         ! (or time-interpolated value) and the A and B parameters
         !
         ! (see SET_DRY_SURFACE_PRESSURE). It is not derived from the
        ! where,
         ! wet edge pressures. This distinction is important for
        !
         ! mass conservation.
        !  Ratio of dry     mass density of dry air
         State_Met%AD(I,J,L) = (State_Met%DELP_DRY(I,J,L) * G0_100)
         ! to wet air  =  --------------------------
     &                        * AREA_M2
        !  mass [.]        mass density of moist air
         !
        ! NOTES:
        !    (1) Air pressure difference across grid box is the same for
         !       wet (DELP) and dry air partial pressures due to
        !        the assumption of constant humidity within the grid box.
         !   (2) Specific humidity in GEOS-Chem (SPHU) is stored in [g/kg]
         !
         !==============================================================
         State_Met%ADMOIST(I,J,L) = State_Met%DELP(I,J,L) * G0_100
    &                              * AREA_M2
        State_Met%AD(I,J,L) = State_Met%ADMOIST(I,J,L)
     &                        * ( 1.0e+0_fp - SPHU_kgkg )
 


--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)
--[[User:Lizzie Lundgren|Lizzie Lundgren]] ([[User talk:Lizzie Lundgren|talk]]) 16:32, 18 May 2015 (UTC)
Line 615: Line 586:
In the new implementation there are two air densities stored in State_Met:
In the new implementation there are two air densities stored in State_Met:


* <tt>MAIRDEN</tt> - wet air density (includes water vapor mass)
* <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)
* <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]


Both air densities are calculated using the ideal gas law with weighted partial pressures.
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 updates include calculation of (1) the moist air density (MAIRDEN) from the ideal gas law with weighted dry air and water vapor partial pressures, and (2) the dry air partial density (AIRDEN)n from the ideal gas law using the dry air partial pressure only. The mass density of dry air in the grid box is less than the moist air density by up to around 2% at the surface.
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).  
 
The updated density values 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). In the future, we will adjust all instancescalculation using the updated dry air mass (AD) divided by the updated volume (AIRVOL).


Parameters in <tt>Headers/CMN_GTCM_mod.F</tt>:
Parameters in <tt>Headers/CMN_GTCM_mod.F</tt>:
Line 662: Line 632:
     &                                / RSTARG / State_Met%T(I,J,L)
     &                                / RSTARG / State_Met%T(I,J,L)
   
   
        !===========================================================
         ! Set grid box dry air density [kg/m3]
         ! Set grid box dry air density [kg/m3]  
         !
         !==========================================================
        ! NOTE: Air density is derived from dry surface pressure
         State_Met%AIRDEN(I,J,L) = State_Met%PMEAN_DRY(I,J,L) * PCONV
        ! following implementation of the moisture fix, and therefore
     &                            * AIRMW * MCONV
        ! its equation must be dry air mass divided by volume.
    &                            / RSTARG / State_Met%T(I,J,L)
        ! 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




Latest revision as of 17:35, 6 December 2016

Problem Description

GCST member 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-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 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-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 species units, discussed separately at the 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:
  1. air number density unit conversions
  2. gas partial pressure calculations
  • Units of tracer concentration (State_Chm%Tracers), previously in [v/v], are now:
  1. [kg/kg dry air] in convection and advection
  2. [v/v dry air] everywhere else
  • New State_Met variables:
  1. PMID_DRY - dry air partial pressure at pressure mid-point [hPa]
  2. PEDGE_DRY - dry air partial pressure at level lower edge [hPa]
  3. MAIRDEN - moist air density (previously AIRDEN) [kg/m3]
  4. AIRNUMDEN - dry air density [molec/cm3]
  5. TV - virtual temperature [K]
  • Pre-existing State_Met variables that have modified definitions and algorithms:
  1. AD - now defined as dry air mass (previously total air mass), and derived from surface dry pressure
  2. 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:
  1. RH - now uses Nordquist 1973 formulation for water vapor saturation pressure
  2. BXHEIGHT - now uses virtual temperature in the hypsometric equation
  3. 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:
  1. AIRVOL - impacted by change in box height calculation
  • Moist air mass and density replace dry air mass and density when:
  1. MOISTQ (tendency in specific humidity [kg/kg/s]) is used with air quantities for unit conversion
  2. 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:

  1. Air number density conversions (ucx_mod, mercury_mod, gasconc)
  2. Partial pressure calculations (ucx_mod, sulfate_mod)
  3. Surface pressure diagnostic ND31 (per request of Alex Turner for CH4
  4. 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:

  1. Pressure level diagnostic outputs
  2. Cloud top height calculation
  3. Potential temperature calculation
  4. Volume cloud fraction calculation
  5. Pressure level thresholding (i.e. height proxy)
  6. Barometric formula
  7. Hypsometric equation (now uses virtual temperature in stat_chem_mod and vdiff_mod)
  8. Mean free path calculations

--Lizzie Lundgren (talk) 15:54, 18 May 2015 (UTC)

Moist Air Pressure

Legacy implementation

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

New implementation

There are now several moist air pressures to choose from:

  1. PEDGE - wet air pressure at bottom edge of grid box
  2. PMID - wet air pressure at grid box pressure mid-point, defined as arithmetic mean of edge pressures
  3. 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 -> deg C    
REAL(fp), PARAMETER   :: PCONV   = 1.00e+2_fp      ! hPa -> Pa      
REAL(fp), PARAMETER   :: RHCONV  = 1.00e-2_fp      ! % -> fraction     
REAL(fp), PARAMETER   :: MCONV   = 1.00e-3_fp      ! g -> 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 )

--Lizzie Lundgren (talk) 16:43, 29 September 2016 (UTC)

Water Vapor Saturation Pressure

Legacy implementation

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 )

New implementation

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 )

--Lizzie Lundgren (talk) 16:32, 18 May 2015 (UTC)

Water Vapor Partial Pressure and Relative Humidity

Legacy implementation

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

New implementation

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


--Lizzie Lundgren (talk) 16:32, 18 May 2015 (UTC)

Dry Air Partial Pressure

Legacy implementation

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

New implementation

There are now several dry air pressures to choose from:

  1. PEDGE_DRY - dry air partial pressure at bottom edge of grid box
  2. PMID_DRY - dry air partial pressure at grid box mid-point, defined as arithmetic mean of edge pressures
  3. 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


--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 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.

Legacy implementation

Parameters in Headers/CMN_GCTM_mod.F: 

     ! 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

New implementation

Parameters in Headers/CMN_GCTM_mod.F: 

     ! 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: 

        ! 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

--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.

Legacy implementation

Parameters in Headers/CMN_GCTM_mod.F: 

     ! 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 GeosCore/dao_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 )

        !===========================================================
        ! 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 )

New implementation

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 Headers/CMN_GCTM_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

     ! 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 GeosCore/dao_mod.F: 

        !=============================================================
        ! 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 )


--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).

Legacy implementation

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

New implementation

Same as legacy.

--Lizzie Lundgren (talk) 16:32, 18 May 2015 (UTC)

Air Density

Legacy implementation

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 GeosCore/dao_mod.F: 

        !===========================================================
        ! 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)

New implementation

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

  • MAIRDEN - wet air density (includes water vapor mass), in [kg/m3]
  • AIRDEN - dry air density (does not include water vapor mass), in [kg/m3]
  • AIRNUMDEN - 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 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

     ! RSTARG : Universal gas constant [J/K/mol] 
     REAL(fp), PARAMETER :: RSTARG   = 8.31450e+0_fp

Parameters in GeosCore/dao_mod.F: 

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

Implementation in GeosCore/dao_mod.F: 

        !===========================================================
        ! 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


--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.

Legacy implementation

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 GeosCore/dao_mod.F: 

     ! 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)

New implementation

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 GeosCore/dao_mod.F: 

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

Implementation in GeosCore/dao_mod.F: 

        !=============================================================
        ! 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 ) )


--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.

  1. The routine calculates the H2O tracer mass from specific humidity by multiplying by mass. Since State_Met%AD is now dry air mass, the updates include changing State_Met%AD to State_Met%ADMOIST in this routine.
  2. 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.
  3. 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 for at least one grid box. A flag is used to tracer whether SPHU and RH are reset.

Legacy implementation 

        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

New implementation 

        IF ( READ_SPHU ) THEN
           STT(I,J,L,IDTH2O) = State_Met%SPHU(I,J,L) * 1.e-3_fp *
    &                          State_Met%ADMOIST(I,J,L)
        ELSE
           ! Calculate specific humidity [g/kg] from transported H2O
           SPHU = STT(I,J,L,IDTH2O) * 1.e+3_fp 
    &               / State_Met%ADMOIST(I,J,L)

           State_Met%SPHU(I,J,L) = SPHU

           ... (new RH calculate here) ...

           State_Met%RH(I,J,L) = ( Ev_mid / Esat ) * 100e+0_fp 

           IF ( .not. RESET_MOISTURE ) THEN
              RESET_MOISTURE = .TRUE.
           ENDIF
        ENDIF

     IF ( RESET_MOISTURE ) THEN
        CALL AIRQNT( am_I_Root, State_Met, RC )
     ENDIF

--Lizzie Lundgren (talk) 16:32, 18 May 2015 (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.

Legacy implementation example 

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

New implementation example 

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

--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. While units are [v/v dry air], comments in the code indicate moist mixing ratio is the assumed tracer concentration unit for the convection scheme.

In the new Implementation:

  1. Convert tracer units from [v/v dry air] to [kg/kg total air] (wet) at start of DO_CONVECTION and back to [v/v dry] at end of DO_CONVECTION. The conversions are routines Convert_DryVV_to_DryMMR, Convert_DryMMR_to_MoistMMR, Convert_MoistMMR_to_DryMMR, and Convert_DryMMR_to_DryVV defined in GeosUtil/unitconv_mod.
  2. Remove TCVV from diagnostics and mercury calculations since there it is no longer necessary to convert between moles and mass.
  3. Use ADMOIST instead of AD for diagnostic unit conversion within DO_CONVECTION.

Changes made in DO_CONVECTION are as follows:

Legacy implementation 

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

New implementation 

! DO_CONVECT begins here!
     ! Convert State_Chm%Tracers from [v/v dry air] to 
     ! [kg/kg total air] prior to convection (ewl, 6/23/15)
     CALL Convert_DryVV_to_DryMMR( am_I_Root, N_TRACERS, Input_Opt,
    &                        State_Chm, RC)
     CALL Convert_DryMMR_to_MoistMMR( am_I_Root, N_TRACERS, State_Met, 
    &                                 State_Chm, RC )
      ...
               TMP3D = TMP3D * State_Met%ADMOIST
      ....
     ! Convert tracer concentration from kg/kg total air to 
     ! mol/mol dry air after convection (ewl, 6/23/15)
     CALL Convert_MoistMMR_to_DryMMR( am_I_Root, N_TRACERS, State_Met, 
    &                               State_Chm, RC ) 
     CALL Convert_DryMMR_to_DryVV( am_I_Root, N_TRACERS, Input_Opt,
    &                        State_Chm, RC)
     END SUBROUTINE DO_CONVECTION

--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] and TCVV values ( dry air molecular wt / tracer molecular wt ) are used to calculate mass flux diagnostics. During advection, 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 <-> mole conversion for wet air, error is introduced in the advection mass balance.

To preserve mass conservation, tracer concentrations are adjusted within a residual mass correction following advection. The correction uses the function GET_AIR_MASS to retrieve the before and after air mass for each grid box. However, it returns the moist air mass rather than dry air mass.

In the new implementation:

  1. 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.
  2. 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.
  3. 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.
  4. TCVV is removed from all diagnostic unit conversions in tpcore routines since there is no longer a need to convert from moles to mass.

Note: Nested grid boundary conditions are set when tracer concentrations are in units of v/v dry air. This is acceptable since all values are converted to kg/kg total air just prior to advection.

An example of changes made within transport_mod, not including removal of residual mass correction, are shown in the following example from subroutine GEOS4_GEOS5_GLOBAL_ADV.

Legacy implementation 

     CALL TPCORE_FVDAS( ... )
      ... 
     CALL SET_FLOATING_PRESSURE( P_TP2 )

New implementation 

     ! Convert tracer concentration from mol/mol dry air to 
     ! kg/kg total air before advection (ewl, 6/15/15)
     CALL Convert_DryVV_to_DryKgKg( am_I_Root, N_TRACERS, 
    &                                Input_Opt, State_Chm, RC)
     CALL Convert_DryKgKg_to_MoistKgKg( am_I_Root, N_TRACERS,  
    &                                    State_Met, State_Chm, RC )

     CALL TPCORE_FVDAS( ... ) 

     ! Convert tracer concentration from kg/kg total air to 
     ! mol/mol dry air after advection (ewl, 6/15/15)
     CALL Convert_MoistKgKg_to_DryKgKg( am_I_Root, N_TRACERS,  
    &                                    State_Met, State_Chm, RC ) 
     CALL Convert_DryKgKg_to_DryVV( am_I_Root, N_TRACERS, 
    &                                Input_Opt, State_Chm, RC)
     ...
     CALL SET_FLOATING_PRESSURE( P_TP2 )
     CALL AIRQNT( am_I_Root, State_Met, RC )

--Lizzie Lundgren (talk) 15:50, 24 June 2015 (UTC)

Retrieved from "https://wiki.seas.harvard.edu/geos-chem/index.php?title=Air_Quantity_Updates_for_v11-01&oldid=29134"