ESCOMP
MPAS physics-dynamics coupling issue
What happened?
CAM physics uses a mass coordinate (hydrostatic pressure), while MPAS employs a constant height formulation, where pressure is non-hydrostatic and treated as a diagnostic variable. Therefore, pressure must be computed from the MPAS prognostic state and passed to CAM physics.
In the current coupling approach, hydrostatic interface pressure is computed from density, water variables, and heights. This calculation is unambiguous. However, mid-level pressure can be determined in multiple ways. The chosen method ensures that mid-level pressure remains consistent with MPAS mid-level heights when diagnosed in CAM physics (see write-up by Lauritzen available on request).
However, this mid-level pressure is not hydrostatic and is not guaranteed to be bounded by the hydrostatic interface pressures. At high horizontal resolutions (~3.75 km), this discrepancy can result in pmid exceeding pint, causing CAM-MPAS to crash during the computation in CAM physics:
src/chemistry/mozart/mo_drydep.F90:
!-------------------------------------------------------------------------------------
! height of 1st level
!-------------------------------------------------------------------------------------
z(i) = - r/grav * air_temp(i) * (1._r8 + .61_r8*spec_hum(i)) * log(p/pg)
For the lowest model level (i=pver) it can happen that p > pg, and then log(p/pg) is positive, which will make z(i) negative. Setting:
pmid(1, iCell) = 0.5*(pint(1,iCell)+pint(2,iCell))
in dp_coupling.F90 solves the problem for this particular "instability" in CAM physics.
This is likely the issue seen in https://github.com/ESCOMP/CAM/issues/442
What are the steps to reproduce the bug?
Run ~3.75 L58 CAM-MPAS
What CAM tag were you using?
Any
What machine were you running CAM on?
CISL machine (e.g. cheyenne)
What compiler were you using?
Intel
Path to a case directory, if applicable
No response
Will you be addressing this bug yourself?
Yes
Extra info
Tagging users familiar with this issue
@briandobbins @adamrher @skamroc @kuanchihwang
What are the steps to reproduce the bug? Run ~3.75 L58 CAM-MPAS
Lol, easy enough!
Background info:
There is some flexibility in CAM physics regarding how pmid is calculated. Currently,pmid is computed such that the height z diagnosed in CAM physics (which is mass/pressure-based) matches the height in MPAS. However, this approach does not guarantee that pmid is bounded by the hydrostatic interface pressures pint provided to CAM physics.
To assess how often pmid (the non-hydrostatic mid-level pressure) falls outside the hydrostatic interface pressures pint, @briandobbins ran a ~3.75 km simulation with additional diagnostics enabled:
dp_epsilon = dp(k)*epsilon
if (pmid(k,iCell)>pint(k,iCell)-dp_epsilon.or.pmid(k,iCell)<pint(k+1,iCell)+dp_epsilon) then
write(*,*) "dbg ",k,pmid(k, iCell),pint(k,iCell),pint(k+1,iCell)
pmid(k, iCell) = 0.5*(pint(k,iCell)+pint(k+1,iCell))
end if
where epsilon=0.05. In other words, the check is triggered if pmid lies within 0.05 * dp of either bounding pint level or complete out of bounds.
In a 2-day simulation, the check was triggered at a total of 845 grid points. In 208 of these cases, pmid was found to be outside the bounds set by pint.
Note: CAM-MPAS only fails because pmid in the first level is out-of-bounds whereas the other levels does not cause a model crash; current model runs have been using:
pmid(1, iCell) = 0.5*(pint(1,iCell)+pint(2,iCell))
Conclusion: This is a rare occurance but obviously not physically acceptable. I suggest to use the code above (without the write statement) so that the mid-level pressure is correct where pmid is close or completely out of bounds but elsewhere the non-hydrostatic pmid is used so that in those points z diagnosed in CAM physics coincides with MPAS's native z coordinate.
Hi Peter, thanks for looking into this. For my understanding, the out-of-bonds pmid can occur at any vertical level, but we were only made aware of this occurring in the first level b/c it triggers the mo_drydep.F90 error? Also, what does "Column 1 Value" refer to in the histogram?