Add new scheme Li2024 of ozone plant damage #2302

lifang0209 · 2023-12-28T18:02:06Z

Description of changes

Changes include:
OzoneMod.F90: add an option for scheme of Li2024:
Include add: (i) CalcOzoneUptakeLi2024OnePoint to calculate ozone uptake for a single point, and change the name of old module for Lombardozzi2015 and Falk schemes as CalcOzoneUptakeLFOnePoint. (ii) add CalcOzoneStressLi2024 and CalcOzoneStressLi2024OnePoint to calculate ozone stress.

Change CanopyFluxesMod.F90, OzoneBaseMod.F90, OzoneOffMod.F90 to read the variable sabv that CalcOzoneUptake use.
subroutine CalcOzoneUptake(this, bounds, num_exposedvegp, filter_exposedvegp, &
forc_pbot, forc_th, rssun, rssha, rb, ram, tlai, forc_o3, sabv)

In bld/namelist_files/namelist_definition_ctsm.xml: add stress_li2024 as the new option

Specific notes

Contributors other than yourself, if any: @slevis-lmwg @dlawrenncar @danicalombardozzi @ekluzek

CTSM Issues Fixed (include github issue #):

Are answers expected to change (and if so in what way)?
Yes, when using the new scheme

Any User Interface Changes (namelist or namelist defaults changes)?
new namelist option: stress_li2024

Testing performed, if any:
@lifang0209 ran I2000Sp and I2000BgcCrop

slevis-lmwg

@lifang0209 thank you for this nice update to the ozone code.

After we talked, I reviewed the code more carefully, posted a few comments, and requested a few improvements.
Scientifically the most pressing is to resolve the inconsistency between the formulas appearing in the code versus those documented in Li et al. (in prep).

slevis-lmwg · 2024-01-17T00:31:48Z

src/biogeophys/OzoneMod.F90

+     if (tlai > lai_thresh) then
+      ! o3 uptake decay
+      if (pftcon%evergreen(pft_type) == 1) then
+       leafturn = 1._r8/(pftcon%leaf_long(pft_type)*365._r8*24._r8)


Please replace 365 with dayspyr and 3600 with secsphr. I don't think that we have a parameter for 24, but please replace if we do.

Could we change leafturn to units of seconds^-1 by multiplying the denominator by secspyear (instead of by 365 * 24) and then get rid of dtimeh?

slevis-lmwg · 2024-01-17T00:48:30Z

src/biogeophys/OzoneMod.F90

+
+   ! calculate o3 flux per timestep
+     if(sabv > 0._r8)then  !daytime
+       o3fluxperdt = o3fluxcrit * dtime * 0.000001_r8


Could we avoid the unit conversion here by setting o3flux and o3_flux_threshold to the units that we need?

You already convert a precursor to o3flux in line 495, so you could embed this conversion there.

The latter (o3_flux_threshold) could be done in lines 511 to 523.

The units of o3flux, o3_flux_threshold, and o3uptake are commonly used, especially in observational community, so I don't want to change them

slevis-lmwg · 2024-01-17T00:58:21Z

src/biogeophys/OzoneMod.F90

+        o3coefv = max(0._r8, min(1._r8, 0.943_r8 * exp(-0.0085*o3uptake)))
+        o3coefg = max(0._r8, min(1._r8, 0.943_r8 * exp(-0.0058*o3uptake)))
+      end if
+      if(pft_type >= 9 .and. pft_type <= 11)then !Shrub


Please replace these numbers with is_shrub that you can find in pftconMod.F90.

Similarly change the other numbers in this section with parameters that already exist in pftconMod.F90.

slevis-lmwg · 2024-01-17T01:07:09Z

src/biogeophys/OzoneMod.F90

+       ! Determine parameter values for this pft
+      if(pft_type >= 1 .and. pft_type <= 3)then  !Needleleaf tree
+        o3coefv = max(0._r8, min(1._r8, 1.005_r8 - 0.0064_r8 * o3uptake))
+        o3coefg = max(0._r8, min(1._r8, 0.965_r8 * o3uptake ** (-0.041)))


These right-hand-sides do not match what I see in Li et al. (in prep), which says

1.014-0.010 POD0.5 0.975-0.037 ln(POD0.5)

formulas in the code are consistent with the Eqs.2-3 in our submitted MS that I sent you on Jan 17. You used an old MS version. The difference is in the data pre-processing which changed the sample of observations we used to generate the response functions shown in Eqs. 2-3. According to suggestion of ecophysiologist Felicity Hayes (our coauthor and the chair of LRTAP ozone impacts on vegetation in observations), in pre-processing, we conducted linear regression using the relative values and corresponding PODY for individual plant species in a study, and data with intercept falling outside the range of 0.9 to 1.1 are removed based on LRTAP convention to ensure a usable response relationship.

slevis-lmwg · 2024-01-17T01:11:11Z

src/biogeophys/OzoneMod.F90

+      end if
+      if(pft_type >= 9 .and. pft_type <= 11)then !Shrub
+        o3coefv = max(0._r8, min(1._r8, 1.000_r8-0.074_r8 * log(o3uptake)))
+        o3coefg = max(0._r8, min(1._r8, 0.991_r8-0.060_r8 * log(o3uptake)))


These right-hand-sides also do not match what I see in Li et al. (in prep), which says

1.142 POD5**(-0.140) 1.075 POD5**(-0.104)

see above reply

slevis-lmwg · 2024-01-17T01:12:42Z

src/biogeophys/OzoneMod.F90

+      end if
+      if (pft_type >= 15)then !Crop
+        o3coefv = max(0._r8, min(1._r8, 0.909_r8 - 0.028_r8 * log(o3uptake)))
+        o3coefg = max(0._r8, min(1._r8, 1.005_r8 - 0.169_r8 * tanh(o3uptake)))


Grasses and crops seem to match what I see in Li et al. (in prep).

see above reply

slevis-lmwg · 2024-01-17T01:19:12Z

src/biogeophys/PhotosynthesisMod.F90

               rs_z(p,iv) = min(1._r8/gs, rsmax0)
-               rs_z(p,iv) = rs_z(p,iv) / o3coefg(p)


@lifang0209 explained to me that moving the effect from the resistance to the conductance, also moves it to the current time step rather than the following time step. I have not checked the validity of that statement.

However @lifang0209 you said that you had tested and found the effect of this change on the simulation to be very small. Could you post some sort of documentation here (possibly a figure) that would make this point clear?

I have added comment in the PhotosynthesisMod.F90:
! Fang Li revised o3 impact on gs
! now: multiply gs response factor (o3coefg) with gs
! clm5.0: skipped the effect on gs and had rs divided by o3coefg,
! which led to a one timestep delay in the o3 impact on gs,
! although the impact on an annual scale is quite small

slevis-lmwg · 2024-01-17T01:27:44Z

src/biogeophys/OzoneMod.F90

+    !-----------------------------------------------------------------------
+
+    ! convert o3 from mol/mol to nmol m^-3
+    o3concnmolm3 = forc_ozone * 1.e9_r8 * (forc_pbot/(forc_th*SHR_CONST_RGAS*0.001_r8))


When you finalize this conversion (see my other comment), then replace any remaining hardwiring with parameter names. If the conversion parameters already exist elsewhere, please reuse them, and otherwise add them locally.

add as local variables

slevis-lmwg · 2024-01-17T01:35:13Z

src/biogeophys/OzoneMod.F90

+        if (pftcon%evergreen(pft_type) == 1) then
+         lai_thresh=0._r8 ! evergreens grow year-round
+       else  ! for deciduous vegetation
+        if(pft_type == 10)then !temperate shrub


As you did with "evergreen" above, please replace 10 with the shrub name from pftconMod.F90.

Similarly replace numbers below with pft names and use is_shrub, is_grass for these types in particular.

For 15 I think you can use nc3crop. All of these are in pftconMod.F90.

ekluzek · 2024-01-19T17:33:29Z

@lifang0209 thanks so much for your work here and this contribution! There is some more work to be done on this, so I'm converting it to a draft. That doesn't change the importance of this, it just marks it more clearly as there's more work to do here....

ekluzek

@lifang0209 thanks so much for the contribution. I have a few suggestions. Mostly that @slevis-lmwg and I will discuss and work on. I had one question for you about dtimeh though...

ekluzek · 2024-01-19T17:36:15Z

src/biogeophys/OzoneMod.F90

    real(r8) :: decay          ! o3uptake decay rate based on leaf lifetime (mmol m^-2)

-    character(len=*), parameter :: subname = 'CalcOzoneUptakeOnePoint'
+    ! o3:h2o resistance ratio defined by Sitch et al. 2007
+    real(r8), parameter :: ko3 = 1.67_r8


@lifang0209 should any of these be added to the parameter file, so they can be modified rather than hardcoded here?

Hi Erik, thanks for your time and suggestions.
The correct (observational) value for ko3 should be
real(r8), parameter :: ko3 = 1.51_r8
which used in CalcOzoneUptakeLi2024OnePoint in the same file (OzoneMod.F90).
I'm not sure if it's good to set different values of a parameter for different schemes in the parameter file.

ekluzek · 2024-01-19T17:41:26Z

src/biogeophys/OzoneMod.F90

@@ -494,21 +640,21 @@ subroutine CalcOzoneUptakeOnePoint( &
       endif

       if (pftcon%evergreen(pft_type) == 1) then
-          leafturn = 1._r8/(pftcon%leaf_long(pft_type)*365._r8*24._r8)
+          leafturn = 1._r8/(pftcon%leaf_long(pft_type)*avg_dayspyr*secspday)


I very much like using functions and constants to get things like the number of days in a year, and the seconds per day. That's an important improvement.

ekluzek · 2024-01-19T17:42:42Z

src/biogeophys/OzoneMod.F90

       else
          leafturn = 0._r8
       endif

       ! o3 uptake decay based on leaf lifetime for evergreen plants
-       decay = o3uptake * leafturn * dtimeh
+       decay = o3uptake * leafturn * dtime


Is this OK to do? Because dtimeh is in hours and dtime is in seconds, won't this change answers for the older versions?

According to Sam's suggestion, leafturn is changed from per hour to per second. Correspondingly, we should use dtime (in seconds) instead of dtimeh (in hours).

ekluzek · 2024-01-19T17:45:00Z

src/biogeophys/OzoneMod.F90

@@ -424,9 +444,125 @@ subroutine CalcOzoneUptake(this, bounds, num_exposedvegp, filter_exposedvegp, &
    end associate

  end subroutine CalcOzoneUptake
+    !-----------------------------------------------------------------------
+  subroutine CalcOzoneUptakeLi2024OnePoint( &


I think this should actually be spun off into it's own file. That would be better in keeping with the OO nature of how this is built. But, @slevis-lmwg and I will discuss, and I need to look at it more closely.

Thanks. That's not my area of expertise; I'll leave it to you and Sam to decide.
If you do this, don't forget move CalcOzoneUptakeLFOnePoint.
CalcOzoneUptakeLi2024OnePoint and CalcOzoneUptakeLFOnePoint are used to calculate the ozone uptake for a single point. The former is used for the new scheme (li2024), and the latter for the lombardozzi2015 and falk schemes.

lifang0209 · 2024-01-19T18:47:01Z

@lifang0209 thanks so much for the contribution. I have a few suggestions. Mostly that @slevis-lmwg and I will discuss and work on. I had one question for you about dtimeh though...

Add new scheme Li2024 of ozone plant damage

facbac6

slevis-lmwg assigned lifang0209, slevis-lmwg and ekluzek Dec 28, 2023

slevis-lmwg added PR status: work in progress labels Dec 28, 2023

Fang Li and others added 2 commits December 28, 2023 11:28

add comment about work to be done

3cb6d61

update Y in PODY and response functions in Li2024

47309ca

lifang0209 marked this pull request as ready for review January 15, 2024 01:58

slevis-lmwg requested changes Jan 17, 2024

View reviewed changes

Revise code according to Sam Levis's suggestions

61cc20f

lifang0209 requested a review from slevis-lmwg January 19, 2024 09:04

ekluzek marked this pull request as draft January 19, 2024 17:33

ekluzek requested changes Jan 19, 2024

View reviewed changes

lifang0209 closed this Jan 19, 2024

lifang0209 reopened this Jan 19, 2024

lifang0209 marked this pull request as ready for review January 19, 2024 19:03

ekluzek marked this pull request as draft January 20, 2024 07:06

ekluzek mentioned this pull request Jun 20, 2024

Update the drv_flds namelist definition file in CTSM to work with latest CAM/cmeps #2170

Closed

samsrabin added enhancement new capability or improved behavior of existing capability science Enhancement to or bug impacting science and removed enh - new science labels Aug 8, 2024

samsrabin removed the PR status: work in progress label Sep 12, 2024

ekluzek added this to the ctsm7.0.0 milestone Feb 4, 2025

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add new scheme Li2024 of ozone plant damage #2302

Add new scheme Li2024 of ozone plant damage #2302

lifang0209 commented Dec 28, 2023

slevis-lmwg left a comment

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024 •

edited

Loading

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

slevis-lmwg Jan 17, 2024

lifang0209 Jan 19, 2024

ekluzek commented Jan 19, 2024

ekluzek left a comment

ekluzek Jan 19, 2024

lifang0209 Jan 19, 2024

ekluzek Jan 19, 2024

lifang0209 Jan 19, 2024

ekluzek Jan 19, 2024

lifang0209 Jan 19, 2024

ekluzek Jan 19, 2024

lifang0209 Jan 19, 2024

lifang0209 commented Jan 19, 2024

		rs_z(p,iv) = min(1._r8/gs, rsmax0)
		rs_z(p,iv) = rs_z(p,iv) / o3coefg(p)

Add new scheme Li2024 of ozone plant damage #2302

Are you sure you want to change the base?

Add new scheme Li2024 of ozone plant damage #2302

Conversation

lifang0209 commented Dec 28, 2023

Description of changes

Specific notes

slevis-lmwg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

slevis-lmwg Jan 17, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ekluzek commented Jan 19, 2024

ekluzek left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lifang0209 commented Jan 19, 2024

slevis-lmwg Jan 17, 2024 •

edited

Loading