Mismatch of vertical dimensions for ak and bk

[climbfuji]

Description

This bug only affects model runs with different radiation and physics levels. I noticed it when using the newly developed debugging capability in CCPP that checks the dimensions of arrays in the metadata against the actual dimensions of the array.

Variables ak and bk in GFS_typedefs.F90 are provided and used as follows, which at first glance suggests that the dimensionality of these arrays is vertical_interface_dimension_for_radiation = Model%levr + 1:

    allocate(Model%ak(1:size(ak)))
    allocate(Model%bk(1:size(bk)))
    Model%ak               = ak
    Model%bk               = bk
    ...
    allocate(Model%si(Model%levr+1))
    !--- Define sigma level for radiation initialization
    !--- The formula converting hybrid sigma pressure coefficients to sigma coefficients follows Eckermann (2009, MWR)
    !--- ps is replaced with p0. The value of p0 uses that in http://www.emc.ncep.noaa.gov/officenotes/newernotes/on461.pdf
    !--- ak/bk have been flipped from their original FV3 orientation and are defined sfc -> toa
    Model%si = (ak + bk * con_p0 - ak(Model%levr+1)) / (con_p0 - ak(Model%levr+1))

The ak and bk variables are calculated in subroutine atmosphere_etalvls in the FV3 dycore:

 subroutine atmosphere_etalvls (ak, bk, flip)
   real(kind=kind_phys), pointer, dimension(:), intent(inout) :: ak, bk
   logical, intent(in) :: flip !< control vertical index flipping

   allocate(ak(npz+1))
   allocate(bk(npz+1))

   if (flip) then
     ak(1:npz+1) = Atm(mygrid)%ak(npz+1:1:-1)
     bk(1:npz+1) = Atm(mygrid)%bk(npz+1:1:-1)
   else
     ak(1:npz+1) = Atm(mygrid)%ak(1:npz+1)
     bk(1:npz+1) = Atm(mygrid)%bk(1:npz+1)
   endif
 end subroutine atmosphere_etalvls

Here, npz is the number of vertical model levels that are the vertical levels for all physics except radiation, npz = Model%levs. Thus, the correct dimension in the CCPP metadata would be vertical_interface_dimension. This can be changed easily, but I want to make sure the following:

For model runs with fewer radiation levels (levr) than physics levels (levs), is it true that the radiation levels are exactly the same as the physics levels up to the point levr and that the levels above are simply ignored by the radiation?

If so, should the calculation of Model%si be written explicitly as follows?

    Model%si(1:Model%levr+1) = (ak(1:Model%levr+1) + bk(1:Model%levr+1) * con_p0 - ak(Model%levr+1)) / (con_p0 - ak(Model%levr+1))

And lastly, can we safely assume that there will never be any model runs where there are more radiation levels (excluding the extra layers at top LTP, those are treated differently), i.e. should we add a guard that prevents Model%levr > Model%levs?

@DusanJovic-NOAA @junwang-noaa @yangfanglin @SMoorthi-emc @SajalKar-NOAA can you clarify please?

[yangfanglin]

It is correct "there will never be any model runs where there are more radiation levels (excluding the extra layers at top LTP, those are treated differently)"

Yes, for model runs with fewer radiation levels (levr) than physics levels (levs), it is true that the radiation levels are exactly the same as the physics levels up to the point levr and that the levels above are simply ignored by the radiation.

Model%si(1:Model%levr+1) should be defined for the entire atmosphere for all layers, not only for layers where radiation operates.

[climbfuji]

It is correct "there will never be any model runs where there are more radiation levels (excluding the extra layers at top LTP, those are treated differently)"

Yes, for model runs with fewer radiation levels (levr) than physics levels (levs), it is true that the radiation levels are exactly the same as the physics levels up to the point levr and that the levels above are simply ignored by the radiation.

Model%si(1:Model%levr+1) should be defined for the entire atmosphere for all layers, not only for layers where radiation operates.

Thanks, Fanglin, this is helpful. I will make those changes in a forthcoming PR.

[SMoorthi-emc]

Just to be clear, pressure corresponding to "Model%levr+1" should be the model top. So, for the radiation calculation, the top layer goes from Model%levr to model top i.e. Model%levs+1. Moorthi

…

On Mon, Nov 22, 2021 at 1:03 PM Dom Heinzeller ***@***.***> wrote: It is correct "there will never be any model runs where there are more radiation levels (excluding the extra layers at top LTP, those are treated differently)" Yes, for model runs with fewer radiation levels (levr) than physics levels (levs), it is true that the radiation levels are exactly the same as the physics levels up to the point levr and that the levels above are simply ignored by the radiation. Model%si(1:Model%levr+1) should be defined for the entire atmosphere for all layers, not only for layers where radiation operates. Thanks, Fanglin, this is helpful. I will make those changes in a forthcoming PR. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#430 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALLVRYWX7R62QZTB4Y77VMTUNKAWXANCNFSM5IRRGOMQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

-- Dr. Shrinivas Moorthi Research Meteorologist Modeling and Data Assimilation Branch Environmental Modeling Center / National Centers for Environmental Prediction 5830 University Research Court - (W/NP23), College Park MD 20740 USA Tel: (301)683-3718 e-mail: ***@***.*** Phone: (301) 683-3718 Fax: (301) 683-3718

[climbfuji]

@SMoorthi-emc This is something that needs to be checked - I don't think this is the case, but I could be wrong

[climbfuji]

Fixed in #431