Minor modification and bug fix in GFS cumulus convection schemes #232

rhaesung · 2024-10-31T16:58:07Z

The code of this PR was provided by @JongilHan66.

Modified prognostic updraft fraction (sigmab) calculation in 'progsigma_calc.f90' which is physically more sound:
a) moisture convergence calculation: integrate from the convection source level rather than from the cloud base
b) 2D advection of sigmab: sigmab advection averaged over the cloud layers rather than taking a maximum sigmab advection from k=2 to the model top
c) To suppress unrealistically large reflectivity in the model first time step, minimum sigmab at the first time step is set to zero
Fix in missing vertical transport of turbulent kinetic energy (TKE) when aerosol transport is turned on (samfdeepcnv.f & samfshalcnv.f)
Introduce TKE at model layer interfaces (tkeh) for use in convection schemes (GFS_typedefs.F90, GFS_typedefs.meta, satmedmfvdifq.F, satmedmfvdifq.meta, samfdeepcnv.f, samfdeepcnv.meta, samfshalcnv.f, and samfshalcnv.meta)
Vertical transports of hydrometeor variables are currently not allowed in the convection schemes. But vertical transports of number concentrations of only cloud water and ice are mistakenly allowed, which is fixed in this update (CCPP_typedefs.F90)
All the modifications and bug fixes above had neutral impacts on the GFS forecasts

…hen aerosol transport is turned on

JongilHan66 · 2024-10-31T17:04:28Z

@rhaesung Please add one more under item 1. in the PR description:
c) To suppress unrealistically large reflectivity in the model first time step, minimum sigmab at the first time step is set to zero

rhaesung · 2024-10-31T17:06:30Z

@rhaesung Please add one more under item 1. in the PR description: c) To suppress unrealistically large reflectivity in the model first time step, minimum sigmab at the first time step is set to zero

@JongilHan66 Done!

lisa-bengtsson

Looks good!

dustinswales

@rhaesung These changes all look good to me.
I have a few tiny suggestions for the precision to make the code more portable.

physics/CONV/progsigma_calc.f90

grantfirl

Looks good to me, @rhaesung

climbfuji

Looks ok to me.

yangfanglin · 2025-01-06T17:15:14Z

What is the schedule for merging this PR ?

grantfirl · 2025-01-06T17:19:32Z

What is the schedule for merging this PR ?

@yangfanglin I think that the plan was to try to merge the remaining PRs from the RRFS branch and then do this one. If this is a higher priority than the RRFS ones, we can switch up. FYI @jkbk2004

grantfirl · 2025-01-29T22:11:07Z

physics/PBL/SATMEDMF/satmedmfvdifq.meta

+  standard_name = vertical_turbulent_kinetic_energy_at_interface
+  long_name = vertical turbulent kinetic energy at model layer interfaces
+  units = m2 s-2
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)


See https://github.com/NOAA-EMC/fv3atm/pull/885/files#r1934724405. This should be vertical_interface_dimension. Since tkeh is currently being initialized in this scheme up to km1 on line 358, it's actually leaving the top value of the array uninitialized (gets filled with garbage every time this scheme is called). I don't know if this is causing the reproducibility issue because subsequent use seems to be limited to only accessing through km1, but I think it is still a bug worth fixing. Sorry for not catching this in my initial review.

Whatever the vertical dimension of tkeh ends up being, since it is intent(out) it should at least be initialized over its entire domain in this scheme. Uninitialized values always seem to cause problems at some point.

rhaesung · 2025-01-29T23:04:36Z

@grantfirl Thanks for catching this! I'll coordinate with @JongilHan66 to update the scheme based on your feedback.

…sics into conv_hr5update

grantfirl · 2025-02-05T19:56:41Z

physics/CONV/SAMF/samfdeepcnv.f

@@ -1286,6 +1285,22 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
             enddo
           enddo
         enddo
+         kk = ntk -2


What happens when this scheme is called when ntk is initialized to -99 (i.e. when sgs_tke isn't available or being advected)? Are you sure that these code changes only get executed when ntk is > 3? I would think that it would result in a seg fault rather than random behavior if kk ends up being negative though.

I noticed that these changes exist within a 'if (.not. hwrf_samdeep) and if (do_aerosols). Is this consistent with the tests that are non reproducible? If so, I'd think that the culprit is in this change (and the similar one in samfshalcnv.F).

@grantfirl Good finding!!!!!

@rhaesung Please modify them as:

line 1288-1303 in samfdeepcnv.f:

kk = ntk -2 do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem = cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *

& (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
ercko(i,k,kk) = ecko(i,k,kk)
endif
endif
enddo
enddo

=>
if(ntk > 2) then
kk = ntk -2
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kb(i) .and. k < kmax(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem = cq * tem
factor = 1. + tem
ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
& (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
ercko(i,k,kk) = ecko(i,k,kk)
endif
endif
enddo
enddo
endif

line 1095-1110 in samfshalcnv.f:

kk = ntk -2 do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem = cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *

& (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
ercko(i,k,kk) = ecko(i,k,kk)
endif
endif
enddo
enddo

=>

if(ntk > 2) then kk = ntk -2 do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem = cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor ercko(i,k,kk) = ecko(i,k,kk) endif endif enddo enddo endif

Thanks!!

@grantfirl Thanks! @JongilHan66 Sure, I'll test these changes.

@grantfirl @JongilHan66 No luck with this change. Something else is happening. FYI, I tested both with this change and by removing tkeh. The baseline comparison case that failed is: cpld_restart_gfsv17 intel.

rhaesung · 2025-02-10T15:20:11Z

physics/CONV/progsigma_calc.f90

-      if(flag_init .and. .not. flag_restart)then
-         do i = 1,im
-            if(cnvflg(i))then
-               sigmab(i)=0.03


@JongilHan66 It looks like this part is needed to reproduce the baseline for at least the cpld_restart_gfsv17 Intel test. Did you intend to remove this?

@rhaesung Yes, it is intended. sigmab is always initialized as zero in the beginning of progsigma_calc.f90, so I don't think it causes the reproducibility problem. Although it may takes more time, could you test one by one for this PR change from the codes before this PR change?

@JongilHan66 Thank you for the clarification. I understand that sigmab is initialized to zero in progsigma_calc.f90, but given that the cpld_restart_gfsv17 test passed with this part included, it seems like it may still be playing a role in ensuring correct behavior. I will follow your suggestion and test the PR changes step-by-step, comparing with the previous version, to identify what specifically is affecting the reproducibility. I’ll keep you updated with the progress.

rhaesung · 2025-02-10T19:32:30Z

physics/CONV/progsigma_calc.f90

-         enddo
-      endif
+         endif!cnvflg
+      enddo

      do k=1,km
         do i=1,im


@JongilHan66 In line 214, shouldn't sigmab(i) =MAX(sigmind, sigmab(i)) be changed to sigmab(i) = MAX(sigmind_new, sigmab(i))?

@rhaesung You are right. Please change it.

@JongilHan66 does this mean you change the minimum sigma value only at the first time-step?

@rhaesung You are right. Please change it.

@JongilHan66 Okay, this correction has allowed most of the failed test cases to pass, but there are still issues with the reproducibility of the cpld_control_gfsv17 and cpld_restart_gfsv17 tests. I still suspect the changes in progsigma_calc.f90. I'll keep you updated.

You may need: if(flag_init .and. .not. flag_restart)then, because in the case of restart you don't want to use the initial value, but rather the value read in from the restart file.

You may need: if(flag_init .and. .not. flag_restart)then, because in the case of restart you don't want to use the initial value, but rather the value read in from the restart file.

@lisa-bengtsson Thank you for the suggestion! I agree with your point. I’ll incorporate the if(flag_init .and. .not. flag_restart) condition to ensure that we only use the initial value during the initialization phase and not during a restart.

JongilHan66 · 2025-02-11T02:17:06Z

@lisa-bengtsson Yes. Initial value of sigmab=0.03 tended to produce unrealistically large reflectivity at the initial time especially in the RRFS.

lisa-bengtsson · 2025-02-11T04:19:23Z

@lisa-bengtsson Yes. Initial value of sigmab=0.03 tended to produce unrealistically large reflectivity at the initial time especially in the RRFS.

@JongilHan66 I recall now... it's taken some time for this PR to get through the queue.

rhaesung · 2025-02-13T14:57:42Z

@JongilHan66 @grantfirl I confirmed that the changes in progsigma_calc.f90 in PR #252 prevent the completion of the cpld_restart_gfsv17 test when compared to the new baseline, with the following error:

120: FATAL from PE 120: NaN in input field of mpp_reproducing_sum(_2d), this indicates numerical instability

@JongilHan66 Could you please review the logic for sigmind_new to ensure it doesn’t result in invalid values during the simulation? I think we need something like the if(flag_init .and. .not. flag_restart) condition in line 66, since the restart run is failing?

JongilHan66 · 2025-02-13T15:23:39Z

@rhaesung I don't think sigmind_new causes the instability as it is a local constant variable (minimum updraft fraction; sigmind=0.01). Have you run it with the original progsigma_calc.f90 before the modification?

rhaesung · 2025-02-13T15:58:26Z

@rhaesung I don't think sigmind_new causes the instability as it is a local constant variable (minimum updraft fraction; sigmind=0.01). Have you run it with the original progsigma_calc.f90 before the modification?

@JongilHan66 Yes, I recall that the develop branch worked fine. Please review the change in Jili’s PR 252. Since this change prevents the cpld_restart_gfsv17 test from completing with the above error when compared to the new baseline, I suspect the issue is related to the proper handling of sigmind_new and sigmind during the simulation.

JongilHan66 · 2025-02-13T16:13:07Z

@rhaesung The changes associated with sigmind_new & sigmind follow Jili's changes. And I find the current changes are exactly the same as Jili's PR252. In the restart, sigmind_new will be sigmind=0.01. I don't get it why the constant sigind_new causes a problem in restart. How about setting sigmind_new=0.0 before line 66 (i.e., if (flag_init) then)?

JongilHan66 · 2025-02-13T17:05:39Z

@rhaesung @lisa-bengtsson I thought that if flag_restart=.true., flag_init automatically becomes .false, so sigmind_new=sigmind=0.01. Am I wrong?

lisa-bengtsson · 2025-02-13T17:58:15Z

@rhaesung @lisa-bengtsson I thought that if flag_restart=.true., flag_init automatically becomes .false, so sigmind_new=sigmind=0.01. Am I wrong?

I thought so too, but it turned out flag_restart can be true and the flag_init can be true at the same time - indicating that it is the first time-step under restart. That is why I put the .not. flag_restart condition in there.

JongilHan66 · 2025-02-13T18:15:33Z

@lisa-bengtsson Thanks for the clarification!! @rhaesung Please change the 66 line to if(flag_init .and. .not. flag_restart). But I don't think this change can avoid the numerical instability in restart.

rhaesung · 2025-02-13T19:25:35Z

@JongilHan66 @lisa-bengtsson @grantfirl Okay, I confirmed that if(flag_init .and. .not. flag_restart) is at least resolving the issue in PR 252. I’ll re-test the failed tests with this change for this PR.

rhaesung added 4 commits October 25, 2024 18:25

modified prognostic updraft fraction (sigmab) calcultation

14c2f52

fix in missing vertical transport of turbulent kinetic energy (TKE) w…

36e7b63

…hen aerosol transport is turned on

introduce TKE at model layer interfaces (tkeh)

1ace72d

fix compilation error

b325da7

rhaesung requested review from JongilHan66, grantfirl, Qingfu-Liu, dustinswales, lisa-bengtsson and haiqinli as code owners October 31, 2024 16:58

This was referenced Oct 31, 2024

Minor modification and bug fix in GFS cumulus convection schemes NOAA-EMC/fv3atm#885

Open

Minor modification and bug fix in GFS cumulus convection schemes ufs-community/ufs-weather-model#2487

Open

rhaesung requested a review from RuiyuSun October 31, 2024 17:02

lisa-bengtsson approved these changes Oct 31, 2024

View reviewed changes

JongilHan66 approved these changes Oct 31, 2024

View reviewed changes

haiqinli approved these changes Oct 31, 2024

View reviewed changes

dustinswales reviewed Nov 1, 2024

View reviewed changes

physics/CONV/progsigma_calc.f90 Outdated Show resolved Hide resolved

physics/CONV/progsigma_calc.f90 Outdated Show resolved Hide resolved

replace the hardcoded precision constants with the appropriate precision

f6f1892

dustinswales approved these changes Nov 12, 2024

View reviewed changes

grantfirl approved these changes Nov 15, 2024

View reviewed changes

grantfirl requested a review from climbfuji November 15, 2024 16:05

Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update

a1aaf1a

climbfuji approved these changes Nov 21, 2024

View reviewed changes

Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update

39e00e2

Qingfu-Liu approved these changes Dec 9, 2024

View reviewed changes

Fixing the incorrectly modified precision constant

80609a9

rhaesung added 2 commits January 13, 2025 15:50

Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update

c241b99

Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update

2ba9a65

grantfirl reviewed Jan 29, 2025

View reviewed changes

rhaesung added 3 commits February 4, 2025 01:00

Merge remote-tracking branch 'upstream' into conv_hr5update

daeabfa

Merge remote-tracking branch 'upstream' into conv_hr5update

d63e3f1

Merge branch 'conv_hr5update' of https://github.com/rhaesung/ccpp-phy…

485a7c7

…sics into conv_hr5update

grantfirl reviewed Feb 5, 2025

View reviewed changes

rhaesung commented Feb 10, 2025

View reviewed changes

Minor modification and bug fix in GFS cumulus convection schemes #232

Are you sure you want to change the base?

Minor modification and bug fix in GFS cumulus convection schemes #232

Conversation

rhaesung commented Oct 31, 2024 • edited Loading

JongilHan66 commented Oct 31, 2024

rhaesung commented Oct 31, 2024

lisa-bengtsson left a comment

Choose a reason for hiding this comment

dustinswales left a comment

Choose a reason for hiding this comment

grantfirl left a comment

Choose a reason for hiding this comment

climbfuji left a comment

Choose a reason for hiding this comment

yangfanglin commented Jan 6, 2025

grantfirl commented Jan 6, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rhaesung commented Jan 29, 2025

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

JongilHan66 commented Feb 11, 2025

lisa-bengtsson commented Feb 11, 2025

rhaesung commented Feb 13, 2025 • edited Loading

JongilHan66 commented Feb 13, 2025

rhaesung commented Feb 13, 2025

JongilHan66 commented Feb 13, 2025

JongilHan66 commented Feb 13, 2025

lisa-bengtsson commented Feb 13, 2025

JongilHan66 commented Feb 13, 2025

rhaesung commented Feb 13, 2025

rhaesung commented Oct 31, 2024 •

edited

Loading

rhaesung commented Feb 13, 2025 •

edited

Loading