Update MSIS used for WACCM upper boundary to NRLMSIS 2.1

[fvitt]

What is the feature/what would you like to discuss?

The proposal here is to update the MSIS empirical model used for WACCM upper boundary conditions with NRLMSIS 2.1. This has a Fortran90 interface which is more supported in CESM/CAM than the Fortran77 interface in msise00.F90, in addition to scientific updates.

See: https://doi.org/10.1029/2020EA001321 https://doi.org/10.1029/2022JA030896

Is there anyone in particular you want to be part of this conversation?

@dan800 et al.

Will this change (regression test) answers?

Yes

Will you be implementing this enhancement yourself?

Yes

[dan800]

This seems like a good update and well overdue. Has the interface to MSIS changed much? I presume we start with not changing what we are asking MSIS to provide.

[fvitt]

The difference I see is the new MSIS wants altitude rather than pressure and the date format is different.

[rgarc1a]

Francis: I agree with Dan that we probably should be using the new MSIS. Having to address the model in terms of geometric altitude instead of pressure will have to be handled carefully. WACCM lev is pressure everywhere in the middle atmosphere, so the old MSIS was perfect for that. We will instead have to compute altitude (and altitude at a given pressure will change as the middle atmosphere contracts because of CO2 cooling. So, before we rush in to make change, let's have a discussion about how this might work.

[dan800]

let's have a discussion about how this might work.

Well, looking at the code, some initial thougths:

• Is it worth checking with NRL as to if a version capable of returning variables on pressure is available? I recall the old code had something like GHP as well as GTD, so both pressure and alt were options.
• It seems better to call MSISCALC() directly, rather than rely on the GTD8D "Legacy wrapper" which reconfigures the output to match the old version of MSIS.
• Consider setting lzalt_type = .false., then Z is treated as geopotential height rather than geodetic (hopefully makes for one less conversion process in the integration into WACCM).
• We should consider having the option of using MSIS NO density rather than NOEM - if it improves things then it reduces the number of files input and code complexity. It might be able to cover a broader range of solar variability.
• Is there a standard routine in CESM that moves from variables between gph (Z3) and pressure? Or do we have to code up an interpolator? Perhaps @PeterHjortLauritzen can comment?
• Is Z3 a state variable that is always available or do we need to calculate it to set the UBC. Might want to move setting UBC to where it's first calculated?
• Will this cause restart issues since the UBC now varies depending on model state, whereas before it only depended on date/UT/f107/ap.

[dan800]

Digging a little deeper - if I am reading this correctly (from nrlmsise00.F), it seems the GHP routine "converted" to pressure the output of GTD by simply making a relatively hard wired estimate of the altitude of a given pressure level.

C        PRESS - PRESSURE LEVEL(MB)
PL=ALOG10(PRESS)
C      Initial altitude estimate
      IF(PL.GE.-5.) THEN
         IF(PL.GT.2.5) ZI=18.06*(3.00-PL)
         IF(PL.GT..75.AND.PL.LE.2.5) ZI=14.98*(3.08-PL)
         IF(PL.GT.-1..AND.PL.LE..75) ZI=17.8*(2.72-PL)
         IF(PL.GT.-2..AND.PL.LE.-1.) ZI=14.28*(3.64-PL)
         IF(PL.GT.-4..AND.PL.LE.-2.) ZI=12.72*(4.32-PL)
         IF(PL.LE.-4.) ZI=25.3*(.11-PL)
         IDAY=MOD(IYD,1000)
         CL=GLAT/90.
         CL2=CL*CL
         IF(IDAY.LT.182) CD=1.-IDAY/91.25
         IF(IDAY.GE.182) CD=IDAY/91.25-3.
         CA=0
         IF(PL.GT.-1.11.AND.PL.LE.-.23) CA=1.0
         IF(PL.GT.-.23) CA=(2.79-PL)/(2.79+.23)
         IF(PL.LE.-1.11.AND.PL.GT.-3.) CA=(-2.93-PL)/(-2.93+1.11)
         Z=ZI-4.87*CL*CD*CA-1.64*CL2*CA+.31*CA*CL
      ENDIF
      IF(PL.LT.-5.) Z=22.*(PL+4.)**2+110

The Z calculated appears to calculated as a second order polynomial with latitude and a correction for day of year. So depending on the pressure level at the lid where the UBC is set and the latitude and day of year, you always getting the same Z fed to GTD. Above a certain level there's no correction for DOY or lat.

For WACCM, which I think has a top at 6e-6 hPa, alog10(PRESS) is -5.22 and Z=142.8 (For CAM-MT I can't recall if it still uses MSIS - I think we went with zero flux.)

So WACCM UBC was always GTD fed with an altitude of 142.8 km.

Seems there's room for improvement there!

[rgarc1a]

"Is Z3 a state variable that is always available": I don't know but I do know Z3 is always output on h0 files (by default?). It is generally a 3 D variable, so one would need to figure out how to use it to translate pressure (lev) into geometric altitude aa a function of latitude.

"Consider setting lzalt_type = .false., then Z is treated as geopotential height rather than geodetic (hopefully makes for one less conversion process in the integration into WACCM)" Yes, that would seem to be the way to go for WACCM.

Anyway, I suppose the first thing to do is ascertain whether Z3 is plays available.

[fvitt]

Yes Z3 is always available -- it is derived from state%zm which is updated each time step:

    do k = 1, pver
      z3(:ncol,k) = state%zm(:ncol,k) + state%phis(:ncol)*rga
    end do
    call outfld('Z3      ',z3,pcols,lchnk)

The state includes geopotential height at model layer midpoints and interfaces:

          zm        ! geopotential height above surface at midpoints (m)
          zi         ! geopotential height above surface at interfaces (m)

These are set updated via a call to geopotential_t in d_p_coupling. See https://github.com/ESCOMP/CAM/blob/cam_development/src/physics/cam/geopotential.F90

So, one should be able to use state%zi(icol,1) for the geopotential height of icol column at the upper boundary.

[rgarc1a]

Thanks, Francis. If this is the case, then we would need to figure out how to use it for the upper b.c. As I recall, the UBC in WACCM6 is a function of latitude only, so I presume this means we would need to calculate the zonal mean Z3 at the upper boundary (?)

[fvitt]

@rgarc1a No, the ubc in WACCM is not a function of latitude only -- msis_ubc is set as a function of column.

See https://github.com/ESCOMP/CAM/blob/c97e39c1b97562fe84fd4a16bdc6042c599727aa/src/chemistry/utils/mo_msis_ubc.F90#L170C1-L194C1

It seems we could set the upper BC for each column using state%zi of each column.

[rgarc1a]

OK. If so, then this should be pretty simple.

[dan800]

OK. If so, then this should be pretty simple.

Now you've cursed it... at least you didn't use the t word. Still might have issues on restarts or perhaps the setting of the UBC is always done after a calculation of Z3 and that the Z3 saved in a restart is the same one within a tilmestep that used during the execution of the run.

[fvitt]

Nothing is simple in this model. However, we should be able to deal with restart issues.

[rgarc1a]

"we should be able to deal with restart issues". Yes, I am sure this is trivial... ;)