3DTKE EDMF GFS PBL related updates

[BinLiu-NOAA]

Description

Scale-aware 3DTKE EDMF GFS PBL related updates/modifications from @JongilHan66, @zhup01, @samuelkyfung.

• Update and fix the index from ix to im for the sa3dtke related variables (deform_[1-3]) in driver/UFS/atmosphere.F90. (from @JongilHan66, @zhup01, @samuelkyfung)
• Add deallocate for the sa3dtke related variables in model/fv_arrays.F90. (from @zhup01, @samuelkyfung)
• Additional sa3dtke related modifications from @JongilHan66 in dyn_core.F90 regarding computations of def_1(deform_1), def_2(deform_2), def_3(defrom_3); zl: model layers where ua, va, w are located; zm: model interface, zm(k)=0.5*(zl(k)+zl(k+1))
  • def_1(square of vertical shear) is now computed at zm, consistent with shr2 in satmedmfvdifq.F
  • def_2(square_of_horizontal_shear) is computed at zl where tke is also located. Horizontal gradients are computed using (i+1)-(i-1) & (j+1)-(j-1) rather than the current (i+1)-(i) & (j+1)-(j) to include gradients in west and south sides of 'ua' points.
  • def_3(rate of horizontal TKE transfer and pressure correlation) is computed at zl with correction of horizontal gradients of tke(e) & dku3d_e(Kq), e.g., d(Kq(de/dx))/dx is computed as
    d(Kq(de/dx))/dx=[0.5(Kq(i+1)+Kq(i))(e(i+1)-e(i))/dx-0.5(Kq(i-1)+Kq(i))(e(i)-e(i-1))/dx]/dx.
    dku3d_e(kq) & dku3d_h are now computed at zl in satmedmfvdifq.F

Fixes # (issue)

• addresses Update 3DTKE EDMF GFS PBL and combine with TTE-EDMF ufs-community/ufs-weather-model#2815

Related pull requests:

• Scale-aware 3DTKE EDMF and TTE-EDMF GFS PBL related updates ufs-community/ufs-weather-model#2806
• Scale-aware 3DTKE EDMF and TTE-EDMF GFS PBL related modifications NOAA-EMC/ufsatm#993
• Update scale-aware 3DTKE EDMF and introduce TTE-EDMF in GFS PBL scheme ufs-community/ccpp-physics#303

How Has This Been Tested?

Tested under ufs-weather-model and with HAFS and GFS applications.

Checklist:

Please check all whether they apply or not

• My code follows the style guidelines of this project
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• Any dependent changes have been merged and published in downstream modules

[gaokun227]

The shear production terms should be colocated with TKE, which is defined at the layer centers. Therefore, the changes related to du/dz, dv/dz and def_1, which move these variables to the layer interfaces, are unlikely to be meaningful upgrades, as these terms will need to be converted to layer-center values later anyway. In fact, the original method (before this PR) for calculating du/dz and dv/dz is more accurate.

[BinLiu-NOAA]

The shear production terms should be colocated with TKE, which is defined at the layer centers. Therefore, the changes related to du/dz, dv/dz and def_1, which move these variables to the layer interfaces, are unlikely to be meaningful upgrades, as these terms will need to be converted to layer-center values later anyway. In fact, the original method (before this PR) for calculating du/dz and dv/dz is more accurate.

Thanks, @gaokun227! @JongilHan66 for your information/consideration with Kun's comment/suggestion.

[JongilHan66]

The shear production terms should be colocated with TKE, which is defined at the layer centers. Therefore, the changes related to du/dz, dv/dz and def_1, which move these variables to the layer interfaces, are unlikely to be meaningful upgrades, as these terms will need to be converted to layer-center values later anyway. In fact, the original method (before this PR) for calculating du/dz and dv/dz is more accurate.

@gaokun227 Since u & v are defined at the layer center (same as TKE), the computation of du/dz & dv/dz at the layer interfaces would be more accurate. To first interpolate u & v into the layer interfaces and then compute du/dz & dv/dz at the model layer would weaken the actual vertical shear. Also, computation of def_1 at the layer interfaces is consistent with def_1 computation in the 1D TKE scheme.

[lharris4]

The original code used a high-order interpolation of u and v from layer-means (not layer centers) to layer interfaces to trivially compute dudz and dvdz on layers, similar to how high-order reconstructions are done in the Piecewise-Parabolic Method and other high-order finite-volume schemes. This is consistent with how FV3 computes w for vertical straining terms and how it computes advective winds for gz on interfaces in the nonhydrostatic solver.

The new code computing dudz on interfaces assumes that u and v represent layer-centroid values and then interpolates the interface dudz using a simple two-point average, which is low accuracy.

[lharris4]

Please check closely where dwdx, etc. are computed and whether a full-grid-length centered difference is needed here. This will compute dwdx at cell centers instead of cell interfaces, and necessitates the additional interpolation of dwdx below on lines 2338--2357.

[JongilHan66]

Please check closely where dwdx, etc. are computed and whether a full-grid-length centered difference is needed here. This will compute dwdx at cell centers instead of cell interfaces, and necessitates the additional interpolation of dwdx below on lines 2338--2357.

@lharris4 yes, dwdx is computed at cell centers (vertically at layers) and then 2 layer vertical average of dwdx is used for deform_1 computation which is defined at the interfaces

[lharris4]

Why are the branches for the top and bottom layers needed?

[JongilHan66]

Why are the branches for the top and bottom layers needed?

@lharris4 Actually, the deform_1 value at k=1 (npz in the CCPP 3D TKE scheme) is not used in the 3D TKE scheme, but we still assign a reasonable deform_1 value at k=1 for a safety.

[JongilHan66]

The original code used a high-order interpolation of u and v from layer-means (not layer centers) to layer interfaces to trivially compute dudz and dvdz on layers, similar to how high-order reconstructions are done in the Piecewise-Parabolic Method and other high-order finite-volume schemes. This is consistent with how FV3 computes w for vertical straining terms and how it computes advective winds for gz on interfaces in the nonhydrostatic solver.

The new code computing dudz on interfaces assumes that u and v represent layer-centroid values and then interpolates the interface dudz using a simple two-point average, which is low accuracy.

@lharris4 In the new code, we don't interpolate dudz which is defined on interfaces. The deform_1~dudz^2 is also defined on interfaces and transferred into CCPP 3DTKE PBL scheme and multiplied by Km (eddy diffusivity also defined on interfaces). Then, taking average of Kmdeform_1 at 2 interface levels (i.e., 0.5(Km(k)*deform_1(k)+km(k-1)*deform_1(k-1))) is used for the TKE shear production term as TKE is defined on layers. This is consistent with what we did in the 1DTKE PBL scheme.

[gaokun227]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.

However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.

We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

[JongilHan66]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.

However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.

We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

[gaokun227]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

[JongilHan66]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

@gaokun227 How much order the damp_t=damp_vt associated damping is? Is it higher than 2nd-order?

[gaokun227]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

@gaokun227 How much order the damp_t=damp_vt associated damping is? Is it higher than 2nd-order?

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

[JongilHan66]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

@gaokun227 How much order the damp_t=damp_vt associated damping is? Is it higher than 2nd-order?

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

@gaokun227 Thanks!! I find that nord=2 is set in the current GFS workflow. This means that the current GFS is using the 2nd order damping. Right?

[gaokun227]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

@gaokun227 How much order the damp_t=damp_vt associated damping is? Is it higher than 2nd-order?

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

@gaokun227 Thanks!! I find that nord=2 is set in the current GFS workflow. This means that the current GFS is using the 2nd order damping. Right?

@JongilHan66 nord=2 means 4th order. That means GFS and HAFS use the same option.

[JongilHan66]

Just leaving a general note: the updates in this PR look fine, assuming you've tested them and confirmed they produce reasonable results.
However, as we stated in the previous PR, the algorithms used here for the strain-rate tensor element calculations resemble the finite-difference method, and are not entirely consistent with the finite volume method used in FV3.
We have been making upgrades to the 3D TKE scheme at GFDL and incorporating the 3D TKE scheme entirely into the FV3. The upgrades include using finite-volume methods for strain-rate tensor element calculations and introducing the TKE-based diffusion as 2nd-order diffusive fluxes to nearly all variables (including TKE itself). I am leading a documentation paper now and should be able to release the code soon.

@gaokun227 If I understand correctly, there are a vorticity damping (damping coefficient, damp_t=damp_vt) for temperature (T), moisture (q), and cloud condensate (I am not sure whether the other tracers including TKE have this damping or not). Does your TKE-based 2nd-order diffusion imply that the damping coefficient (i.e., damp_t) is replaced by a horizontal turbulent eddy diffusivity (i.e., dku3d_e in the code) or that you add new turbulent eddy diffusion such as a form of d (dku3d_e dT/dx) /dx for all the T, q, and tracers in addition to the damping?

@JongilHan66 Lucas introduced a new 2nd-order diffusive flux component, which is independent from the higher-order diffusive flux you mentioned. We use TKE to compute the diffusion coefficient for the 2nd-order diffusive flux component and the high-order damping code remain untouched (we could turn it off via namelist parameters).

@gaokun227 How much order the damp_t=damp_vt associated damping is? Is it higher than 2nd-order?

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

@gaokun227 Thanks!! I find that nord=2 is set in the current GFS workflow. This means that the current GFS is using the 2nd order damping. Right?

@JongilHan66 nord=2 means 4th order. That means GFS and HAFS use the same option.

@gaokun227 Thanks for clarification!!

[lharris4]

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

@gaokun227 Thanks!! I find that nord=2 is set in the current GFS workflow. This means that the current GFS is using the 2nd order damping. Right?

@JongilHan66 nord=2 means 4th order. That means GFS and HAFS use the same option.

@gaokun227 Thanks for clarification!!

One correction: nord=2 means 6th order for both divergence and vorticity damping (more appropriately, damping applied to dynamical fluxes, not including tracers). If you choose nord=3, divergence damping will be 8th order; due to halo size limitations the vorticity/flux damping remains 6th order.

[gaokun227]

It depends on the choice of nord. We can set it to 2nd order but the most common options are 4th order (e.g., HAFS) and 6th order (GFDL SHiELD suite)

@gaokun227 Thanks!! I find that nord=2 is set in the current GFS workflow. This means that the current GFS is using the 2nd order damping. Right?

@JongilHan66 nord=2 means 4th order. That means GFS and HAFS use the same option.

@gaokun227 Thanks for clarification!!

One correction: nord=2 means 6th order for both divergence and vorticity damping (more appropriately, damping applied to dynamical fluxes, not including tracers). If you choose nord=3, divergence damping will be 8th order; due to halo size limitations the vorticity/flux damping remains 6th order.

Sorry for my misinformation. I forgot about the +1 rule.

[XiaqiongZhou-NOAA]

Reviewed. Looks great to me.

[gaokun227]

I was under the impression that once the 3D TKE scheme is introduced, the TKE-based coefficient should replace the default second-order divergence damping coefficient.

@zhup01 please confirm.

[BinLiu-NOAA]

I was under the impression that once the 3D TKE scheme is introduced, the TKE-based coefficient should replace the default second-order divergence damping coefficient.

@zhup01 please confirm.

@gaokun227, This is a new update/change from @JongilHan66. With this change (damp2 = [Smagorinsky damping] + dt*Kh ), his test from GFS side provides improved skills. Meanwhile, we will test this change from HAFS side as well.

@JongilHan66 and @zhup01, please feel free to provide additional comments/insights. Thanks!

P.S., Another alternative is damp2 = max( [Smagorinsky damping], dt*Kh), for which we will also test from HAFS side.

[zhup01]

Yes, once the 3D TKE scheme is activated, the TKE-based coefficient should replace the default damping coefficient for second-order divergence applied to momentum.

[JongilHan66]

I was under the impression that once the 3D TKE scheme is introduced, the TKE-based coefficient should replace the default second-order divergence damping coefficient.
@zhup01 please confirm.

@gaokun227, This is a new update/change from @JongilHan66. With this change (damp2 = [Smagorinsky damping] + dt*Kh ), his test from GFS side provides improved skills. Meanwhile, we will test this change from HAFS side as well.

@JongilHan66 and @zhup01, please feel free to provide additional comments/insights. Thanks!

P.S., Another alternative is damp2 = max( [Smagorinsky damping], dt*Kh), for which we will also test from HAFS side.

@gaokun227 Changes of damp2 to either damp2 = [Smagorinsky damping] + dtKh or damp2 = max( [Smagorinsky damping], dtKh) (one of which may be determined by HAFS test results) are mainly because the current replacement of Smagorinsky damping with dt*Kh causes the GFS model run to occasionally fail. With the changes, there are no GFS run failures.

[lharris4]

@JongilHan66 Ideally, the best solution would be damp2 = max( [Smagorinsky damping], dtKh) as this would allow us to specify a minimum linear background diffusion to ensure numerical stability, without adding artificial viscosity to the more physical diagnosed damping.

In either case, we could then reduce the background linear diffusion to be as little as we need. The best option would have to be determined through testing.

[JongilHan66]

@JongilHan66 Ideally, the best solution would be damp2 = max( [Smagorinsky damping], dtKh) as this would allow us to specify a minimum linear background diffusion to ensure numerical stability, without adding artificial viscosity to the more physical diagnosed damping.

In either case, we could then reduce the background linear diffusion to be as little as we need. The best option would have to be determined through testing.

@lharris4 Yes, damp2 = max( [Smagorinsky damping], dtKh) looks physically more sound to me than the current damp2 = [Smagorinsky damping]+dtKh. But the latter shows slightly better forecast skill in the GFS tests. I think HAFS team is already testing both options to see if which one gives better HAFS track and intensity forecasts. By the way, Weiguo told me HWRF, HMON, and NMM are also using Smagorinsky type horizontal damping which is proportional to '[deformation+TKE]' (not 'max[deformation, TKE]').

[BinLiu-NOAA]

@lharris4, @gaokun227, @zhup01, An update on this 3DTKE option related damp2 choice in sw_core.F90: We conducted two experiments from HAFS side for some selected 2023 and 2024 NATL storms, comparing using damp2=max([Smagorinsky damping], dtKh) vs damp2=[Smagorinsky damping]+dtKh. Results show that using damp2=max([Smagorinsky damping], dtKh) produces overall slightly better intensity (vmax skills) for forecast lead times of days 4 and 5, along with overall better intensity (vmax) biases, than using damp2=[Smagorinsky damping]+dtKh.

After communicating/discussing with @JongilHan66, he also agrees and recommends using the choice of damp2=max([Smagorinsky damping], dt*Kh).

With that, the related source code has been updated in this pull request branch. Thanks!

[vithikashah001]

@gaokun227 @lharris4 @JongilHan66 this needs to be reviewed

[lharris4]

If the purpose of the modifications to the computations of the deformation terms and the vertical derivatives are to ensure consistency with the 1D PBL schemes, I would like this to be clearly marked. Please note that our forthcoming 3D in-line TKE scheme may not use the same computations.

[vithikashah001]

@gaokun227 @JongilHan66 could you please review this?

[jkbk2004]

all tests are done ok at ufs-community/ufs-weather-model#2806. @laurenchilutti @vithikashah001 can you merge this pr?