Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Handmerge feature/pcolarco/#1-gocart2g-updates into develop #68

Merged

Conversation

mathomp4
Copy link
Member

@mathomp4 mathomp4 commented Aug 4, 2021

This is a PR to merge in feature/pcolarco/#1-gocart2g-updates into develop. This was done as a hand-merge to be careful about it (rather than depending on GitHub).

Note, this possibly could affect the UFS work so I'm pinging @rmontuoro as this is what we want all future development based on. I think I kept things working so that your UFS code still works.

As far as I know how to test GOCART2G this is zero diff to feature/pcolarco/#1-gocart2g-updates. I did a c24 run with OPS GOCART2G in 2015 and...it was zero-diff?

For @pcolarco's benefit, develop "won" in the merge only in a few places as seen here:

https://github.com/GEOS-ESM/GOCART/compare/feature/pcolarco/%231-gocart2g-updates..feature/mathomp4/merge-develop-pcolarco-%231-gocart2g-updates

One main change that was on develop is that the scheme in dust is now specified by a string and not an integer. We can always undo this as needed, but it allowed for an easier merge. So:

emission_scheme: 1 # Ginoux

becomes:

emission_scheme: ginoux

The other main change is that in the Dust State Specs we went from:

 DU_SAND    | 1      | xy   | N   |      | sand_fraction
 DU_SILT    | 1      | xy   | N   |      | silt_fraction
 DU_CLAY    | 1      | xy   | N   |      | clay_fraction

to:

 DU_SAND    | 1      | xy   | N   |      | volume_fraction_of_sand_in_soil
 DU_SILT    | 1      | xy   | N   |      | volume_fraction_of_silt_in_soil
 DU_CLAY    | 1      | xy   | N   |      | volume_fraction_of_clay_in_soil

I think this is because @rmontuoro and UFS need a CF compliant standard name in the LONG NAME and not our...not so compilant names.

gmao-esherman and others added 30 commits March 9, 2021 09:07
removed diagnostic print statements
@mathomp4 mathomp4 added 0 diff The changes in this pull request have verified to be zero-diff with the target branch. Non 0-diff The changes in this pull request are non-zero-diff labels Aug 4, 2021
@mathomp4 mathomp4 self-assigned this Aug 4, 2021
@mathomp4 mathomp4 requested review from tclune, weiyuan-jiang and a team as code owners August 4, 2021 16:16
@mathomp4 mathomp4 requested a review from rmontuoro August 4, 2021 16:16
@mathomp4
Copy link
Member Author

mathomp4 commented Aug 4, 2021

@rmontuoro It would be good if you can test this branch with UFS and make sure nothing got broken on your end.

Copy link
Contributor

@rmontuoro rmontuoro left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These change will break UFS-Aerosols and require further code development on the NOAA side.

DU_CLAY | 1 | xy | N | scheme == 'fengsha' | volume_fraction_of_clay_in_soil
DU_SAND | 1 | xy | N | scheme == 'fengsha' | volume_fraction_of_sand_in_soil
DU_SILT | 1 | xy | N | scheme == 'fengsha' | volume_fraction_of_silt_in_soil
DU_SRC | 1 | xy | N | | erod - dust emissions
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These DU_* fields should be associated with a conditional statement (scheme == name) according to early discussions with @tclune.
Is there a particular reason why the COND column should be empty in this update?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro Well how do you do it if there are two schemes that need it? Both k14 and ginoux use du_src? I removed it because of that.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As for the others, DU_Z0, DU_GVF are k14 only, so those could be scheme == 'k14' protected.

But DU_SAND, DU_SILT, and DU_CLAY look to be used by both k14 and fengsha.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mathomp4 - thanks for the explanation. This is a question for @tclune. I believe something like:

scheme == 'ginoux' .or. scheme == 'k14'

could work.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, it does look like scheme == 'ginoux' .OR. scheme == 'k14' could work here. So maybe we can protect the IMPORT pointers that are affecting you a la:

 DU_SRC     | 1      | xy   | N   | scheme == 'ginoux' .OR. scheme == 'k14'     | erod - dust emissions
 DU_Z0      | 1      | xy   | N   | scheme == 'k14'     | aerodynamic_surface_roughness_for_aeolian_processes
 DU_GVF     | 1      | xy   | N   | scheme == 'k14'     | GVF
 DU_SAND    | 1      | xy   | N   | scheme == 'fengsha' .OR. scheme == 'k14'     | volume_fraction_of_sand_in_soil
 DU_SILT    | 1      | xy   | N   | scheme == 'fengsha' .OR. scheme == 'k14'     | volume_fraction_of_silt_in_soil
 DU_CLAY    | 1      | xy   | N   | scheme == 'fengsha' .OR. scheme == 'k14'     | volume_fraction_of_clay_in_soil
 DU_RDRAG   | m-1    | xy   | N   | scheme == 'fengsha' | drag_partition
 DU_SSM     | 1      | xy   | N   | scheme == 'fengsha' | sediment_supply_map
 DU_UTHRES  | m s-1  | xy   | N   | scheme == 'fengsha' | surface_dry_threshold_velocity
 FRSNOW     | 1      | xy   | N   | scheme == 'fengsha' | surface_snow_area_fraction
 SLC        | 1      | xy   | N   | scheme == 'fengsha' | liquid_water_content_of_soil_layer
 DU_TEXTURE | 1      | xy   | N   | scheme == 'k14'     | soil_texture
 DU_VEG     | 1      | xy   | N   | scheme == 'k14'     | vegetation_type
...
 U10N       | m s-1  | xy   | N   | scheme == 'k14'     | equivalent_neutral_10-meter_eastward_wind
 V10N       | m s-1  | xy   | N   | scheme == 'k14'     | equivalent_neutral_10-meter_northward_wind
...
 WCSF       | m3 m-3 | xy   | N   | scheme == 'k14'     | water_surface_layer

Though TROPP is harder as it is scheme agnostic.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's doable:

category: IMPORT
#----------------------------------------------------------------------------------------
#  VARIABLE            | DIMENSIONS  |          Additional Metadata
#----------------------------------------------------------------------------------------
     NAME   | UNITS    | DIMS | VLOC | COND                                     | LONG NAME
#----------------------------------------------------------------------------------------
 DU_SRC     | 1        | xy   | N    | scheme == 'ginoux' .OR. scheme == 'k14'  | erod - dust emissions
 DU_Z0      | 1        | xy   | N    | scheme == 'k14'                          | aerodynamic_surface_roughness_for_aeolian_processes
 DU_GVF     | 1        | xy   | N    | scheme == 'k14'                          | GVF
 DU_SAND    | 1        | xy   | N    | scheme == 'fengsha' .OR. scheme == 'k14' | volume_fraction_of_sand_in_soil
 DU_SILT    | 1        | xy   | N    | scheme == 'fengsha' .OR. scheme == 'k14' | volume_fraction_of_silt_in_soil
...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro I can supply you with this file if you like so you can test with it before I push to this branch?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mathomp4 - sure.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro here you go. I had to add ".txt" at the end so GitHub would upload it, so you'll have to un-do that.

DU2G_StateSpecs.rc.txt

Copy link
Contributor

@rmontuoro rmontuoro Aug 6, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @mathomp4. I've built and run this update within the UFS as a preliminary test after some integration work, and suggested a few changes in my latest code review.

I would also like to add the changes in my PR #65 as it looks like they are still relevant.

WET1 | 1 | xy | N | | surface_soil_wetness
LWI | 1 | xy | N | | land-ocean-ice_mask
FRLAKE | 1 | xy | N | | fraction_of_lake
TROPP | Pa | xy | N | | tropopause_pressure_based_on_blended_estimate
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These new imports are concerning. They will break the UFS model and require further development work that cannot occur at this time, since we are just starting our evaluation runs.

@pcolarco - could you please clarify how satisfy/compute the following imports?

  • TROPP
  • U10N/V10N
  • WCSF

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Of these, U10N, V10N, and WCSF seem to be K14 related.

TROPP though, seems to be GEOS endemic and it seems to come from tropovars in https://github.com/GEOS-ESM/GMAO_Shared/blob/main/GEOS_Shared/tropovars.F90. I think TropP3 is eventually TROPP_BLENDED

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mathomp4 - Thank you, this is helpful. We could implement tropovars in UFS-Aerosols to compute TROPP inline, but it requires potential vorticity as an input and this field is not available at this time.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm. Well, you are beyond my ken. Maybe @amdasilva or someone with science knowledge could help come up with a tropopause pressure that could be used as an approximate proxy? (Also, I have no idea if that is actually where TROPP_BLENDED comes from, but I think it might from my search.)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand it correctly, a "blended" tropopause pressure is determined "blending" the tropopause k-levels estimated by looking at temperature (ktherm) and potential vorticity (kept) profiles:

kblend = max( ktherm,kepv )

We could use other methods to compute the pressure at the tropopause that do not require the potential vorticity (see, for example, Reichler et al., 2003). However, I don't know how they would compare with the blended method and this seems to introduce an uncomfortably high degree of uncertainty at this stage.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely needs a comment from @amdasilva

Copy link
Member Author

@mathomp4 mathomp4 Aug 4, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro There is another possibility: maybe you don't care? If I look at Process_Library/GOCART2G_Process.F90, tropp is only used to calculate stexttau and stscatau, e.g.,:

if( present(stexttau) .and. associated(stexttau) ) then
if (ple(i,j,k) .le. tropp(i,j)) then
stexttau(i,j,w) = stexttau(i,j,w) + tau
elseif(ple(i,j,k) .gt. tropp(i,j) .and. ple(i,j,k-1) .lt. tropp(i,j)) then
stexttau(i,j,w) = stexttau(i,j,w) + log(tropp(i,j)/ple(i,j,k-1))/log(ple(i,j,k)/ple(i,j,k-1))*tau
endif
endif

So, if UFS doesn't need to know about "stratospheric ext. AOT at 550 nm" or "stratospheric sct. AOT at 550 nm" maybe it's a throwaway? Just fill it with something that won't cause a "divide-by-zero" and move along?

Copy link
Contributor

@rmontuoro rmontuoro Aug 4, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mathomp4 - good catch! It is hard to make a decision on this right now since we are still discussing how aerosol-radiation feedback will be introduced in the UFS using GOCART. We might need to provide stexttau to the radiation scheme running in the atmospheric model. Also, I would like to preserve the ability to compute these diagnostics in UFS-Aerosols. Nevertheless, your suggestion is good and could be a starting point.
This also means that the accuracy of TROPP may not play such a big role here.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One other option would be to add a procedure computing the atmospheric pressure at the tropopause to the process library. This would help achieve consistency across different applications.

@@ -26,7 +21,7 @@ hydrophobic_fraction: 0.5
fscav: 0.0 0.4

# particle radius
particle_radius_microns: 0.0212 0.0212
particle_radius_microns: 0.35 0.35
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some of the shifts in parameter values are disconcerting. This one in particular which is > 10x.

tclune
tclune previously approved these changes Aug 5, 2021
Copy link
Contributor

@tclune tclune left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, its different. Hmm.

@mathomp4
Copy link
Member Author

mathomp4 commented Aug 5, 2021

Okay. I pushed up the latest Dust State Specs with the scheme logic as it seems to be zero-diff for ginoux and k14. I also added back a print I mistakenly removed along with an assert in SetServices (to make sure you supply a supported scheme), and put in the three fengsha settings that were in develop.

Of course, fengsha fails. I tried it with the ginoux settings and you of course get:

        EXTDATA: INFO: ExtData could not satisfy import DU_RDRAG
        EXTDATA: INFO: ExtData could not satisfy import DU_SSM
        EXTDATA: INFO: ExtData could not satisfy import DU_UTHRES
        EXTDATA: INFO: ExtData could not satisfy import FRSNOW
        EXTDATA: INFO: ExtData could not satisfy import SLC
pe=00080 FAIL at line=00851    MAPL_ExtDataGridCompMod.F90              <Found   5 unfullfilled imports in extdat>

because we don't have all the necessary imports in GEOS yet. (And I need to fix that error message in MAPL.)

Copy link
Contributor

@rmontuoro rmontuoro left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Additional changes are required to make this update work with the UFS.

@@ -150,8 +147,20 @@ subroutine SetServices (GC, RC)
call ESMF_ConfigGetAttribute (cfg, self%Ch_DU_res, label='Ch_DU:', __RC__)
call ESMF_ConfigGetAttribute (cfg, self%rlow, label='radius_lower:', __RC__)
call ESMF_ConfigGetAttribute (cfg, self%rup, label='radius_upper:', __RC__)
call ESMF_ConfigGetAttribute (cfg, self%clayFlag, label='clayFlag:', __RC__)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The statements on line 150-153 need to be guarded as they are only required for the k14 dust scheme.
I suggest consolidating the _ASSERT on line 159 in a select case block. See patch below:

--- a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90
+++ b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90
@@ -147,10 +147,6 @@ contains
     call ESMF_ConfigGetAttribute (cfg, self%Ch_DU_res,  label='Ch_DU:', __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%rlow,       label='radius_lower:', __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%rup,        label='radius_upper:', __RC__)
-    call ESMF_ConfigGetAttribute (cfg, self%clayFlag,   label='clayFlag:', __RC__)
-    call ESMF_ConfigGetAttribute (cfg, self%f_swc,      label='soil_moisture_factor:', __RC__)
-    call ESMF_ConfigGetAttribute (cfg, self%f_scl,      label='soil_clay_factor:', __RC__)
-    call ESMF_ConfigGetAttribute (cfg, self%uts_gamma,  label='uts_gamma:', __RC__)
 
     call ESMF_ConfigGetAttribute (cfg, emission_scheme, label='emission_scheme:', default='ginoux', __RC__)
     self%emission_scheme = ESMF_UtilStringLowerCase(trim(emission_scheme), __RC__)
@@ -171,11 +167,21 @@ contains
 
 !   read FENGSHA-specific parameters
 !   --------------------------------
-    if (self%emission_scheme == 'fengsha') then
-      call ESMF_ConfigGetAttribute (cfg, self%alpha,  label='alpha:', __RC__)
-      call ESMF_ConfigGetAttribute (cfg, self%gamma,  label='gamma:', __RC__)
-      call ESMF_ConfigGetAttribute (cfg, self%kvhmax, label='vertical_to_horizontal_flux_ratio_limit:', __RC__)
-    end if
+    select case (self%emission_scheme)
+      case ('fengsha')
+        call ESMF_ConfigGetAttribute (cfg, self%alpha,  label='alpha:', __RC__)
+        call ESMF_ConfigGetAttribute (cfg, self%gamma,  label='gamma:', __RC__)
+        call ESMF_ConfigGetAttribute (cfg, self%kvhmax, label='vertical_to_horizontal_flux_ratio_limit:', __RC__)
+      case ('k14')
+        call ESMF_ConfigGetAttribute (cfg, self%clayFlag,   label='clayFlag:', __RC__)
+        call ESMF_ConfigGetAttribute (cfg, self%f_swc,      label='soil_moisture_factor:', __RC__)
+        call ESMF_ConfigGetAttribute (cfg, self%f_scl,      label='soil_clay_factor:', __RC__)
+        call ESMF_ConfigGetAttribute (cfg, self%uts_gamma,  label='uts_gamma:', __RC__)
+      case ('ginoux')
+        ! nothing to do
+      case default 
+      _ASSERT_RC(.false., "Unallowed emission scheme: "//trim(self%emission_scheme)//". Allowed: ginoux, k14, fengsha", ESMF_RC_NOT_IMPL)
+    end select

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro Done and tested as zero-diff to ginoux and k14. Please re-test with fengsha

ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_StateSpecs.rc Outdated Show resolved Hide resolved
@mathomp4
Copy link
Member Author

mathomp4 commented Aug 6, 2021

This looks like we are close. If @rmontuoro can agree with this PR. We can pull it in as well as #71 (which doesn't affect ginoux or k14) and then make our pre-release.

@mathomp4 mathomp4 requested a review from rmontuoro August 6, 2021 13:34
Copy link
Contributor

@rmontuoro rmontuoro left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This revision runs with UFS but changes are recommended.

@@ -2392,23 +3549,28 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch
end if !present(extcoef)...

if( (present(exttau) .and. associated(exttau)) .or. &
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@tclune - is this valid Fortran?

if (present(a) .and. associated(a)) then

Should this be:

if (present(a)) then
   if (associated(a)) then

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro I talked with @tclune and he agrees this is an issue. For now we will approve (since Intel is at least "working") and merge what is here but I'll make an issue to fix this. I don't want to try and untangle this:

if( (present(exttau) .and. associated(exttau)) .or. &
(present(stexttau) .and. associated(stexttau)) .or. &
(present(scatau) .and. associated(scatau)) .or. &
(present(stscatau) .and. associated(stscatau)) .or. &
(present(exttau25) .and. associated(exttau25)) .or. &
(present(exttaufm) .and. associated(exttaufm)) .or. &
(present(scatau25) .and. associated(scatau25)) .or. &
(present(scataufm) .and. associated(scataufm)) ) then

on a Friday afternoon!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Issue #72 created

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, this is not quite as bad as I first thought. Can be converted to 2 conditionals:

if (present(exttau) .and. present(stexttau) .and. ...) then
   if (associated(exttau) .and. associated(stexttau) .and. ...) then
       ...


! Set wavelengths in universal config

call MAPL_ConfigSetAttribute (cf, self%wavelengths_profile, label='wavelengths_for_profile_aop_in_nm:', __RC__)
call MAPL_ConfigSetAttribute (cf, self%wavelengths_vertint, label='wavelengths_for_vertically_integrated_aop_in_nm:', __RC__)
call MAPL_ConfigSetAttribute (cf, wavelengths_diagmie, label='aerosol_monochromatic_optics_wavelength:', __RC__)
call MAPL_ConfigSetAttribute (cf, wavelengths_diagmie, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

wavelengths_diagmie is still declared as integer:

integer, allocatable, dimension(:) :: wavelengths_diagmie

which is inconsistent with self%wavelengths_profile and self%wavelengths_vertint (real).

I would still recommend including changes from PR #65 into this branch to address this issue and streamline the code. For example, conversion from nm to m should be done only once (see https://github.com/GEOS-ESM/GOCART/pull/65/files#diff-682724d291af0ad213c4e1b2ffa97156de4989543290b60bf362adfa98e3d8da):

! Convert input wavelengths from nm to m for internal use
self%wavelengths_profile = 1.e-09 * self%wavelengths_profile
self%wavelengths_vertint = 1.e-09 * self%wavelengths_vertint
wavelengths_diagmie = 1.e-09 * wavelengths_diagmie
! Set wavelengths in universal config

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmontuoro I guess they are declared integers because they are integers in the file:

aerosol_monochromatic_optics_wavelength_in_nm_from_LUT: 470 550 670 870

But as you said, the others aren't. For now, we might accept this as is and then you or I can make a (new) PR to fix it once we have a good "baseline" to compare to.

@mathomp4 mathomp4 merged commit e392831 into develop Aug 6, 2021
@mathomp4 mathomp4 deleted the feature/mathomp4/merge-develop-pcolarco-#1-gocart2g-updates branch August 6, 2021 15:54
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
0 diff The changes in this pull request have verified to be zero-diff with the target branch. Non 0-diff The changes in this pull request are non-zero-diff
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants