Get Deposition Resistance Outputs

kiarash_farzad · February 22, 2024, 9:02pm

Hi,

I was wondering if there is any quick way to get outputs from the following or similar variables from CMAQ:

     Real            :: Ra      ! aerodynamic resistance (s/m)
     Real            :: Rb      ! quasi-laminar soil/water resistance (s/m)
     Real            :: Rb_leaf ! quasi-laminary leaf resistance (s/m)
     Real            :: Rst     ! Stomatal resistance (s/m)
     Real            :: Rcut    ! cuticlular resistance (s/m)
     Real            :: Rgc     ! under canopy soil resistance (s/m)
     Real            :: Rg      ! soil resistance (s/m)
     Real            :: Rwat    ! surface water resistance (s/m)

Thanks,
Kiarash

hogrefe.christian · February 23, 2024, 12:23pm

Hi Kiarash,

no, to my knowledge there is no quick way to output these component resistances used in the deposition velocity calculations. Doing so would require modifications to the M3DRY or STAGE dry deposition code, depending on which option you are using - I am guessing STAGE due to the Rb_leaf term. Please note that in STAGE all of these component resistances are being calculated for each land use (LU) type for use in the calculation of LU specific deposition velocities and fluxes, and that the resulting LU specific deposition velocities and fluxes can be output by setting CTM_MOSAIC to Y. If you wanted to make code modifications to output these component resistances for each LU type and species, keep in mind that the resulting files would become quite large.

Also note that the grid-scale (not LU specific) inverted aerodynamic resistance RADYNI is available from MCIP and in CMAQv5.4 STAGE is used to re-normalize the LU specific Ra terms such that the LU weighted aggregate of the LU specific Ra terms matches the grid-scale values from MCIP (see subroutine RA_WRF in stage/MOSAIC_MOD/F).

Christian

jbash · February 23, 2024, 1:08pm

Hi Kiarash,

When using the STAGE option in CMAQ v5.4, you can get land use specific, Ra, Rst, LAI, Ustar, Z0, and Vegetation fraction by setting the CTM_MOSAIC environmental variable to Y in the run script. It is possible to add other variables to this output by adding them the the MOSAIC_Type derived data type and allocating them in ASX_DATA_MOD.F. These can then be populated in STAGE_MOD.F and added to the write statements in the STAGE_OUTPUT.F modules. If you wish to do this, I will be happy to help.

Take care,
Jesse

kiarash_farzad · February 29, 2024, 2:23pm

Hi Jesse,

I am getting back to you to see whether I can make the changes or not.

First I have made the following changes:

STAGE_OUTPUT.F:

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RA',
     &                            DATE, TIME, MOSAIC_DATA%RA ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, MOSAIC_DATA%RB ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, MOSAIC_DATA%RB_L ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RST',
     &                            DATE, TIME, Mosaic_Data%RSTW ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, Mosaic_Data%RCUT ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, Mosaic_Data%RGC ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, Mosaic_Data%RG ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, '??',
     &                            DATE, TIME, Mosaic_Data%RWAT ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

ASX_DATA_MOD.F:

!> Sub grid cell resistances
         Real, Allocatable :: RA    ( :,:,: )    ! aerodynamic resistance [s/m]
         Real, Allocatable :: RB    ( :,:,: )    ! quasi-laminar soil/water resistance [s/m]
         Real, Allocatable :: RB_L  ( :,:,: )    ! quasi-laminary leaf resistance [s/m]
         Real, Allocatable :: RSTW  ( :,:,: )    ! Stomatal Resistance of water [s/m]
         Real, Allocatable :: RCUT  ( :,:,: )    ! cuticlular resistance [s/m]
         Real, Allocatable :: RGC   ( :,:,: )    ! under canopy soil resistance [s/m]
         Real, Allocatable :: RG    ( :,:,: )    ! soil resistance [s/m]
         Real, Allocatable :: RWAT  ( :,:,: )    ! surface water resistance [s/m]
      End Type MOSAIC_Type

I have two questions here:
1- Where should I populate these variables, I think they get populated them selves.
2- Will there be any changes in MOSAIC_MOD.F?
I also wage on “MOSAIC_Type” part you mentioned.

I appreciate if you could help me out with it.

Thanks,
Kiarash

jbash · March 6, 2024, 2:42pm

Hi Kiarash,

There should be no changes to MOSAIC_MOD.F. The best area to populate the variables is in STAGE_MOD.F where they are used. The resistances are split between water and land RB will have to be set for both conditions. The code should look something like.
In the if ( Water( l ) )) conditional:

Mosaic_data%RWAT( c,r,l ) = Rwat
Mosaic_data%RB( c,r,l ) = Rb

In the if ( lai .Gt. 0.0 ) conditional
Mosaic_data%RB_L( c,r,l ) = Rb_leaf
Mosaic_data%RCUT( c,r,l ) = Rcut

After the if( lai .Gt. 0.0) conditional
Mosaic_data%RG( c,r,l ) = Rg
Mosaic_data%RGC( c,r,l ) = Rgc

This should populate your variables and they should be saved in ASX_DATA_MOD.F and available to write to disk in STAGE_OUTPUT.F. Make sure after allocating the new arrays to initialize to 0.0 as these values are only calculated if need based on the land cover.

Take care,
Jesse

kiarash_farzad · March 12, 2024, 7:50pm

Hi Jesse,

Thanks a lot for your guidance.
I was able to successfully output the resistances after a few tries.
Appreciate your help a lot!

I am sharing the changes that I made here, so that if anyone else in the community wants, they can just use it.

CMAQv5.4
ASX_DATA_MOD.F (3 changes):

under “Sub grid cell resistances”, line 190
under “Allocation”, line 548
just after step (2), initiate with 0.0

(1)
!> Sub grid cell resistances
         Real, Allocatable :: RA    ( :,:,: )    ! aerodynamic resistance [s/m]
         Real, Allocatable :: RB    ( :,:,: )    ! quasi-laminar soil/water resistance [s/m]
         Real, Allocatable :: RB_L  ( :,:,: )    ! quasi-laminary leaf resistance [s/m]
         Real, Allocatable :: RSTW  ( :,:,: )    ! Stomatal Resistance of water [s/m]
         Real, Allocatable :: RCUT  ( :,:,: )    ! cuticlular resistance [s/m]
         Real, Allocatable :: RGC   ( :,:,: )    ! under canopy soil resistance [s/m]
         Real, Allocatable :: RG    ( :,:,: )    ! soil resistance [s/m]
         Real, Allocatable :: RWAT  ( :,:,: )    ! surface water resistance [s/m]

(2)
         ALLOCATE( Mosaic_Data%USTAR   ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%LAI     ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%DELTA   ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%VEG     ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%Z0      ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RA      ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RB      ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RB_L      ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RSTW    ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RCUT    ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RGC    ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RG    ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%RWAT    ( NCOLS,NROWS,Tile_Data%n_lufrac ),
     &             Mosaic_Data%NAME    ( Tile_Data%n_lufrac ),
     &             Mosaic_Data%LU_Type ( Tile_Data%n_lufrac ),
     &             STAT = ALLOCSTAT )

(3)
         Mosaic_Data%USTAR   = 0.0
         Mosaic_Data%LAI     = 0.0
         Mosaic_Data%DELTA   = 0.0
         Mosaic_Data%VEG     = 0.0
         Mosaic_Data%Z0      = 0.0
         Mosaic_Data%RSTW    = 1.0e6
         Mosaic_Data%RA      = 1.0e6
         Mosaic_Data%RB      = 1.0e6
         Mosaic_Data%RB_L    = 0.0
         Mosaic_Data%RCUT    = 0.0
         Mosaic_Data%RGC     = 0.0
         Mosaic_Data%RG      = 0.0
         Mosaic_Data%RWAT    = 0.0
         Mosaic_Data%NAME    = Tile_Data%name_lu
         Mosaic_Data%LU_Type = Tile_Data%cat_lu

STAGE_MOD.F:

!-------------------------------------------------------------------------------------------------
! Quasi-laminar boundary layer resistance
!-------------------------------------------------------------------------------------------------
               Rb  = Calc_Rbw( ustar, scc_pr_23 )

!-------------------------------------------------------------------------------------------------
! Air-water exchange
!-------------------------------------------------------------------------------------------------
               If( Water( l ) ) Then
!-------------------------------------------------------------------------------------------------
! Resistance to surface water
!-------------------------------------------------------------------------------------------------
                  Rwat = Calc_Rwater( ustar, q_2m, temp_2m, temp_g, tw,  LeBasM( s ), heff_wat, O3_hit,
     &                                Hg_hit, sea_ice  )

                  If(c_wat .Gt. 0.0 ) Then
                     LU_Emis = c_wat / ( Ra + Rb + Rwat )
                  End If

                  Tile_Data%Lu_Vd( c,r,n,l ) = 1.0 / ( Ra + Rb + Rwat )
                  f_wat  = f_wat  + frac_lu * ( c_wat - c_atm ) / ( Ra + Rb + Rwat )

                  Mosaic_Data%RWAT( c,r,l ) = Rwat
                  Mosaic_Data%RB( c,r,l ) = Rb

               Else
!-------------------------------------------------------------------------------------------------
! Air-land exchange
!-------------------------------------------------------------------------------------------------
! Resistance to air-canopy exchange
!-------------------------------------------------------------------------------------------------
C Calculate Rst
                  If( lai .Gt. 0.0 ) Then
!-------------------------------------------------------------------------------------------------
! Quazi Laminar Resistance to leaf
!-------------------------------------------------------------------------------------------------
                     Rb_leaf = Calc_Rb_leaf( k_vis, dif_T, ustar, l_leaf, lai )

!-------------------------------------------------------------------------------------------------
! Resistance to air-stomatal exchange
!-------------------------------------------------------------------------------------------------
                     Rst = Calc_Rst( Mosaic_Data%RSTW( c,r,l ), dwat_T, dif_T, heff_ap, f0( s ), lai )

!-------------------------------------------------------------------------------------------------
! Resistance to air-cuticle exchange
!-------------------------------------------------------------------------------------------------
                     Rcut = Calc_Rcut( temp_g, rel_rx(s), molwt_all( s ), M_ac( s ), heff, dif_T,
     &                            a_cut( l ), snow, no_snow, dry, wet, rh, lai, O3_hit, NH3_hit, ABFLUX )

                     Mosaic_Data%RCUT( c,r,l ) = Rcut
                     Mosaic_Data%RB_L( c,r,l ) = Rb_leaf

!-------------------------------------------------------------------------------------------------
! Resistance to air-canopy covered soil exchange
!-------------------------------------------------------------------------------------------------

                  Else ! LAI = 0.0
                     Rst     = 1.0e6
                     Rcut    = 1.0e6
                     Rb_leaf = 1.0e6
                  End If ! LAI

!-------------------------------------------------------------------------------------------------
! Resistance to air-base soil exchange
!-------------------------------------------------------------------------------------------------
                  If( ABFLUX .And. NH3_Hit ) Then
                     Call Get_NH3_Comp( c_stom, c_ll, c_grnd, dif_T, r, c, l, l_ag )                                                                       ! Qauzi Laminar Boundary Layer resistance
                  End If ! ABFLUX and NH3

                  Rgc = Calc_Rg( dif_T, k_vis, rel_rx( s ), Ra, ustar, lai, temp_g, molwt_all( s ), M_ac( s ),
     &                           heff, sm_v1cm,sm_v5cm, sm_vsat, sm_vfc, sm_vwlt, sm_vres, sm_bslp, zsoil1,
     &                           ir_frac, dry, wet, snow, no_snow, O3_Hit, NH3_Hit, ABFLUX )

                  Rg  = Calc_Rg( dif_T, k_vis, rel_rx( s ), Ra, ustar, 0.0, temp_g, molwt_all( s ), M_ac( s ),
     &                           heff, sm_v1cm,sm_v5cm, sm_vsat, sm_vfc, sm_vwlt, sm_vres, sm_bslp, zsoil1,
     &                           ir_frac, dry, wet, snow, no_snow, O3_Hit, NH3_Hit, ABFLUX )

                  Mosaic_Data%RGC( c,r,l ) = Rgc
                  Mosaic_Data%RG( c,r,l ) = Rg

STAGE_OUTPUT.F (3 changes):

increase number of NVARS3D from + 7 to + 13
under “Set output file characteristics based on GRIDDESC and open the dep velocity dignostic file” line 146
under “STAGE Diagnostic Output” line 258

(1)
NVARS3D = Tile_Data%n_Vd + 13

(2)
N = 0
            N = N + 1 ! 1
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RA'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'aerodynamic resistance for land use category'

            N = N + 1 ! 2
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'LUFRAC'
            UNITS3D( N ) = 'dimensionless'
            VDESC3D( N ) = 'Fractional land use'

            N = N + 1 ! 3
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'LAI'
            UNITS3D( N ) = 'dimensionless'
            VDESC3D( N ) = 'leaf area index for land use category'

            N = N + 1 ! 4
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'USTAR'
            UNITS3D( N ) = 'm s-1'
            VDESC3D( N ) = 'friction velocity for land use category'

            N = N + 1 ! 5
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'Z0'
            UNITS3D( N ) = 'm'
            VDESC3D( N ) = 'surface roughness for land use category'

            N = N + 1 ! 6
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RST'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'Stomatal resistance to water vapor'

            N = N + 1 ! 7
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'VEG'
            UNITS3D( N ) = 'ratio'
            VDESC3D( N ) = 'Vegetation coverage'

            N = N + 1 ! 8
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RB'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'quasi-laminar soil/water resistance'

            N = N + 1 ! 9
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RB_L'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'quasi-laminary leaf resistance'

            N = N + 1 ! 10
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RCUT'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'cuticlular resistance'

            N = N + 1 ! 11
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RGC'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'under canopy soil resistance'

            N = N + 1 ! 12
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RG'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'soil resistance'

            N = N + 1 ! 13
            VTYPE3D( N ) = M3REAL
            VNAME3D( N ) = 'RWAT'
            UNITS3D( N ) = 's m-1'
            VDESC3D( N ) = 'surface water resistance'

(3)
               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'LUFRAC',
     &                            DATE, TIME, Tile_Data%LUFRAC ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF


               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RA',
     &                            DATE, TIME, MOSAIC_DATA%RA ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RB',
     &                            DATE, TIME, Mosaic_Data%RB ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RB_L',
     &                            DATE, TIME, Mosaic_Data%RB_L ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RST',
     &                            DATE, TIME, Mosaic_Data%RSTW ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RCUT',
     &                            DATE, TIME, Mosaic_Data%RCUT ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RGC',
     &                            DATE, TIME, Mosaic_Data%RGC ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RG',
     &                            DATE, TIME, Mosaic_Data%RG ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

               IF ( .NOT. WRITE3( CTM_DEPV_MOS, 'RWAT',
     &                            DATE, TIME, Mosaic_Data%RWAT ) ) THEN
                   XMSG = 'Could not write ' // CTM_DEPV_MOS // ' file'
                  CALL M3EXIT ( PNAME, DATE, TIME, XMSG, XSTAT1 )
               END IF

Best,
Kiarash

Topic		Replies	Views
About HNO3 dry deposition velocity Chemistry and Aerosols cmaq	2	359	April 5, 2023
About particle deposition scheme used in CMAQ CMAQ	2	408	October 19, 2020
How to modify ozone dry deposition velocity? CMAQ	1	241	September 30, 2023
Finding CMAQ species dry/wet deposition files in the Data Warehouse Data Warehouse cmaq	1	344	August 23, 2021
Can CMAQ output Jvalue of NO2? Meteorology	1	264	April 18, 2023

Get Deposition Resistance Outputs

Related Topics