Skip to content

Commit

Permalink
Merge branch 'hotfix-v7.2'
Browse files Browse the repository at this point in the history 
This merge addresses several issues in the MPAS-Atmosphere model.

Specifically, this merge introduces the following fixes:

- Correct the use of uninitialized memory in the init_atm_case_squall_line
  routine by initializing the qvb array to zero before its first use.

- Fix a bug in the vertical extrapolation of relative humidity and specific
  humidity to model levels below the lowest first-guess level when first-guess
  levels are given in top-to-bottom order in the input intermediate file.

- Fix reproducibility issues in several fields within the Noah LSM over land-ice
  points when running with different MPI task counts; however, only one of these
  fields -- smstav, the surface moisture availability field -- persists outside
  of the physics driver and is written to MPAS-Atmosphere restart files.

- Correct the units and description attributes for the GWDO fields var2d, con,
  oa1, oa2, oa3, oa4, ol1, ol2, ol3, and ol4 in both the init_atmosphere and
  atmosphere core Registry.xml files.

- Add code that had inadvertently been omitted for computing dtheta_dt_mp, the
  potential temperature heating rate from microphysics. Prior to this change,
  the dtheta_dt_mp field would always contain a constant zero value when written
  to model output files.

- Correct a check on the availability of the dtheta_dt_mix variable when
  computing depv_dt_mix in the PV diagnostics module.

- Correct the parallel computation of iLev_DT (and other fields that depend on
  it, including u_pv, v_pv, theta_pv, vort_pv, depv_dt_diab_pv, and
  depv_dt_fric_pv) in the PV diagnostics module.

* hotfix-v7.2:
  Correct units and descriptions for GWDO fields var2d, con, oa[1-4], ol[1-4]
  Halo exchange inTropo in flood fill to find DT.
  Fix parallel reproducibility error in 'smstav' field
  Correct check on availability of dtheta_dt_mix when computing depv_dt_mix
  Fix bug in vertical interp of r.h. and s.h. when levels are given top-to-bottom
  Initialize qvb to zero before its first use in init_atm_case_squall_line
  Remove redundant mpas_pool_get_array calls for dtheta_dt_mp computation
  Use dt_dyn rather than dt_microp when computing dtheta_dt_mp
  Remove trailing whitespace and adjust indentation of dtheta_dt_mp computation
  Obtain 'rvord' from mpas_constants for computation of dtheta_dt_mp
  Increment version number to 7.2
  Calculate dtheta_dt_mp by finite difference around microphysics call.
  • Loading branch information
@mgduda
mgduda committed Feb 14, 2022
2 parents 7557b47 + 0206851 commit 9096116
Show file tree
Hide file tree
Showing 13 changed files with 176 additions and 107 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
MPAS-v7.1
MPAS-v7.2
====

The Model for Prediction Across Scales (MPAS) is a collaborative project for
Expand Down
26 changes: 13 additions & 13 deletions src/core_atmosphere/Registry.xml
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
<?xml version="1.0"?>
<registry model="mpas" core="atmosphere" core_abbrev="atm" version="7.1">
<registry model="mpas" core="atmosphere" core_abbrev="atm" version="7.2">

<!-- **************************************************************************************** -->
<!-- ************************************** Dimensions ************************************** -->
Expand Down Expand Up @@ -3067,35 +3067,35 @@
<!-- ... PARAMETERIZATION OF GRAVITY WAVE DRAG OVER OROGRAPHY: -->
<!-- ================================================================================================== -->

<var name="var2d" type="real" dimensions="nCells" units="m^{2}"
description="variance of orography"/>
<var name="var2d" type="real" dimensions="nCells" units="m"
description="standard deviation of subgrid-scale orography"/>

<var name="con" type="real" dimensions="nCells" units="m^{2}"
description="convexity of orography"/>
<var name="con" type="real" dimensions="nCells" units="unitless"
description="orographic convexity"/>

<var name="oa1" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="asymmetry of subgrid-scale orography for westerly flow"/>

<var name="oa2" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="asymmetry of subgrid-scale orography for southerly flow"/>

<var name="oa3" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="asymmetry of subgrid-scale orography for south-westerly flow"/>

<var name="oa4" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="asymmetry of subgrid-scale orography for north-westerly flow"/>

<var name="ol1" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="effective orographic length for westerly flow"/>

<var name="ol2" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="effective orographic length for southerly flow"/>

<var name="ol3" type="real" dimensions="nCells" units="unitles"
description="directional asymmetry function of orography"/>
description="effective orographic length for south-westerly flow"/>

<var name="ol4" type="real" dimensions="nCells" units="unitless"
description="directional asymmetry function of orography"/>
description="effective orographic length for north-westerly flow"/>
</var_struct>

<var_struct name="tend_iau" time_levs="1" packages="iau">
Expand Down
3 changes: 3 additions & 0 deletions src/core_atmosphere/diagnostics/Registry_pv.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,5 +62,8 @@
<var name="iLev_DT" type="integer" dimensions="nCells Time" units="-"
description="Lowest vertical level at or above dynamic tropopause (.lt.1 if 2 PVU below column; .gt.nLevels if 2PVU above column)"/>

<var name="inTropo" type="integer" dimensions="nVertLevels nCells Time" units="0/1"
description="1 if within troposphere based on EPV flood fill"/>

</var_struct>

127 changes: 77 additions & 50 deletions src/core_atmosphere/diagnostics/pv_diagnostics.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -591,29 +591,41 @@ subroutine floodFill_tropo(mesh, diag, pvuVal)
! (2) flood fill troposphere (<2pvu) from troposphere seeds near surface.
!Somewhat paradoxically, the bottom of the stratosphere is lower than the top of the troposphere.
!Originally, it was assumed that each (MPI) domain would have >0 cells with "right" DT found by flood filling.
!However, for "small" domains over the Arctic say during winter, the entire surface can be capped by high PV.
!So, we need to communicate between domains during the flood fill or else we find the DT at the surface.
!The extreme limiting case is if we had every cell as its own domain; then, it's clear that there has to be communication.

!The "output" is iLev_DT, which is the vertical index for the level >= pvuVal. If >nVertLevels, pvuVal above column. If <2, pvuVal below column.
!Communication between blocks during the flood fill may be needed to treat some edge cases appropriately.

use mpas_pool_routines, only : mpas_pool_get_dimension, mpas_pool_get_array
use mpas_pool_routines, only : mpas_pool_get_dimension, mpas_pool_get_array, mpas_pool_get_field
use mpas_dmpar, only : mpas_dmpar_max_int,mpas_dmpar_exch_halo_field
use mpas_derived_types, only : dm_info, field2DInteger

implicit none

type (mpas_pool_type), intent(in) :: mesh
type (mpas_pool_type), intent(inout) :: diag
real(kind=RKIND), intent(in) :: pvuVal
integer :: iCell, k, nChanged, iNbr, iCellNbr, levInd
integer, pointer :: nCells, nVertLevels

integer :: iCell, k, nChanged, iNbr, iCellNbr, levInd, haloChanged, global_haloChanged
integer, pointer :: nCells, nVertLevels, nCellsSolve
integer, dimension(:), pointer :: nEdgesOnCell, iLev_DT
integer, dimension(:,:), pointer :: cellsOnCell
integer, dimension(:,:), pointer :: cellsOnCell, inTropo

type (field2DInteger), pointer :: inTropo_f

real(kind=RKIND) :: sgnHemi, sgn
real(kind=RKIND),dimension(:),pointer:: latCell
real(kind=RKIND), dimension(:,:), pointer :: ertel_pv

integer, dimension(:,:), allocatable :: candInTropo, inTropo !whether in troposphere
type (dm_info), pointer :: dminfo

integer, dimension(:,:), allocatable :: candInTropo !whether in troposphere

call mpas_pool_get_dimension(mesh, 'nCells', nCells)
call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
Expand All @@ -622,9 +634,9 @@ subroutine floodFill_tropo(mesh, diag, pvuVal)
call mpas_pool_get_array(diag, 'ertel_pv', ertel_pv)
!call mpas_pool_get_array(diag, 'iLev_DT_trop', iLev_DT)
call mpas_pool_get_array(diag, 'iLev_DT', iLev_DT)
call mpas_pool_get_array(diag, 'inTropo', inTropo)

allocate(candInTropo(nVertLevels, nCells+1))
allocate(inTropo(nVertLevels, nCells+1))
candInTropo(:,:) = 0
inTropo(:,:) = 0
!store whether each level above DT to avoid repeating logic. we'll use cand as a isVisited marker further below.
Expand All @@ -645,56 +657,71 @@ subroutine floodFill_tropo(mesh, diag, pvuVal)
do k=1,levInd
if (candInTropo(k,iCell) .GT. 0) then
inTropo(k,iCell) = 1
candInTropo(k,iCell) = 0
!candInTropo(k,iCell) = 0
nChanged = nChanged+1
end if
end do
end do

!flood fill from the given seeds. since I don't know enough fortran,
!we'll just brute force a continuing loop rather than queue.
do while(nChanged .GT. 0)
nChanged = 0
do iCell=1,nCells
do k=1,nVertLevels
!update if candidate and neighbor in troposphere
if (candInTropo(k,iCell) .GT. 0) then
!nbr below
if (k .GT. 1) then
if (inTropo(k-1,iCell) .GT. 0) then
inTropo(k,iCell) = 1
candInTropo(k,iCell) = 0
nChanged = nChanged+1
cycle
end if
end if
!side nbrs
do iNbr = 1, nEdgesOnCell(iCell)
iCellNbr = cellsOnCell(iNbr,iCell)
if (inTropo(k,iCellNbr) .GT. 0) then
inTropo(k,iCell) = 1
candInTropo(k,iCell) = 0
nChanged = nChanged+1
cycle
call mpas_pool_get_field(diag, 'inTropo', inTropo_f)
dminfo => inTropo_f % block % domain % dminfo
global_haloChanged = 1
do while(global_haloChanged .GT. 0) !any cell in a halo has changed, to propagate to other domains
global_haloChanged = 0 !aggregate the number of changed cells w/in the loop below
do while(nChanged .GT. 0)
nChanged = 0
do iCell=1,nCells !should we look for neighbors of hallo cells?
!do iCell=1,nCellsSolve !should we look for neighbors of hallo cells?
do k=1,nVertLevels
!update if candidate and neighbor in troposphere
if ((candInTropo(k,iCell) .GT. 0) .AND. (inTropo(k,iCell).LT.1) ) then
!nbr below
if (k .GT. 1) then
if (inTropo(k-1,iCell) .GT. 0) then
inTropo(k,iCell) = 1
!candInTropo(k,iCell) = 0
nChanged = nChanged+1
cycle
end if
end if
end do
!nbr above
if (k .LT. nVertLevels) then
if (inTropo(k+1,iCell) .GT. 0) then
inTropo(k,iCell) = 1
candInTropo(k,iCell) = 0
nChanged = nChanged+1
cycle

!side nbrs
do iNbr = 1, nEdgesOnCell(iCell)
iCellNbr = cellsOnCell(iNbr,iCell)
if (inTropo(k,iCellNbr) .GT. 0) then
inTropo(k,iCell) = 1
!candInTropo(k,iCell) = 0
nChanged = nChanged+1
exit
end if
end do

!nbr above
if (k .LT. nVertLevels) then
if (inTropo(k+1,iCell) .GT. 0) then
inTropo(k,iCell) = 1
!candInTropo(k,iCell) = 0
nChanged = nChanged+1
cycle
end if
end if
end if
end if !candIn
end do !levels
end do !cells
!here's where a communication would be needed for edge cases !!!
end do !while

end if !candIn
end do !levels
end do !cells
global_haloChanged = global_haloChanged+nChanged
end do !while w/in domain
!communicate to other domains for edge case where a chunk of a block hasn't gotten to fill
nChanged = global_haloChanged
call mpas_dmpar_max_int(dminfo, nChanged, global_haloChanged)
if (global_haloChanged .GT. 0) then !communicate inTropo everywhere
call mpas_dmpar_exch_halo_field(inTropo_f)
end if
nChanged = global_haloChanged !so each block will iterate again if anything changed
end do !while haloChanged
deallocate(candInTropo)
!Fill iLev_DT with the lowest level above the tropopause (If DT above column, iLev>nVertLevels. If DT below column, iLev=0.
do iCell=1,nCells
Expand Down Expand Up @@ -1457,7 +1484,7 @@ subroutine calc_pvBudget(state, time_lev, diag, mesh, tend, tend_physics)
depv_dt_mp(k,iCell) = 0.0_RKIND
end if
if (associated(dtheta_dt_mp)) then
if (associated(dtheta_dt_mix)) then
call calc_grad_cell(gradtheta, &
iCell, k, nVertLevels, nEdgesOnCell(iCell), verticesOnCell, kiteAreasOnVertex, &
cellsOnCell, edgesOnCell, cellsOnEdge, dvEdge, edgeNormalVectors, &
Expand Down
Loading

0 comments on commit 9096116

Please sign in to comment.