3DTKE EDMF GFS PBL scheme related changes and add an option to use liquid potential temperature in local mixing for GFSPBL PR #288

[BinLiu-NOAA]

Description

Bring in the scale-aware 3DTKE related changes for GFS TKE EDMF PBL scheme developed by Prof. Ping Zhu (@zhup01) from FIU and his group (e.g., @samuelkyfung).

• addresses Bring in the scale-aware 3DTKE EDMF GFS PBL related changes ufs-weather-model#2732

Combining with PR #288: Add an option to use liquid potential temperature in local mixing for GFSPBL

Add a namelist option in GFS TKE EDMF PBL scheme to apply an approximate way to represent local mixing/diffusion process in T-equation using liquid/ice water potential temperature. (from @WeiguoWang-NOAA)

• address Add an option in GFS TKE EDMF PBL scheme to use liquid/ice potential temperature fluxes ufs-weather-model#2762

[grantfirl]

@mkavulich The variables added to this scheme are all new and their standard names should be reviewed. Can you help with this? I will put my suggestions here, but please chime in whether you agree/not.

[grantfirl]

"square_of_vertical_shear_due_to_dynamics"

[grantfirl]

"square_of_horizontal_shear_due_to_dynamics"

[grantfirl]

"transfer_rate_of_tke_due_to_dynamics"

[grantfirl]

do_scale_aware_3d_tke

[grantfirl]

horizontal_atmosphere_momentum_diffusivity_for_dyanmics

[grantfirl]

I'm confused what this variable represents. Can you clarify? I understand that it is a diffusivity (K) used in the 3D TKE scheme. Is it used for the diffusion of TKE itself? In this case, why put "momentum" in there?

[zhup01]

dku3d_e is the horizontal eddy diffusivity for scalar (TKE). The standard_name is a mistake. We will double-check all the standard names and make necessary changes. Thanks.

[samuelkyfung]

Thank you for the comments! We have updated the standard names for the variables def_1, def_2, def_3, dku3d_h, and dku3d_e.

[grantfirl]

may consider saving 1/3 as variable for performance reasons

[samuelkyfung]

Thanks! We have added a variable inv3 as 1/3 for this calculation.

[JongilHan66]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

[zhup01]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

[JongilHan66]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

Eq.7e in Deardorff, Km (eddy diffusivity for momentum) is used with the factor of 2 rather than Ke (eddy diffusivity for TKE, which is set to be same as Kh [eddy diffusivity for heat and moisture] and is larger than Km in unstable condition; i.e., Ke=Kh=Km/prnum, prnum is Prantl number and is smaller than 1 in unstable condition). In all the text books (e.g., Eq.13.37 in Arya,, 1988) and Journal papers I have read, TKE transfer and pressure correlation was parameterized as -Ke d(TKE)/dz. Here, note that Ke is used, NOT Km.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

I've never seen 3D shear is used for the gradient Richardson number (Ri) calculation. Note that bf is vertical buoyancy and does not include horizontal buoyancy. 3D buoyancy may have to be used for Ri calculation if 3D shear should be used. By the way, is there any reference where 3D shear is used if you know?

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

I still argue that def_2 should be removed. Larger 3D TKE would increase the vertical mixing length (ml). Vertical wind shear can tilt vertical parcel movement horizontally and consequently, reduce the vertical mixing length (e.g., see Han et al, 2024, Updates in the NCEP GFS PBL and Convection Models with Environmental Wind Shear Effect and Modified Entrainment and Detrainment Rates and Their Impacts on the GFS Hurricane and CAPE Forecasts. Weather and Forecasting, Vol 39, 679-688). Physically, horizontal shear may reduce the horizontal mixing length.

[zhup01]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

Eq.7e in Deardorff, Km (eddy diffusivity for momentum) is used with the factor of 2 rather than Ke (eddy diffusivity for TKE, which is set to be same as Kh [eddy diffusivity for heat and moisture] and is larger than Km in unstable condition; i.e., Ke=Kh=Km/prnum, prnum is Prantl number and is smaller than 1 in unstable condition). In all the text books (e.g., Eq.13.37 in Arya,, 1988) and Journal papers I have read, TKE transfer and pressure correlation was parameterized as -Ke d(TKE)/dz. Here, note that Ke is used, NOT Km.

Yes, in Deardorff (1980) it used Km, but TKE was not treated as a passive tracer in his paper. In HAFS, however, it is treated as a passive tracer. But how can one justify replacing 2Km with Kq is correct, but not 2kq? To me, doubling eddy diffusivity here is appropriate because it considers both TKE transport and pressure diffusion. After all, this is a parameterization and I don't think there is any observation available to back up which is correct for parameterizing the TKE transport and pressure diffusion. It just people's preference. I don't think eithe Deardorff or Arya had observations to back up their use of 2km or kq.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

I've never seen 3D shear is used for the gradient Richardson number (Ri) calculation. Note that bf is vertical buoyancy and does not include horizontal buoyancy. 3D buoyancy may have to be used for Ri calculation if 3D shear should be used. By the way, is there any reference where 3D shear is used if you know?

By definition, Richardson number is the ratio of buoyancy production to shear production in the TKE budget equation (Stull 1988). For typical boundary layer problems, only 1D TKE budget equation is considered under horizontal homogeneous assumption in which horizontal shear is neglected in the TKE shear production. But Richardson number is a measure of static stability, horizontal shear certainly will affect static stability when it is not negligible. AI answer is that " In summary, while the standard Richardson number calculation prioritizes vertical shear production, horizontal shear might be included in specialized studies or models depending on the specific phenomenon being analyzed. "

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

I still argue that def_2 should be removed. Larger 3D TKE would increase the vertical mixing length (ml). Vertical wind shear can tilt vertical parcel movement horizontally and consequently, reduce the vertical mixing length (e.g., see Han et al, 2024, Updates in the NCEP GFS PBL and Convection Models with Environmental Wind Shear Effect and Modified Entrainment and Detrainment Rates and Their Impacts on the GFS Hurricane and CAPE Forecasts. Weather and Forecasting, Vol 39, 679-688). Physically, horizontal shear may reduce the horizontal mixing length.

Since we blend eddy diffusivities any way, it OK we reduce back to shr2 in the 1D framework. I don't think it will affect much.

[samuelkyfung]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

Thank you for the reminder, we have added sa3dtke as a flag in both ‘mfpblq.f’ and ‘mfscuq.f’. When sa3dtke is true, sigma will be set to 0.

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

The index order are reversed now, the order is: k do-loop then i do-loop.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

Eq.7e in Deardorff, Km (eddy diffusivity for momentum) is used with the factor of 2 rather than Ke (eddy diffusivity for TKE, which is set to be same as Kh [eddy diffusivity for heat and moisture] and is larger than Km in unstable condition; i.e., Ke=Kh=Km/prnum, prnum is Prantl number and is smaller than 1 in unstable condition). In all the text books (e.g., Eq.13.37 in Arya,, 1988) and Journal papers I have read, TKE transfer and pressure correlation was parameterized as -Ke d(TKE)/dz. Here, note that Ke is used, NOT Km.

Yes, in Deardorff (1980) it used Km, but TKE was not treated as a passive tracer in his paper. In HAFS, however, it is treated as a passive tracer. But how can one justify replacing 2Km with Kq is correct, but not 2kq? To me, doubling eddy diffusivity here is appropriate because it considers both TKE transport and pressure diffusion. After all, this is a parameterization and I don't think there is any observation available to back up which is correct for parameterizing the TKE transport and pressure diffusion. It just people's preference. I don't think eithe Deardorff or Arya had observations to back up their use of 2km or kq.

We are still keeping the calculation of def_3 unchanged at this moment.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

I've never seen 3D shear is used for the gradient Richardson number (Ri) calculation. Note that bf is vertical buoyancy and does not include horizontal buoyancy. 3D buoyancy may have to be used for Ri calculation if 3D shear should be used. By the way, is there any reference where 3D shear is used if you know?

By definition, Richardson number is the ratio of buoyancy production to shear production in the TKE budget equation (Stull 1988). For typical boundary layer problems, only 1D TKE budget equation is considered under horizontal homogeneous assumption in which horizontal shear is neglected in the TKE shear production. But Richardson number is a measure of static stability, horizontal shear certainly will affect static stability when it is not negligible. AI answer is that " In summary, while the standard Richardson number calculation prioritizes vertical shear production, horizontal shear might be included in specialized studies or models depending on the specific phenomenon being analyzed. "

We also keep using the 3D shear for the gradient Richardson number (Ri) calculation at this moment.

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

Thank you for the comment. We have made the corresponding changes.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

Thank you for the stylistic comment. We have made the modification for ‘.lt.’ and .gt..

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

I still argue that def_2 should be removed. Larger 3D TKE would increase the vertical mixing length (ml). Vertical wind shear can tilt vertical parcel movement horizontally and consequently, reduce the vertical mixing length (e.g., see Han et al, 2024, Updates in the NCEP GFS PBL and Convection Models with Environmental Wind Shear Effect and Modified Entrainment and Detrainment Rates and Their Impacts on the GFS Hurricane and CAPE Forecasts. Weather and Forecasting, Vol 39, 679-688). Physically, horizontal shear may reduce the horizontal mixing length.

Since we blend eddy diffusivities any way, it OK we reduce back to shr2 in the 1D framework. I don't think it will affect much.

Thank you for point this out! We have switched back to use shr2 for calculating the vertical mixing length in the 1D framework.

[BinLiu-NOAA]

PR #288 now has been merged into this PR #279. Also synced with the latest dev/ufs branch as of 06/20/2025.

[JongilHan66]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

[JongilHan66]

1. Double counting scale-aware (SA) nonlocal mass flux (MF) mixing: the current scheme introduces a partition function, ‘pfnl’, to have SA (i.e., the nonlocal MF approaches zero when the model approaches LES resolution). However, the current nonlocal scheme already has a SA parameterization, i.e., M’=(1-sigma)^2 M, where M is MF, M’ is updated MF with SA, and sigma is the updraft fraction (please see Eqs. 40 & 41 in Han and Bretherton, 2019). When the model approaches LES resolution, sigma is close to 1.0 and the nonlocal MF approaches zero. So, not to double count the SA, I recommend to remove ‘pfnl’. Or, if you prefer your SA parameterization with ‘pfnl’, you need to remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’.

You are right. We missed this part. We will remove ‘sigma’ in the MF subroutines of ‘mfpblq.f’ and ‘mfscuq.f’

2. Index order of i and k in do-loop (L1556-1562; L3087-3094): about 20 years ago, I experienced the reverse order of do-loop (i.e., do i; do k rather than the current do k;do i) caused significantly larger computing time (I guess it was an IBM machine). After that, I’ve never used the reverse do-loop order. Although I am not sure the reverse order may also increase the do-loop computing time for the current machines, I suggest changing the reverse i-k do-loop order to k-i order if you can.

The index order of i and k in do-loop does not appear to be an issue now. We’ve run it on WCOSS2, Orion, Jet, and Hercules. We don’t see any problems. For operational 2 km resolution, there is virtually no difference in computing time for 126-hour forecast between the default and modified HAFS. Even at 800 m resolution, both default and modified HAFS can complete 75-hour simulations in 8-hour walltime. But certainly we can switch the order of i and k.

3. Note that eddy diffusivities, def_1 and def_2 should be defined at model interfaces (zm levels in the code; zl is model layer height [zl(1)~10; zl(2)~30], zi is model interface height [zi(1)=0, zi(2)~20], and zm(k)=zi(k+1) [zm(1)~20]). Since all the prognostic variables except w are located at zl layers, please confirm that def_1 & def_2 are defined at zm levels in dycore.

First of all, def_1 and def_2 are not the eddy diffusivities. They are equivalent to shr2 in 1D framework. The calculation of def_1 and def_2 in dycore uses A-grid winds, which are the same velocities transferring to CCPP. So, they are consistent with the calculation of shr2 in the 1D EDMF PBL scheme in CCPP. They are defined at zm levels.

It seems def_2(square of horizontal shear) is computed in dyn_core.F90 & sw_core.F90. But it looks like def_2 (which involves dudx, dvdy, dudy, dvdx, etc.) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. Am I wrong?

4. L1872: as I understand, def_3 should be read as “rate of ‘horizontal’ TKE transfer and pressure correlation” because the vertical rate is used for implicit tridiagonal solver. Also, it should be computed at TKE layers (i.e., zl layers), unlike def_1 & def_2 (which are computed at zm levels). By the way, in L1872 why the factor 2 is multiplied? Eq.(5) in Han and Bretherton (2019) does not have the factor 2 multiplication.

We are using the formula of Deardorff (1980), i.e., Eq. 7e in his paper. It should have a factor of 2 there. This is because the equation includes both TKE transport and pressure diffusion. If splitting it into two separate components and each follows the down-gradient transport, you will have a factor of 2 there.

5. L1405 & 1612: ri is gradient Richardson number, defined as buoyancy divided by square of vertical wind shear. So ri=bf/def_1, not ri=bf/(def_1+def_2).

No, Richardson number is defined as the ratio of buoyancy production to shear production. In 1D case, only vertical shear is considered. But in 3D case, static stability is affected by 3D shear represented by def_1+def_2 (equivalent to shr2 in 1D framework). Def_1 is only a component of 3D shear production. Please see Eq. 11 in the manuscript that we submitted.

6. L1534-1536: thetal is defined at zl layers, so dz should be zl(i,k+1)-zl(i,k), not zi(i,k+1)-zi(i,k). 1/dz, i.e., 1/(zl(i,k+1)-zl(i,k)) can be replaced by rdzt(i,k) for saving computing time. Also, change grav/theta(i,1) to gotvx(i,k).

I agree. We can make the change as suggested.

7. L1618-1619, L3092: although dkq(diffusivity for TKE)=dkt(diffusivities for other scalar fields) in the current scheme, i.e., dkq=prtkedkt (prtke=1.0), it is better to also distinguish horizontal diffusivities for TKE and other scalar fields. So, I suggest to introduce dkt_h & dku3d_t and define dkt_h=dku_h/prnum, dkq_h=prtke*dkt_h, and dku3d_t=dkt_h

For the current version of the scheme, we may not need to introduce dkt_h to distinguish it from dkq_h because in the dycore the potential temperature equation is solved still using vorticity damping. So, there is no need to define horizontal eddy diffusivity for heat that is not used in the dycore. But in future when we apply the 3D-SA TKE-based damping to all physical quantities, we will make the change accordingly.

8. Change ‘.lt.’=> ‘<’ and ‘.gt.’ => ‘>’ for consistency.

Yes, we will do this.

9. L1199 & L1225: remove ‘def_2’ because it computes vertical wind shear effect on vertical mixing length.

def_1+def_2 in 3D framework is equivalent to shr2 in 1D framework. In the real world, all components of shear tensor will have an effect on vertical mixing. This is similar to Richardson number. In 1D framework, only vertical shear can be considered. But in 3D framework, all components in shear tensor can affect stability. I’d like to keep it because when resolution becomes sufficiently low, the horizontal shear will be negligible. So, it will not affect much at low resolution.

[JongilHan66]

In addition, dku3d_h & dku3d_e are defined at zm levels in satmedmfvdifq.F. But it seems that they are used assuming defined at zl layers in sw_core.F90 & dyn_core.F90. If I am right, they should be defined at zl layers in satmedmfvdifq.F

[zhup01]

In addition, dku3d_h & dku3d_e are defined at zm levels in satmedmfvdifq.F. But it seems that they are used assuming defined at zl layers in sw_core.F90 & dyn_core.F90. If I am right, they should be defined at zl layers in satmedmfvdifq.F

dku3d_h & dku3d_e are parameterized based on tke. So, they are defined at the same level of tke. I am travelling now. I'll look into this when I am back next week.

[JongilHan66]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.

Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

[zhup01]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.

Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

[BinLiu-NOAA]

Thanks, @zhup01! I know you might be on travel now. Wondering if you or @samuelkyfung could make this related changes and commit them back soon (saying today or tomorrow), or you might need additional time to make these changes. Currently, this PR is in the commit queue and almost ready for ufs-weather-model level multiple platforms regression tests. Thanks!

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

[JongilHan66]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2dwdz2+(dudz2+dvdz**2)+(dwdxdudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.

2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdxdudz+dwdydvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.

3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.

4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).

5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

[zhup01]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2_dwdz2+(dudz2+dvdz**2)+(dwdx_dudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.
2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdx_dudz+dwdy_dvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.
3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.
4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).
5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

Jongil, thanks for the comments. dku3d_e and dku3d_h are the eddy diffusivities blended by dku (LSMS extreme) and dku_les (LES extreme), so they should be defined at the same vertical levels as dku (LSMS). This is how they are currently calculated. As for how to appropriately calculate horizontal gradients using A-grid winds is something that we need to further discuss with the FV3 dycore developer at GFDL. Since the grid-box structure is complicated, we need to determine what is the best way to calculate horizontal gradients to be compatible with 1D structure of satmedmfvdifq.F. Kun previously also suggested using D-grid winds for calculating 3D gradients for 3D TKE scheme. The problem becomes complicated. I am in China now. I suggest that we schedule a separate meeting to discuss about it after I am back next week.

[jkbk2004]

@BinLiu-NOAA @zhup01 @grantfirl @JongilHan66 we are trying to schedule ufs-community/ufs-weather-model#2734. If additional patch-up or update is needed, then another quick fix pr can be followed up after @zhup01 comes back. So are we ok to move on for full testing ufs-community/ufs-weather-model#2734 and target to commit by tomorrow?

[yangfanglin]

I'd suggest we move forward with this PR commit, and make future updates in another PR

[JongilHan66]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2_dwdz2+(dudz2+dvdz**2)+(dwdx_dudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.
2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdx_dudz+dwdy_dvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.
3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.
4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).
5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

Jongil, thanks for the comments. dku3d_e and dku3d_h are the eddy diffusivities blended by dku (LSMS extreme) and dku_les (LES extreme), so they should be defined at the same vertical levels as dku (LSMS). This is how they are currently calculated. As for how to appropriately calculate horizontal gradients using A-grid winds is something that we need to further discuss with the FV3 dycore developer at GFDL. Since the grid-box structure is complicated, we need to determine what is the best way to calculate horizontal gradients to be compatible with 1D structure of satmedmfvdifq.F. Kun previously also suggested using D-grid winds for calculating 3D gradients for 3D TKE scheme. The problem becomes complicated. I am in China now. I suggest that we schedule a separate meeting to discuss about it after I am back next week.

As for dku3d_e and dku3d_h, they should be consistent between satmedmfvdifq.F and dycore. Since they are used in zl layers in dycore, they should be re-computed in zl layers in satmedmfvdifq.F (currently they are defined at zm levels). Or, if you don't want to re-compute them in satmedmfvdifq.F, in dycore they should be re-defined as dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) [for k=1, dku3d_h(1)=0.5dku3d_h(1)] and dku3d_e(k)=0.5(dku3d_e(k-1)+dku3d_e(k)).

[BinLiu-NOAA]

I'd suggest we move forward with this PR commit, and make future updates in another PR

Agree with @yangfanglin. Probably we can create a follow-up issue and PR to further improve the related source codes.

@JongilHan66 Wondering if this sounds reasonable to you as well. Thanks!

[zhup01]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2_dwdz2+(dudz2+dvdz**2)+(dwdx_dudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.
2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdx_dudz+dwdy_dvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.
3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.
4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).
5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

Jongil, thanks for the comments. dku3d_e and dku3d_h are the eddy diffusivities blended by dku (LSMS extreme) and dku_les (LES extreme), so they should be defined at the same vertical levels as dku (LSMS). This is how they are currently calculated. As for how to appropriately calculate horizontal gradients using A-grid winds is something that we need to further discuss with the FV3 dycore developer at GFDL. Since the grid-box structure is complicated, we need to determine what is the best way to calculate horizontal gradients to be compatible with 1D structure of satmedmfvdifq.F. Kun previously also suggested using D-grid winds for calculating 3D gradients for 3D TKE scheme. The problem becomes complicated. I am in China now. I suggest that we schedule a separate meeting to discuss about it after I am back next week.

As for dku3d_e and dku3d_h, they should be consistent between satmedmfvdifq.F and dycore. Since they are used in zl layers in dycore, they should be re-computed in zl layers in satmedmfvdifq.F (currently they are defined at zm levels). Or, if you don't want to re-compute them in satmedmfvdifq.F, in dycore they should be re-defined as dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) [for k=1, dku3d_h(1)=0.5_dku3d_h(1)] and dku3d_e(k)=0.5_(dku3d_e(k-1)+dku3d_e(k)).

I am not quite sure about this. Note that A-grid wind in dycore is the same A-grid wind in satmedmfvdifq.F. dku3d_e and dku3d_h are calculated in satmedmfvdifq.F and passed to dycore. So, dku3d_e and dku3d_h are at the same levels in dycore and satmedmfvdifq.F.

dku3d_e is used in dyn_core.F90 for calculating TKE gradient calculation. In satmedmfvdifq.F, dku3d_e and dku3d_h are calculated based on 0.5(tke(i,k)+tke(i,k+1)). To be consistent with this, I've made following changes in dyn_core.F90,
do k=2,npz !pzhu, changed
!-------------------------------------
do j=js,je+2
do i=is,ie+2
tke_1(i,j)=0.5*(q(i,j,k,ntke)+q(i,j,k-1,ntke)) !pzhu, changed.
enddo
enddo
Note that in dycore the vertical level is opposite to CCPP, so level npz in dycore is level 1 in ccpp.
I'll ask Samuel to push this change back if you agree.

dku3d_h is used in sw_core.F90 for damping coefficient calculation. I am not sure if dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) is a simple solution because divergence is calculated using D-grid winds. So, I'd keep this and address the issue in future PR.

[JongilHan66]

I'd suggest we move forward with this PR commit, and make future updates in another PR

Agree with @yangfanglin. Probably we can create a follow-up issue and PR to further improve the related source codes.

@JongilHan66 Wondering if this sounds reasonable to you as well. Thanks!

@BinLiu-NOAA yes, it sounds reasonable to me.

[BinLiu-NOAA]

I'd suggest we move forward with this PR commit, and make future updates in another PR

Agree with @yangfanglin. Probably we can create a follow-up issue and PR to further improve the related source codes.
@JongilHan66 Wondering if this sounds reasonable to you as well. Thanks!

@BinLiu-NOAA yes, it sounds reasonable to me.

Thanks, @JongilHan66!

With that, @jkbk2004, I think this PR and the related ufs-community/ufs-weather-model/pull/2734 are ready for the commit queue now. Thanks

[JongilHan66]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2_dwdz2+(dudz2+dvdz**2)+(dwdx_dudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.
2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdx_dudz+dwdy_dvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.
3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.
4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).
5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

Jongil, thanks for the comments. dku3d_e and dku3d_h are the eddy diffusivities blended by dku (LSMS extreme) and dku_les (LES extreme), so they should be defined at the same vertical levels as dku (LSMS). This is how they are currently calculated. As for how to appropriately calculate horizontal gradients using A-grid winds is something that we need to further discuss with the FV3 dycore developer at GFDL. Since the grid-box structure is complicated, we need to determine what is the best way to calculate horizontal gradients to be compatible with 1D structure of satmedmfvdifq.F. Kun previously also suggested using D-grid winds for calculating 3D gradients for 3D TKE scheme. The problem becomes complicated. I am in China now. I suggest that we schedule a separate meeting to discuss about it after I am back next week.

As for dku3d_e and dku3d_h, they should be consistent between satmedmfvdifq.F and dycore. Since they are used in zl layers in dycore, they should be re-computed in zl layers in satmedmfvdifq.F (currently they are defined at zm levels). Or, if you don't want to re-compute them in satmedmfvdifq.F, in dycore they should be re-defined as dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) [for k=1, dku3d_h(1)=0.5_dku3d_h(1)] and dku3d_e(k)=0.5_(dku3d_e(k-1)+dku3d_e(k)).

I am not quite sure about this. Note that A-grid wind in dycore is the same A-grid wind in satmedmfvdifq.F. dku3d_e and dku3d_h are calculated in satmedmfvdifq.F and passed to dycore. So, dku3d_e and dku3d_h are at the same levels in dycore and satmedmfvdifq.F.

dku3d_e is used in dyn_core.F90 for calculating TKE gradient calculation. In satmedmfvdifq.F, dku3d_e and dku3d_h are calculated based on 0.5(tke(i,k)+tke(i,k+1)). To be consistent with this, I've made following changes in dyn_core.F90, do k=2,npz !pzhu, changed !------------------------------------- do j=js,je+2 do i=is,ie+2 tke_1(i,j)=0.5*(q(i,j,k,ntke)+q(i,j,k-1,ntke)) !pzhu, changed. enddo enddo Note that in dycore the vertical level is opposite to CCPP, so level npz in dycore is level 1 in ccpp. I'll ask Samuel to push this change back if you agree.

dku3d_h is used in sw_core.F90 for damping coefficient calculation. I am not sure if dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) is a simple solution because divergence is calculated using D-grid winds. So, I'd keep this and address the issue in future PR.

In fact, I don't know A or D grid structures. I don't understand why we need to compute tke_1(i,j)=0.5*(q(i,j,k,ntke)+q(i,j,k-1,ntke)). Since def_3 is used in tke layers (i.e., zl layer) in satmedmfvdifq.F, def_3 should also be estimated in zl layers in dycore. As I understand, all the prognostic variables (I am not sure for the vertical velocity) are located in zl layers in dycore like in CCPP physics and so, dku3d_h should be also located in the zl layers in dycore. I think divergence damping or horizontal diffusion may need a horizontal gradient of such as dku3d_h*dudx, so they may need dku3d_h at i+1, i-1, j-1, j+1, not just at i,j, but I am not quite sure.

[zhup01]

In satmedmfvdifq.F, dku3d_h is given as a function of tkeh (tkeh=0.5*(tke(k)+tke(k+1))), which is defined at zm levels.
Also, it seems that def_2 (which is computed based on dudx, dvdy, dudy, dvdx, etc., in dyn_core.F90 & sw_core.F90) is defined at zl layers (where prognostic u & v are located) rather than at zm levels. If it is right, tke production associated with def_2 should be corrected.

def_1 and def_2 are defined between the levels where ua, va, and wa are defined. dku3d_e and dku3d_h are defined at the same levels of dku, which is defined at the level of tkeh. In satmedmfvdifq.F, they are correctly defined . However, in dyn_core.F90, currently tke is used. I'll change it to tkeh=0.5*(tke(k)+tke(k+1)) and make it consistent with satmedmfvdifq.F. We will let you know when it's done.

To be correctly used in satmedmfvdifq.F, def_1, def_2, & def_3 should be defined at the center of horizontal grid box:

1. def_1=2_dwdz2+(dudz2+dvdz**2)+(dwdx_dudz+dwdy*dvdz): as I understand, ua, va, & w are defined at the center of vertical grid box (i.e., zl layer) and ue, ve, & we are defined at the edges of vertical grid box (i.e., model interface (zi) level). So, dudz=(ue(k+1)-ue(k))/dz, dvdz=(ve(k+1)-ve(k))/dz, & dwdz=(we(k+1)-we(k))/dz are defined at zl layer (the same layer where tke is defined). dwdx~(w(i+1)-w(i))/dx & dwdy~(w(j+1)-w(j))/dy are also defined at zl layer. So, def_1 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_1 is defined at zm(k)=zi(k+1) level. On the other hand, it seems that ua, va, & w are defined at the center of horizontal grid box. If so, dwdx & dwdy should be computed proportional to (w(i+1)-w(i-1))/2dx & dwdy~(w(j+1)-w(j-1))/2dy, in order for them to be defined at the center of the center of horizontal grid box where dudz, dvdz, & dwdz are defined.
2. def_2=2*(dudx2+dvdy2)+(dudy+dvdx)2+(dwdx2+dwdy**2)+(dwdx_dudz+dwdy_dvdz):
   dudx~(ua(i+1)-ua(i))/dx, dvdy~(va(j+1)-va(j))/dy, dudy~(ua(j+1)-ua(j))/dy, dvdx~(va(i+1)-va(i))/dx. All of dudx, dvdy, dudy, dvdx, dwdx, dwdy, dudz, dvdz are defined at the center of vertical grid box (zl layer) and thus, def_2 is defined at zl layer, which is inconsistent with satmedmfvdifq.F where def_2 is defined at zm(k)=zi(k+1) level. Again, if ua, va, & w are defined at the center of horizontal grid box, dudx, dvdy, dudy, dvdx should be computed proportional to (ua(i+1)-ua(i-1))/2dx, (va(j+1)-va(j-1))/2dy, dudy~(ua(j+1)-ua(j-1))/2dy, dvdx~(va(i+1)-va(i-1))/2dx, respectively, in order for them to be defined at the center of the center of horizontal grid box.
3. def_3=d(kq dedx)/dx + d(kq dedy)/dy:
   dedx~(e(i+1)-e(i))/dx, dedy~(e(j+1)-e(j))/dy, dkdx~(dku3d_e(i+1)-dku3d_e(i))/dx, dkdy~(dku3d_e(j+1)-dku3d_e(j))/dy. Since e is located at the center of vertical grid box (zl layer), def_3 is defined at zl layer, which is consistent with satmedmfvdifq.F where def_3 is assumed to be defined at zl layer. However, dku3d_e is defined at zm level in satmedmfvdifq.F, which is inconsistent. So, dku3d_e should be computed at the zl layer in satmedmfvdifq.F. Again, since e & dku3d_e are defined at the center of horizontal grid box, dedx, dedy, dkdx, dkdy should be computed proportional to (e(i+1)-e(i-1))/2dx, (e(j+1)-e(j-1))/2dy, (dku3d_e(i+1)-dku3d_e(i-1))/2dx, (dku3d_e(j+1)-dku3d_e(j-1))/2dy, respectively, in order for them to be defined at the center of the center of horizontal grid box.
4. sw_core.F90: to be used in 'damp2' computation, dku3d_h should be computed at the zl layer in satmedmfvdifq.F. Also, since dku3d_h is defined at the center of horizontal grid box, horizontal diffusion associated with damp2 may need dku3d_h information at adjacent grid points (i.e., dku3d_h at i+1, i-1, j+1, j-1).
5. Initial GFS tests indicate that the GFS performance are very sensitive to especially def_1 & def_2. So, estimation of def_1, def_2, and def_3 at correct same points (e.g., center point of cubic volume) would be very important.

Jongil, thanks for the comments. dku3d_e and dku3d_h are the eddy diffusivities blended by dku (LSMS extreme) and dku_les (LES extreme), so they should be defined at the same vertical levels as dku (LSMS). This is how they are currently calculated. As for how to appropriately calculate horizontal gradients using A-grid winds is something that we need to further discuss with the FV3 dycore developer at GFDL. Since the grid-box structure is complicated, we need to determine what is the best way to calculate horizontal gradients to be compatible with 1D structure of satmedmfvdifq.F. Kun previously also suggested using D-grid winds for calculating 3D gradients for 3D TKE scheme. The problem becomes complicated. I am in China now. I suggest that we schedule a separate meeting to discuss about it after I am back next week.

As for dku3d_e and dku3d_h, they should be consistent between satmedmfvdifq.F and dycore. Since they are used in zl layers in dycore, they should be re-computed in zl layers in satmedmfvdifq.F (currently they are defined at zm levels). Or, if you don't want to re-compute them in satmedmfvdifq.F, in dycore they should be re-defined as dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) [for k=1, dku3d_h(1)=0.5_dku3d_h(1)] and dku3d_e(k)=0.5_(dku3d_e(k-1)+dku3d_e(k)).

I am not quite sure about this. Note that A-grid wind in dycore is the same A-grid wind in satmedmfvdifq.F. dku3d_e and dku3d_h are calculated in satmedmfvdifq.F and passed to dycore. So, dku3d_e and dku3d_h are at the same levels in dycore and satmedmfvdifq.F.
dku3d_e is used in dyn_core.F90 for calculating TKE gradient calculation. In satmedmfvdifq.F, dku3d_e and dku3d_h are calculated based on 0.5(tke(i,k)+tke(i,k+1)). To be consistent with this, I've made following changes in dyn_core.F90, do k=2,npz !pzhu, changed !------------------------------------- do j=js,je+2 do i=is,ie+2 tke_1(i,j)=0.5*(q(i,j,k,ntke)+q(i,j,k-1,ntke)) !pzhu, changed. enddo enddo Note that in dycore the vertical level is opposite to CCPP, so level npz in dycore is level 1 in ccpp. I'll ask Samuel to push this change back if you agree.
dku3d_h is used in sw_core.F90 for damping coefficient calculation. I am not sure if dku3d_h(k)=0.5*(dku3d_h(k-1)+dku3d_h(k)) is a simple solution because divergence is calculated using D-grid winds. So, I'd keep this and address the issue in future PR.

In fact, I don't know A or D grid structures. I don't understand why we need to compute tke_1(i,j)=0.5*(q(i,j,k,ntke)+q(i,j,k-1,ntke)). Since def_3 is used in tke layers (i.e., zl layer) in satmedmfvdifq.F, def_3 should also be estimated in zl layers in dycore. As I understand, all the prognostic variables (I am not sure for the vertical velocity) are located in zl layers in dycore like in CCPP physics and so, dku3d_h should be also located in the zl layers in dycore. I think divergence damping or horizontal diffusion may need a horizontal gradient of such as dku3d_h*dudx, so they may need dku3d_h at i+1, i-1, j-1, j+1, not just at i,j, but I am not quite sure.

Jongil, I'll discuss this with you and Samuel in detail when I am back. At least three of us get on the same page first. Thanks!

[jkbk2004]

All tests are done ok at ufs-community/ufs-weather-model#2734. @grantfirl @rhaesung can you merge this pr?