Merge branch 'hotfix-v7.2'

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.
MPAS-Dev · Feb 14, 2022 · 9096116 · 9096116
2 parents 7557b47 + 0206851
commit 9096116
Show file tree

Hide file tree

Showing 13 changed files with 176 additions and 107 deletions.
diff --git a/README.md b/README.md
@@ -1,4 +1,4 @@
-MPAS-v7.1
+MPAS-v7.2
 ====
 
 The Model for Prediction Across Scales (MPAS) is a collaborative project for

diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
@@ -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 ************************************** -->
@@ -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">

diff --git a/src/core_atmosphere/diagnostics/Registry_pv.xml b/src/core_atmosphere/diagnostics/Registry_pv.xml
@@ -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>
 
diff --git a/src/core_atmosphere/diagnostics/pv_diagnostics.F b/src/core_atmosphere/diagnostics/pv_diagnostics.F
@@ -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)
@@ -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.
@@ -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
@@ -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, &