Add support for staggered atmospheric levels #364

phil-blain · 2021-05-26T20:59:30Z

PR checklist

I also refactored the code a bit to make it more easily understandable.

This was implemented in our in-house versions of CICE4 and CICE5 and I'm porting it over to the Consortium as part of our project of updating our systems to CICE6.

The two functions 'psimhu' and 'psixhu' in icepack_atmo compute the stability functions in the unstable case for momentum and scalars (heat and water), respectively. They are both function of X := (1 - 16 \zeta)^(1/4) (see [1], section 8.1.1). In the code X is computed first ('xqq') and passed to 'psimhu' and 'psixhu'. Add two new functions with more explicit names, 'psi_momentum_unstable' and 'psi_scalar_unstable', that take \zeta ('hol') as argument, compute X internally and return the evaluated stability function. A subsequent commit will refactor the code to use them. [1] B.G. Kauffman and W.G. Large. The CCSM coupler, version 5.0.1. 2002. URL: https://github.com/CICE-Consortium/CICE/blob/master/doc/PDF/KL_NCAR2002.pdf.

Use the new functions 'psi_momentum_unstable' and 'psi_scalar_unstable' introduced in the preceding commit to compute the stability function evaluated at 'hol'. Remove duplicated computation of 'xqq'. Also remove the now unused functions 'psimhu' and 'psixhu'.

In a subsequent commit, we will support receiving the momentum and scalar atmospheric variables at different levels. In order to minimize code duplication, factor out the computation of the stability parameter 'hol' and the stability function 'psi[mx]h' in function 'compute_stability_parameter' and subroutine 'compute_stability_function'. 'compute_stability_function' also makes available the unit step function evaluated at 'hol', 'stable', since it is used at the end of 'atmo_boundary_layer' to compute the diagnostic temperature and humidity.

In some applications, we might want the atmospheric input variables (wind, temperature and humidity) to be given at different vertical levels. For example, we might want to receive the winds at some level and the scalar variables (temperature and humidity) at another. To support this use case, add an optional argument 'zlvs' to subroutine 'icepack_atm_boundary', and pass it down to 'atmo_boundary_layer' also as an optional argument, so that existing applications do not have to be changed if they do not want to use this feature. The name 'zlvs' is chosen to mirror the existing variable 'zlvl', but with the suffix 's' for 'scalars', mimicking the existing convention for several variables in 'atmo_boundary_layer'. In 'atmo_boundary_layer', add new local variables 'alzs' and 'hols' to hold quantities computed at zlvs. Also, for clarity, rename 'alz' to 'alzm' and 'hol' to 'holm', where the 'm' suffix stands for 'momentum', again mimicking existing convention in the subroutine. If 'zlvs' is not passed to 'atmo_boundary_layer', it is assumed that all variables are given at the same level (zlvl), so in that case there is no change in the results.

The previous commit added an optional 'zlvs' argument to 'icepack_atm_boundary' to allow atmospheric input variables to be given at different vertical levels. Also add 'zlvs' as an optional argument to 'icepack_step_therm1', to make the functionality available via the Icepack interfaces.

eclare108213

This looks okay to me, just reading through it. Can you add a note about it in the text of the documentation, explaining how it works?

phil-blain · 2021-05-27T14:17:41Z

@eclare108213 yes I'll do that.

phil-blain · 2021-05-27T18:24:02Z

I've just pushed a new version of the last commit with some modifications to the Science guide. See the build here: https://cice-consortium-icepack--364.org.readthedocs.build/en/364/science_guide/sg_boundary_forcing.html#atmosphere

apcraig

There may be two places where an optional argument is passed into a subroutine without checking if that argument was passed into the calling subroutine. This is one of the major issues with the current implementation of icepack. My recommendation has been to create a local variable, set it to a default if the optional variable is not passed, then pass that local variable into other subroutines. Something like

real (kind=dbl_kind) :: lzlvs

if (present(zlvs)) then
  lzlvs = zlvs
else
  lzlvs = zlvl
endif

call subroutine(..., zlvs=lzlvs, ....)

If there are better ways to handle this, I'd love to know! I don't think the "optional" attribute is carried with the argument as you try to pass it down the calling tree, but maybe there is a way to do this?

apcraig · 2021-05-27T19:54:23Z

columnphysics/icepack_therm_vertical.F90

@@ -2539,7 +2541,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
                                        Qa_iso=l_Qa_iso,           &
                                        Qref_iso=Qrefn_iso,      &
                                        uvel=uvel, vvel=vvel,    &
-                                        Uref=Urefn)
+                                        Uref=Urefn, zlvs=zlvs)


It looks like passing zlvs here could be a problem if zlvs is not passed into icepack_step_therm1.

apcraig · 2021-05-27T19:55:13Z

columnphysics/icepack_atmo.F90

@@ -965,7 +969,7 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                   Qa_iso=l_Qa_iso,         &
                                   Qref_iso=l_Qref_iso,     &
                                   uvel=l_uvel, vvel=l_vvel,  &
-                                   Uref=l_Uref)
+                                   Uref=l_Uref, zlvs=zlvs)


It looks like passing zlvs here could be a problem if zlvs is not passed into icepack_atm_boundary.

phil-blain · 2021-05-27T20:09:43Z

Yes, I thought about that myself when writing the code and checked my copy of Modern Fortran Explained. It is said that the optional attribute is carried forward though the call stack. That is why I did not use local variables as you suggest...

Section 5.13 "Keyword and optional arguments" (p .89):

In such a case, an absent optional argument is also regarded as absent in the second-level subprogram. [...] Such absent arguments may be propagated through any number of calls, provided the dummy argument is optional each case

apcraig · 2021-05-27T20:16:37Z

Interesting. I think we've run into problems before not checking optional arguments, but maybe I'm not remembering correct. (See for example Qa_iso in icepack_atmo.F90). If this is robust, that's going to be very helpful moving forward. Maybe Qa_iso creates different problems because it's an array?

phil-blain · 2021-05-27T20:27:30Z

Maybe some compiler implementations you were using are not conforming... might be worth fetching this branch and running the tests suite on more machines ? It builds and runs fine on Intel and Gfortran...

eclare108213 · 2021-05-27T22:39:49Z

doc/source/science_guide/sg_boundary_forcing.rst

@@ -143,6 +144,10 @@ respectively, given the ice velocity :math:`\vec{U}_i`:
   \end{aligned}
   :label: stars

+
+Note that *atmo_boundary_layer* also accepts an optional argument, ``zlvs``, to support receiving scalar quantities from the atmospheric model (humidity and temperature)
+at a different vertical level than the winds. In that case a separate stability :math:`\Upsilon_s` is computed using the same formula as above but substituting :math:`z_o` by :math:`z_{o,s}`.


The "staggered" terminology makes sense here in the way that you use it, but it would be helpful to just add (staggered) after at a different vertical level than the winds, to more directly connect the comment below with this one. Optional.

yes, good idea.

I've just pushed an updated commit that now reads

Note that atmo_boundary_layer also accepts an optional argument, zlvs, to support staggered atmospheric levels, i.e. receiving scalar quantities from the atmospheric model (humidity and temperature)at a different vertical level than the winds.

eclare108213

Was interfaces.include.orig committed intentionally? Looks like it's a personal copy.

phil-blain · 2021-05-28T16:52:00Z

Oups. I did git add doc so it was automatically included... it's not a personal copy, it's generated by the gen_interfaces.sh script and not cleaned up at the end...

Icepack/doc/generate_interfaces.sh

Line 42 in eb59c95

mv ${rstfile} ${rstfile}.orig

I'll re-push without it.

Add an index entry for 'zlvs' and mention it in the 'Atmopshere' section of the science guide.

apcraig · 2021-05-28T17:11:30Z

I ran full test suites of both Icepack and CICE with the staggered atmosphere levels and everything works fine. If we can leverage the fact that the optional attribute is carried thru the code for future developments, that is very helpful. There may be cases where we can't use it, but when we can, we should. Thanks @phil-blain.
https://github.com/CICE-Consortium/Test-Results/wiki/icepack_by_hash_forks#a3824f10e3510615067f91e6c8ee97dddb9c2b24
https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#d71069b2bc51ab5d4fb32ca22a35d453ebedd377

When we compute the diagnostic temperature and humidity (at zTrf) at the end of icepack_atmo::atmo_boundary_layer, we shift the stability parameter (hols) by multiplying it by zTrf/zlvl, but this is incorrect when staggered atmospheric levels are used, as zlvl is the atmospheric level height for momentum, not for scalars (temperature and humidity). This is a bug introduced in 37f2a17 (Add support for staggered atmospheric levels (CICE-Consortium#364), 2021-05-28) but only changes diagnostics if staggered atmospheric levels are used, i.e. if zlvs /= zlvl. Correct that by dividing by zlvs, the atmospheric level height for scalars, when it is present. While at it, adjust the indentation of the description comment for 'zlvs', so it is aligned with the other variables.

In 37f2a17 (Add support for staggered atmospheric levels (CICE-Consortium#364), 2021-05-28), we introduced the subroutine compute_stability_function, which computes the stability function during the iteration in 'icepack_amto::atmo_boundary_layer'. This subroutine also returns the unit step function 'stable' at 'hol', since it is used at the end of the computation to compute the diagnostic temperature and humidity. Note that we overwrite 'stable' in the second call to compute_stability_function, but this is OK since we call it first for momentum and then for scalars, and it is the value for scalars that we want to use later for diagnostic temperature and humdity. In a following commit we will want to compute the diagnostic _velocity_ in a different way than what is done now, and we will need the unit step function computed at the momentum level (which is currently overwritten). In preparation for that change, add a new variable 'stablem' for the unit step function at the momentum level, and rename the existing 'stable' variable to 'stables' (for scalars). Adjust the computation of 'psix2' accordingly.