fix missing or mis-matched variable types in nst_water_prop #126

DeniseWorthen · 2023-11-12T16:04:27Z

fixes missing type-kind in module_nst_water_prop.f90 #125
fixes variable type mis-matches such as

if(z>0) then
       zgamma=z/gamma
       f_aw=(f/z)*((gamma/z)*(1-exp(-zgamma))-exp(-zgamma))

by adding zero and one as parameter(kind_phys) variables.

Additional cleanup of various nst-related source code:

remove column 6 continuation lines in sfc_nst*f files and rename as .f90 files
remove unused variables and make only needed variables, routines or procedures public
move use statements to module level and remove use statements or implict none from included subroutines

For #125, initial testing found two non-reproducing baselines with this change. However retesting at 6397387 has shown all UWM baselines are B4B (see ufs-community/ufs-weather-model#1993 for logs)

* add one,zero and half * fix instances of reals compared to integer and integers used in real expressions

* fix type mis-match in call to int_epn using parameter zero in module_nst_model

XuLi-NOAA

The changes are fine to me.

DeniseWorthen · 2023-11-13T19:17:55Z

@XuLi-NOAA @Qingfu-Liu I'm still finding issues w/ some of the nst modules. I don't know if either of you can explain. In sfc_nst.f, one of the input arguments is sinlat=GFS_Data(cdata%blk_no)%Grid%sinlat.

It is used here as an argument to the function grv in module_nst_water_prop.f90:

https://github.com/DeniseWorthen/ccpp-physics/blob/c2ec9e5e165dac07a991424bd22bb0667cf32017/physics/sfc_nst.f#L372

The function is

function grv(lat)
  real(kind=kind_phys) :: lat
  real(kind=kind_phys) :: gamma,c1,c2,c3,c4,pi,phi,x
  gamma=9.7803267715
  c1=0.0052790414
  c2=0.0000232718
  c3=0.0000001262
  c4=0.0000000007
  pi=3.141593


  phi=lat*pi/180
  x=sin(phi)
  grv=gamma*(1+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
  !print *,'grav=',grv,lat
end function grv

Does this seem right? It appears as the function grv is expecting latitude, not sinlat, since it is converting it to radians and then taking the sine.

DeniseWorthen · 2023-11-13T22:20:14Z

I've found the reference to gravity forumla 1980. If sinlat is the argument passed in to sfc_nst.f, it seems clear that the grv function should be:

real(kind_phys) function grv(x)
  real(kind=kind_phys) :: x  ! sinlat
  real(kind=kind_phys) :: gamma,c1,c2,c3,c4
  gamma=9.7803267715
  c1=0.0052790414
  c2=0.0000232718
  c3=0.0000001262
  c4=0.0000000007

  grv=gamma*(one+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
  !print *,'grav=',grv,lat
end function grv

dustinswales · 2023-11-13T22:49:29Z

physics/module_nst_water_prop.f90

@@ -12,6 +12,8 @@ module module_nst_water_prop
  public :: rhocoef,density,sw_rad,sw_rad_aw,sw_rad_sum,sw_rad_upper,sw_rad_upper_aw,sw_rad_skin,grv,solar_time_from_julian,compjd, &
            sw_ps_9b,sw_ps_9b_aw,get_dtzm_point,get_dtzm_2d

+  integer, parameter :: kp = kind_phys


I know this probably predates your work, but why use two fields for the same thing and reference both throughout the code?
My suggestion(s) would be to modify the include statement:

use machine, only : kind_phys

use machine, only : kp =>kind_phys
and use kp in the code.
OR replace all instances of kp with kind_phys.

Yes, I can definitely do that. I prefer use machine, only : kp =>kind_phys and then using as you say only kp. Mostly I was following what seemed typical style since I don't work on the ccpp side much.

dustinswales · 2023-11-13T22:54:08Z

I've found the reference to gravity forumla 1980. If sinlat is the argument passed in to sfc_nst.f, it seems clear that the grv function should be:
real(kind_phys) function grv(x)
  real(kind=kind_phys) :: x  ! sinlat
  real(kind=kind_phys) :: gamma,c1,c2,c3,c4
  gamma=9.7803267715
  c1=0.0052790414
  c2=0.0000232718
  c3=0.0000001262
  c4=0.0000000007

  grv=gamma*(one+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
  !print *,'grav=',grv,lat
end function grv

This surely looks correct to me.
I'm sure this changes the results though(?), hopefully for the better.

XuLi-NOAA · 2023-11-13T22:57:21Z

I've found the reference to gravity forumla 1980. If sinlat is the argument passed in to sfc_nst.f, it seems clear that the grv function should be:
real(kind_phys) function grv(x)
  real(kind=kind_phys) :: x  ! sinlat
  real(kind=kind_phys) :: gamma,c1,c2,c3,c4
  gamma=9.7803267715
  c1=0.0052790414
  c2=0.0000232718
  c3=0.0000001262
  c4=0.0000000007

  grv=gamma*(one+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
  !print *,'grav=',grv,lat
end function grv

Denise,

Yes, there is a bug here, as reminded by Aron Wang some time ago.
It is in the calculation of the latitude dependent gravity accelaration. I haven't got a chance to deal with it and the impact should be very small (the variation of the gravity accelaration with latitude).

The best way to handle it is to modify the code in sfc_nst.f, as follows:

add "alat" to the declaration:

real(kind=kind_phys) t12,alon,alat,tsea,sstc,dta,dtz

and

then

change

grav = grv(sinlat(i))

to be

      alat      = asin(sinlat)*rad2deg
      grav      = grv(alat)

Denise,

You are revivew it really carefully.

Can you include the above changes into your PR?

Thanks!

Xu

DeniseWorthen · 2023-11-13T23:14:40Z

@XuLi-NOAA Thanks for confirming. Yes, it is small error at low latitude, but at 45N, the error is about 0.26 % (9.806 vs 9.7803). I actually propose this fix, which eliminates the need to convert from sinlat->lat in sfc_nst.f

diff --git a/physics/module_nst_water_prop.f90 b/physics/module_nst_water_prop.f90
index 5d71ce82..34037c85 100644
--- a/physics/module_nst_water_prop.f90
+++ b/physics/module_nst_water_prop.f90
@@ -487,9 +487,9 @@ contains
   !

 !>\ingroup gfs_nst_main_mod
-function grv(lat)
-  real(kind=kind_phys) :: lat
-  real(kind=kind_phys) :: gamma,c1,c2,c3,c4,pi,phi,x
+real(kind_phys) function grv(x)
+  real(kind=kind_phys) :: x  ! sinlat
+  real(kind=kind_phys) :: gamma,c1,c2,c3,c4
   gamma=9.7803267715
   c1=0.0052790414
   c2=0.0000232718
@@ -497,8 +497,8 @@ function grv(lat)
   c4=0.0000000007
   pi=3.141593

-  phi=lat*pi/180
-  x=sin(phi)
+  !phi=lat*pi/180.0_kp
+  !x=sin(phi)
   grv=gamma*(one+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
   !print *,'grav=',grv,lat
 end function grv

XuLi-NOAA · 2023-11-14T05:19:00Z

@XuLi-NOAA Thanks for confirming. Yes, it is small error at low latitude, but at 45N, the error is about 0.26 % (9.806 vs 9.7803). I actually propose this fix, which eliminates the need to convert from sinlat->lat in sfc_nst.f

diff --git a/physics/module_nst_water_prop.f90 b/physics/module_nst_water_prop.f90
index 5d71ce82..34037c85 100644
--- a/physics/module_nst_water_prop.f90
+++ b/physics/module_nst_water_prop.f90
@@ -487,9 +487,9 @@ contains
   !

 !>\ingroup gfs_nst_main_mod
-function grv(lat)
-  real(kind=kind_phys) :: lat
-  real(kind=kind_phys) :: gamma,c1,c2,c3,c4,pi,phi,x
+real(kind_phys) function grv(x)
+  real(kind=kind_phys) :: x  ! sinlat
+  real(kind=kind_phys) :: gamma,c1,c2,c3,c4
   gamma=9.7803267715
   c1=0.0052790414
   c2=0.0000232718
@@ -497,8 +497,8 @@ function grv(lat)
   c4=0.0000000007
   pi=3.141593

-  phi=lat*pi/180
-  x=sin(phi)
+  !phi=lat*pi/180.0_kp
+  !x=sin(phi)
   grv=gamma*(one+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
   !print *,'grav=',grv,lat
 end function grv

Yes, what you proposed means no need to convert sinlat to lat in sfc_nst.f. However, this conversion occurs at multi-places in the code, since only sinlat but not lat is available (not sure why). The code consistency in all the subroutiones needs to be considered as well.

DeniseWorthen · 2023-11-14T13:13:11Z

@XuLi-NOAA The only use of the function grv is within sfc_nst.f so I think it is safe to simply change the grv function to be consistent with the actual argument (which is sinlat).

~/GitHub/fv3atm_emc/ccpp % grep -ir grv .
Binary file ./framework/doc/DevelopersGuide/CCPP_VARIABLES_FV3.pdf matches
./physics/physics/smoke_dust/dust_fengsha_mod.F90:546:    real :: grvsoilm
./physics/physics/smoke_dust/dust_fengsha_mod.F90:554:    grvsoilm = soilMoistureConvertVol2Grav(slc, sand)
./physics/physics/smoke_dust/dust_fengsha_mod.F90:560:    moistureCorrectionFecan = sqrt(1.0 + 1.21 * max(0., grvsoilm - drylimit)**0.68)
./physics/physics/smoke_dust/dust_fengsha_mod.F90:580:    real :: grvsoilm
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-sw-g112-210809.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-lw-g128-210809.nc matches
./physics/physics/sfc_nst.f:164:     &                                 density,rhocoef,compjd,grv       &
./physics/physics/sfc_nst.f:372:          grav      = grv(sinlat(i))
./physics/physics/module_nst_water_prop.f90:12:  public :: rhocoef,density,sw_rad,sw_rad_aw,sw_rad_sum,sw_rad_upper,sw_rad_upper_aw,sw_rad_skin,grv,solar_time_from_julian,compjd, &
./physics/physics/module_nst_water_prop.f90:486:function grv(lat)
./physics/physics/module_nst_water_prop.f90:498:  grv=gamma*(1+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
./physics/physics/module_nst_water_prop.f90:499:  !print *,'grav=',grv,lat
./physics/physics/module_nst_water_prop.f90:500:end function grv

XuLi-NOAA · 2023-11-14T13:34:01Z

@XuLi-NOAA The only use of the function grv is within sfc_nst.f so I think it is safe to simply change the grv function to be consistent with the actual argument (which is sinlat).

~/GitHub/fv3atm_emc/ccpp % grep -ir grv .
Binary file ./framework/doc/DevelopersGuide/CCPP_VARIABLES_FV3.pdf matches
./physics/physics/smoke_dust/dust_fengsha_mod.F90:546:    real :: grvsoilm
./physics/physics/smoke_dust/dust_fengsha_mod.F90:554:    grvsoilm = soilMoistureConvertVol2Grav(slc, sand)
./physics/physics/smoke_dust/dust_fengsha_mod.F90:560:    moistureCorrectionFecan = sqrt(1.0 + 1.21 * max(0., grvsoilm - drylimit)**0.68)
./physics/physics/smoke_dust/dust_fengsha_mod.F90:580:    real :: grvsoilm
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-sw-g112-210809.nc matches
Binary file ./physics/physics/rte-rrtmgp/rrtmgp/data/rrtmgp-data-lw-g128-210809.nc matches
./physics/physics/sfc_nst.f:164:     &                                 density,rhocoef,compjd,grv       &
./physics/physics/sfc_nst.f:372:          grav      = grv(sinlat(i))
./physics/physics/module_nst_water_prop.f90:12:  public :: rhocoef,density,sw_rad,sw_rad_aw,sw_rad_sum,sw_rad_upper,sw_rad_upper_aw,sw_rad_skin,grv,solar_time_from_julian,compjd, &
./physics/physics/module_nst_water_prop.f90:486:function grv(lat)
./physics/physics/module_nst_water_prop.f90:498:  grv=gamma*(1+(c1*x**2)+(c2*x**4)+(c3*x**6)+(c4*x**8))
./physics/physics/module_nst_water_prop.f90:499:  !print *,'grav=',grv,lat
./physics/physics/module_nst_water_prop.f90:500:end function grv

True, grv is only used in sfc_nst.f. I meant sinlat and asin(sinlat).

DeniseWorthen · 2023-11-14T14:28:47Z

@dustinswales Most tests fail comparison if I fix the grv function (only 78 of 226 pass). I assume CMs would prefer that change to go in as a separate PR, correct?

Also, would it be OK if I cleaned up the old "column 6" continuations in the various nst modules for this PR?

grantfirl · 2023-11-15T19:46:30Z

physics/module_nst_parameters.f90

-       ,rvrdm1 => con_fvirt &
-       ,rd => con_rd &
-       ,rocp => con_rocp  &        !< r/cp          
+  use physcons, only:       &


Note to self: this file only has whitespace changes

grantfirl · 2023-11-15T19:47:16Z

physics/physcons.F90

@@ -33,7 +33,7 @@

 !> This module contains some of the most frequently used math and physics


Note to self: whitespace only

grantfirl

Looks fine to me. The whitespace changes improve the formatting IMO, so it's OK.

DeniseWorthen · 2023-11-15T19:58:34Z

@grantfirl Thanks. What do you think about the issue of the column 6 continuations? I'd like to clean that up too, but that is a lot of white-space also. And I still have a fix to make for Dustin. I'm always curious why the values of zero and one are defined over and over. I thought I would instead add them to the module_nst_parameters.f90 for use by the various nst modules.

dustinswales · 2023-11-16T21:53:00Z

@DeniseWorthen
The change in answers are not surprising, and since we know that we are fixing something that's obviously wrong, I think it's fine to include this change. Maybe others have different opinions.
If you'd like, I can help cleaning up the percision issue I brought up?

DeniseWorthen · 2023-11-16T22:19:24Z

Thanks @dustinswales. I think in general for UWM, my preference is to always keep baseline changing PRs as distinct as possible for ease of going back and isolating where something might have been introduced. We wouldn't want to revert the change for the grv function of course, but it does make ensure the type-variable changes do not impact answers. We could always combine this PR w/ any other that does not change baselines.

For the precision, I think I'd just like to remove the _kp and use _kind_phys where needed, I think that is the fewest changes. There doesn't seem to be any preferred standard...in fact just grepping turned up this one:

tridi.f:14: integer, parameter :: one = 1.0_kind_phys

grantfirl · 2023-11-17T15:36:13Z

@grantfirl Thanks. What do you think about the issue of the column 6 continuations? I'd like to clean that up too, but that is a lot of white-space also. And I still have a fix to make for Dustin. I'm always curious why the values of zero and one are defined over and over. I thought I would instead add them to the module_nst_parameters.f90 for use by the various nst modules.

I would applaud efforts to remove line 6 continuations in f/F90 files. As for redefining zero and one, this is done throughout physics and could therefore be put in the host-defined physical parameters (or somewhere else?) and passed in everywhere. I don't know if it's worth the effort to clean all of that up or not. If you feel strongly about doing it for just NSST, then go ahead, I suppose.

DeniseWorthen · 2023-11-17T15:42:23Z

@grantfirl Thanks. I'll limit myself to the NST modules for any changes.

* remove continuation lines in column 6 and rename files as f90 * remove unused variables and make only needed variables, routines or procedures public * move use statements to module level, remove implicit none and use statements in SRs

DeniseWorthen · 2023-11-27T17:06:35Z

@grantfirl I've pushed changes which remove the continuation lines, and I've tested it again and it is B4B. But I changed the names from ".f" to ".f90" and git is not showing the files as renamed. Instead, the .f files show deleted and the .f90 files show added. I'm not sure how to have folks review....I could perhaps temporarily commit the new .f90 files as .f files?

grantfirl · 2023-11-27T17:21:13Z

@grantfirl I've pushed changes which remove the continuation lines, and I've tested it again and it is B4B. But I changed the names from ".f" to ".f90" and git is not showing the files as renamed. Instead, the .f files show deleted and the .f90 files show added. I'm not sure how to have folks review....I could perhaps temporarily commit the new .f90 files as .f files?

Since the changes are B4B, and it has already been reviewed, I'll just go ahead and review commit f8e1601 using a code difference software on my local machine. It's not a big deal, I don't think.

grantfirl

@mkavulich @Qingfu-Liu I reviewed changes from commit f8e1601 by checking out the previous commit from @DeniseWorthen 's branch, saving the old sft_nst_*.f files, checking out the latest commit, and then inspecting changes via the Meld software. The changes can be classified as: 1) moving module use statements to the module level rather than the subroutine level, 2) using zero, one, etc defined in module_nst_parameters rather than locally, 3) removing the column 6 line continutation (mostly in the declaration section, but some in the main subroutine body for long math operations or function calls), 4) cleaning up whitespace so that it is much more consistent for if blocks, loops, etc.

You're welcome to do the same to review the changes, but since it is B4B, it's probably not worth your time to replicate.

zach1221 · 2023-11-29T21:24:28Z

@grantfirl ufs-wm PR #1993 has finished with the regression testing. Are you able to merge this PR for us, please?

DeniseWorthen and others added 3 commits November 6, 2023 10:49

add kind_phys where missing

9d59cca

Merge branch 'ufs-community:ufs/dev' into feature/fix_nstkind

737fced

clean up type mix-matches

6397387

* add one,zero and half * fix instances of reals compared to integer and integers used in real expressions

DeniseWorthen requested review from XuLi-NOAA, grantfirl, Qingfu-Liu and dustinswales as code owners November 12, 2023 16:04

This was referenced Nov 12, 2023

fix type mis-matches in nst water property module ufs-community/ufs-weather-model#1993

Merged

fix type mis-matches in nst water property module NOAA-EMC/fv3atm#719

Merged

add parameter zero and clean up nst_parameters

3714e76

* fix type mis-match in call to int_epn using parameter zero in module_nst_model

Qingfu-Liu approved these changes Nov 13, 2023

View reviewed changes

XuLi-NOAA approved these changes Nov 13, 2023

View reviewed changes

dustinswales reviewed Nov 13, 2023

View reviewed changes

grantfirl reviewed Nov 15, 2023

View reviewed changes

grantfirl approved these changes Nov 15, 2023

View reviewed changes

DeniseWorthen added 2 commits November 25, 2023 10:41

Merge branch 'ufs/dev' into feature/fix_nstkind

b24beb1

rename and update files

f8e1601

* remove continuation lines in column 6 and rename files as f90 * remove unused variables and make only needed variables, routines or procedures public * move use statements to module level, remove implicit none and use statements in SRs

grantfirl approved these changes Nov 27, 2023

View reviewed changes

grantfirl merged commit a77ed16 into ufs-community:ufs/dev Nov 29, 2023
2 of 3 checks passed

grantfirl mentioned this pull request Feb 23, 2024

UFS-dev PR#126 NCAR/ccpp-physics#1054

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix missing or mis-matched variable types in nst_water_prop #126

fix missing or mis-matched variable types in nst_water_prop #126

DeniseWorthen commented Nov 12, 2023 •

edited

Loading

XuLi-NOAA left a comment

DeniseWorthen commented Nov 13, 2023

DeniseWorthen commented Nov 13, 2023 •

edited

Loading

dustinswales Nov 13, 2023

DeniseWorthen Nov 13, 2023

dustinswales commented Nov 13, 2023 •

edited

Loading

XuLi-NOAA commented Nov 13, 2023

DeniseWorthen commented Nov 13, 2023

XuLi-NOAA commented Nov 14, 2023

DeniseWorthen commented Nov 14, 2023

XuLi-NOAA commented Nov 14, 2023

DeniseWorthen commented Nov 14, 2023 •

edited

Loading

grantfirl Nov 15, 2023

grantfirl Nov 15, 2023

grantfirl left a comment

DeniseWorthen commented Nov 15, 2023

dustinswales commented Nov 16, 2023

DeniseWorthen commented Nov 16, 2023

grantfirl commented Nov 17, 2023

DeniseWorthen commented Nov 17, 2023

DeniseWorthen commented Nov 27, 2023

grantfirl commented Nov 27, 2023

grantfirl left a comment

zach1221 commented Nov 29, 2023

		@@ -33,7 +33,7 @@

		!> This module contains some of the most frequently used math and physics

fix missing or mis-matched variable types in nst_water_prop #126

fix missing or mis-matched variable types in nst_water_prop #126

Conversation

DeniseWorthen commented Nov 12, 2023 • edited Loading

XuLi-NOAA left a comment

Choose a reason for hiding this comment

DeniseWorthen commented Nov 13, 2023

DeniseWorthen commented Nov 13, 2023 • edited Loading

dustinswales Nov 13, 2023

Choose a reason for hiding this comment

DeniseWorthen Nov 13, 2023

Choose a reason for hiding this comment

dustinswales commented Nov 13, 2023 • edited Loading

XuLi-NOAA commented Nov 13, 2023

DeniseWorthen commented Nov 13, 2023

XuLi-NOAA commented Nov 14, 2023

DeniseWorthen commented Nov 14, 2023

XuLi-NOAA commented Nov 14, 2023

DeniseWorthen commented Nov 14, 2023 • edited Loading

grantfirl Nov 15, 2023

Choose a reason for hiding this comment

grantfirl Nov 15, 2023

Choose a reason for hiding this comment

grantfirl left a comment

Choose a reason for hiding this comment

DeniseWorthen commented Nov 15, 2023

dustinswales commented Nov 16, 2023

DeniseWorthen commented Nov 16, 2023

grantfirl commented Nov 17, 2023

DeniseWorthen commented Nov 17, 2023

DeniseWorthen commented Nov 27, 2023

grantfirl commented Nov 27, 2023

grantfirl left a comment

Choose a reason for hiding this comment

zach1221 commented Nov 29, 2023

DeniseWorthen commented Nov 12, 2023 •

edited

Loading

DeniseWorthen commented Nov 13, 2023 •

edited

Loading

dustinswales commented Nov 13, 2023 •

edited

Loading

DeniseWorthen commented Nov 14, 2023 •

edited

Loading