GitFlow: Merge develop into main

.circleci/config.yml

@@ -1,7 +1,7 @@
 version: 2.1
 
 # Anchors in case we need to override the defaults from the orb
-#baselibs_version: &baselibs_version v8.9.0
+#baselibs_version: &baselibs_version v8.32.0
 #bcs_version: &bcs_version v12.0.0
 
 orbs:
@@ -20,10 +20,7 @@ workflows:
           #baselibs_version: *baselibs_version
           repo: GEOSgcm
           checkout_fixture: true
-          # V12 code uses a special branch for now.
-          fixture_branch: feature/sdrabenh/gcm_v12
-          # We comment out this as it will "undo" the fixture_branch
-          #mepodevelop: false
+          mepodevelop: true
           persist_workspace: true # Needs to be true to run fv3/gcm experiment, costs extra
       - ci/run_fv3:
           name: run-FV3-on-<< matrix.compiler >>-with-GEOSgcm
@@ -60,8 +57,6 @@ workflows:
           #baselibs_version: *baselibs_version
           repo: GEOSfvdycore
           checkout_fixture: true
-          # V12 code uses a special branch for now.
-          fixture_branch: feature/sdrabenh/gcm_v12
           mepodevelop: false
           persist_workspace: true # Needs to be true to run fv3/gcm experiment, costs extra
       - ci/run_fv3:

.github/workflows/workflow.yml

@@ -23,13 +23,12 @@ jobs:
       matrix:
         compiler: [ifort, gfortran-14, gfortran-15, ifx]
         build-type: [Debug]
+      fail-fast: false
     uses: GEOS-ESM/CI-workflows/.github/workflows/geosgcm_build_tests.yml@project/geosgcm
     with:
       compiler: ${{ matrix.compiler }}
       cmake-build-type: ${{ matrix.build-type }}
-      fail-fast: false
       fixture-repo: GEOS-ESM/GEOSgcm
-      fixture-ref: feature/sdrabenh/gcm_v12
 
   spack_build:
     uses: GEOS-ESM/CI-workflows/.github/workflows/spack_gcc_build.yml@project/geosgcm
@@ -38,5 +37,5 @@ jobs:
       BUILDCACHE_TOKEN: ${{ secrets.BUILDCACHE_TOKEN }}
     with:
       fixture-repo: GEOS-ESM/GEOSgcm
-      fixture-ref: feature/sdrabenh/gcm_v12
+      load-fms: true
 

DynCore_GridCompMod.F90

@@ -51,7 +51,8 @@ Module FVdycoreCubed_GridComp
                            DYN_DEBUG       => DEBUG,                 &
                            HYDROSTATIC     => FV_HYDROSTATIC,        &
                            fv_getUpdraftHelicity, DEBUG_DYN,  DEBUG_ADV, &
-                           ADIABATIC, SW_DYNAMICS, AdvCore_Advection
+                           ADIABATIC, SW_DYNAMICS, AdvCore_Advection, &
+                           fv_getGradT_2D
    use m_topo_remap, only: dyn_topo_remap
    use CubeGridPrototype, only: register_grid_and_regridders
 
@@ -1360,6 +1361,14 @@ Subroutine SetServices ( gc, rc )
        VLOCATION          = MAPL_VLocationNone,          RC=STATUS )
     VERIFY_(STATUS)
 
+    call MAPL_AddExportSpec ( gc,                                  &
+       SHORT_NAME         = 'WSPD_STABLE300M',                     &
+       LONG_NAME          = 'max_wind_speed_in_stable_cold_surface_layer', &
+       UNITS              = 'm s-1',                               &
+       DIMS               = MAPL_DimsHorzOnly,                     &
+       VLOCATION          = MAPL_VLocationNone,          RC=STATUS )
+    VERIFY_(STATUS)
+
     call MAPL_AddExportSpec ( gc,                                  &
        SHORT_NAME         = 'VVEL_UP_100_1000',                    &
        LONG_NAME          = 'max_vertical_velocity_up_between_100_1000_hPa', &
@@ -2992,10 +3001,14 @@ subroutine Run(gc, import, export, clock, rc)
     integer  :: NKE, NPHI
     integer  :: NUMVARS
     integer  :: ifirstxy, ilastxy, jfirstxy, jlastxy
-    integer  :: kend, i, j, K, L, n
+    integer  :: kend, i, j, K, L, n, ii, jj
     integer  :: im_replay,jm_replay
     logical, parameter :: convt = .false. ! Until this is run with full physics
     logical  :: is_shutoff, is_ringing
+    logical  :: is_stable
+
+    integer  :: k_bot, k_top
+    real     :: z_agl, du, dv
 
     real(r8),     pointer :: phisxy(:,:)
     real(kind=4), pointer ::   phis(:,:)
@@ -5120,6 +5133,12 @@ subroutine Run(gc, import, export, clock, rc)
 
 ! Fill Surface and Near-Surface Variables
 ! ----------------------------------------------
+
+    ! Get the height above the surface
+      do k=1,km+1
+         zle(:,:,k) = zle(:,:,k) - zle(:,:,km+1)
+      enddo
+
                    HGT_SURFACE = 50.0
    if (km .eq. 72) HGT_SURFACE =  0.0
    call MAPL_GetResource ( MAPL, HGT_SURFACE, Label="HGT_SURFACE:", DEFAULT=HGT_SURFACE, RC=STATUS)
@@ -5131,11 +5150,6 @@ subroutine Run(gc, import, export, clock, rc)
       VERIFY_(STATUS)
       if(associated(temp2d)) temp2d = HGT_SURFACE
 
-    ! Get the height above the surface
-      do k=1,km+1
-         zle(:,:,k) = zle(:,:,k) - zle(:,:,km+1)
-      enddo
-
       call MAPL_GetPointer(export,temp2d,'PS',  rc=status)
       VERIFY_(STATUS)
       if(associated(temp2d)) temp2d =  vars%pe(:,:,km+1)
@@ -5218,6 +5232,45 @@ subroutine Run(gc, import, export, clock, rc)
        VERIFY_(STATUS)
    end if
 
+   call MAPL_GetPointer(export,temp2d,'WSPD_STABLE300M',rc=status)           
+   VERIFY_(STATUS)
+    if(associated(temp2d)) then
+       tempxy = vars%pt * vars%pkz
+       temp2d = 0.0
+       do j = jfirstxy, jlastxy
+          jj = j - jfirstxy + 1
+          do i = ifirstxy, ilastxy
+             ii = i - ifirstxy + 1  
+             ! 1. Check if surface air is freezing
+             if (tempxy(i,j,km) <= MAPL_TICE) then
+                ! Assume no inversion until proven otherwise
+                is_stable = .false.
+                ! Start max wind tracking with the lowest model level
+                temp2d(ii,jj) = SQRT(ur(i,j,km)**2 + vr(i,j,km)**2)
+                ! 2. Scan the lowest 300m
+                do k = km-1, 1, -1
+                   ! Height AGL
+                   if ( (0.5 * (zle(i,j,k) + zle(i,j,k+1)) - zle(i,j,km+1)) <= 300.0 ) then
+                      ! Track maximum wind speed
+                      temp2d(ii,jj) = MAX(temp2d(ii,jj), SQRT(ur(i,j,k)**2 + vr(i,j,k)**2))
+                      ! 3. Check for Inversion ANYWHERE in the 300m layer
+                      if (tempxy(i,j,k) > tempxy(i,j,km)) then
+                         is_stable = .true.
+                      endif
+                   else
+                      exit ! Reached top of 300m layer
+                   endif
+                end do
+                ! 4. If no inversion was found in the whole 300m, it's not a katabatic zone. 
+                ! Zero out the wind speed.
+                if (.not. is_stable) then
+                   temp2d(ii,jj) = 0.0
+                endif
+             endif
+          end do
+       end do
+    end if
+
    if (.not. HYDROSTATIC) then
    call MAPL_GetPointer(export,temp2d,'VVEL_UP_100_1000',rc=status)
    VERIFY_(STATUS)
@@ -6485,15 +6538,22 @@ subroutine RunAddIncs(gc, import, export, clock, rc)
     real(r8), allocatable ::  logpe(:,:,:)
     real(r8), allocatable ::  logps(:,:)
 
+    real(r4), allocatable :: t_600(:,:), ps(:,:)
+
     real(FVPRC)              :: dt
 
     real(r4), pointer     :: QOLD(:,:,:)
     real(r4), pointer     :: temp3d(:,:,:)
     real(r4), pointer     :: temp2d(:,:  )
 
+    integer isc,iec, jsc,jec
+    integer isd,ied, jsd,jed
     integer ifirstxy, ilastxy
     integer jfirstxy, jlastxy
     integer im,jm,km, iNXQ
+    integer ii,jj
+    real :: p_layer
+
     real(r4), pointer     :: ztemp1(:,:  )
     real(r4), pointer     :: ztemp2(:,:  )
     real(r4), pointer     :: ztemp3(:,:  )
@@ -6532,6 +6592,16 @@ subroutine RunAddIncs(gc, import, export, clock, rc)
     grid  => state%grid   ! direct handle to grid
     dt    =  state%dt     ! dynamics time step (large)
 
+    isc = grid%is
+    iec = grid%ie
+    jsc = grid%js
+    jec = grid%je
+
+    isd = grid%isd
+    ied = grid%ied
+    jsd = grid%jsd
+    jed = grid%jed
+
     ifirstxy = grid%is
     ilastxy  = grid%ie
     jfirstxy = grid%js
@@ -6771,7 +6841,6 @@ subroutine RunAddIncs(gc, import, export, clock, rc)
     VERIFY_(STATUS)
     if(associated(temp3d)) temp3d = 0.5*( zle(:,:,2:) + zle(:,:,:km) )
 
-
 ! Fill Single Level Variables
 ! ---------------------------
 

FV_StateMod.F90

@@ -85,6 +85,7 @@ module FV_StateMod
   public T_TRACERS, T_FVDYCORE_VARS, T_FVDYCORE_GRID, T_FVDYCORE_STATE
   public fv_fillMassFluxes
   public fv_computeMassFluxes
+  public fv_getGradT_2D, fv_getGradT_3D
   public fv_getVerticalMassFlux
   public fv_getPK
   public fv_getOmega
@@ -2820,6 +2821,114 @@ subroutine adjust_fv3_splits(state, isc, iec, jsc, jec, isd, ied, jsd, jed, npz,
 
 end subroutine adjust_fv3_splits
 
+subroutine fv_getGradT_2D(t_in, ps_in, mag_grad_t)
+  ! INPUT: 2D Temp (REAL4) interpolated to pressure level (with halo bounds)
+  real(REAL4), intent(INOUT) :: t_in(FV_Atm(1)%bd%isd:FV_Atm(1)%bd%ied, &
+                                     FV_Atm(1)%bd%jsd:FV_Atm(1)%bd%jed)
+  real(REAL4), intent(INOUT) :: ps_in(FV_Atm(1)%bd%isd:FV_Atm(1)%bd%ied, &
+                                      FV_Atm(1)%bd%jsd:FV_Atm(1)%bd%jed)
+
+  ! OUTPUT: MAPL Pointer (REAL4) (Assumed shape, 1-based indexing: 1:im, 1:jm)
+  real(REAL4), intent(OUT)   :: mag_grad_t(:,:)
+
+  integer :: isc, iec, jsc, jec
+  integer :: i, j, ii, jj
+  real(REAL4) :: dTdx, dTdy, p_surf_min
+  real(REAL8) :: dx, dy
+
+  isc = FV_Atm(1)%bd%isc ; iec = FV_Atm(1)%bd%iec
+  jsc = FV_Atm(1)%bd%jsc ; jec = FV_Atm(1)%bd%jec
+
+  ! 1. Update halos. FV3 mpp_update_domains natively supports REAL4 arrays.
+  call mpp_update_domains(ps_in, FV_Atm(1)%domain, complete=.false.)
+  call mpp_update_domains( t_in, FV_Atm(1)%domain, complete=.true.)
+
+  ! 2. Calculate gradient, avoid UNDEF, and map to 1-based MAPL array
+  !$OMP PARALLEL DO DEFAULT(NONE) &
+  !$OMP SHARED(FV_Atm, t_in, ps_in, mag_grad_t, isc, iec, jsc, jec) &
+  !$OMP PRIVATE(i, j, ii, jj, dx, dy, dTdx, dTdy, p_surf_min)
+  do j = jsc, jec
+     jj = j - jsc + 1   ! Map to MAPL pointer index
+
+     do i = isc, iec
+        ii = i - isc + 1  ! Map to MAPL pointer index
+
+        ! Find the highest terrain in the 5-point stencil
+        p_surf_min = MIN(ps_in(i,j), ps_in(i-1,j), ps_in(i+1,j), ps_in(i,j-1), ps_in(i,j+1))
+
+        ! SHIELD: Only calculate if the whole stencil is safely > 700 hPa
+        ! (Keeps the 600 hPa level well out of the mountain boundary layer)
+        if ( p_surf_min > 70000.0_REAL4 .and. &
+             t_in(i,j)   /= real(MAPL_UNDEF, REAL4) .and. &
+             t_in(i-1,j) /= real(MAPL_UNDEF, REAL4) .and. &
+             t_in(i+1,j) /= real(MAPL_UNDEF, REAL4) .and. &
+             t_in(i,j-1) /= real(MAPL_UNDEF, REAL4) .and. &
+             t_in(i,j+1) /= real(MAPL_UNDEF, REAL4) ) then
+
+           dx = FV_Atm(1)%gridstruct%dxa(i,j)
+           dy = FV_Atm(1)%gridstruct%dya(i,j)
+
+           ! Centered difference
+           dTdx = (t_in(i+1,j) - t_in(i-1,j)) / real(2.0_REAL8 * dx, REAL4)
+           dTdy = (t_in(i,j+1) - t_in(i,j-1)) / real(2.0_REAL8 * dy, REAL4)
+
+           mag_grad_t(ii,jj) = SQRT(dTdx**2 + dTdy**2)
+
+        else
+           ! Stencil hits topography or uninitialized data
+           mag_grad_t(ii,jj) = 0.0_REAL4
+        endif
+
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+return
+end subroutine fv_getGradT_2D
+
+subroutine fv_getGradT_3D(mag_grad_t)
+  ! OUTPUT: Magnitude of the horizontal temperature gradient
+  real(REAL8), intent(OUT) :: mag_grad_t(FV_Atm(1)%bd%isc:FV_Atm(1)%bd%iec, &
+                                         FV_Atm(1)%bd%jsc:FV_Atm(1)%bd%jec, &
+                                         1:FV_Atm(1)%npz)
+
+  integer :: isc, iec, jsc, jec
+  integer :: i, j, k
+  real(REAL8) :: dTdx, dTdy
+  real(REAL8) :: dx, dy
+
+  isc = FV_Atm(1)%bd%isc ; iec = FV_Atm(1)%bd%iec
+  jsc = FV_Atm(1)%bd%jsc ; jec = FV_Atm(1)%bd%jec
+
+  call mpp_update_domains(FV_Atm(1)%pt, FV_Atm(1)%domain, complete=.true.)
+
+  !$OMP PARALLEL DO DEFAULT(NONE) &
+  !$OMP SHARED(FV_Atm, mag_grad_t, isc, iec, jsc, jec) &
+  !$OMP PRIVATE(i, j, k, dx, dy, dTdx, dTdy)
+  do k = 1, FV_Atm(1)%npz
+     do j = jsc, jec
+        do i = isc, iec
+
+           ! Note: Replace dxa/dya with your exact FV3 metric arrays for cell-center distances.
+           ! (Often found in FV_Atm(1)%gridstruct%dxa or similar)
+           dx = FV_Atm(1)%gridstruct%dxa(i,j)
+           dy = FV_Atm(1)%gridstruct%dya(i,j)
+
+           ! Assuming pt is Absolute Temperature (if it's Potential T, multiply by pkz here)
+           dTdx = (FV_Atm(1)%pt(i+1,j,k) - FV_Atm(1)%pt(i-1,j,k)) / (2.0_REAL8 * dx)
+           dTdy = (FV_Atm(1)%pt(i,j+1,k) - FV_Atm(1)%pt(i,j-1,k)) / (2.0_REAL8 * dy)
+
+           ! Magnitude in K/m
+           mag_grad_t(i,j,k) = SQRT(dTdx**2 + dTdy**2)
+
+        enddo
+     enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+return
+end subroutine fv_getGradT_3D
+
 subroutine fv_getDELZ(delz,temp,pe)
   real(REAL8), intent(OUT) :: delz(FV_Atm(1)%bd%isc:FV_Atm(1)%bd%iec,FV_Atm(1)%bd%jsc:FV_Atm(1)%bd%jec,1:FV_Atm(1)%npz)
   real(REAL8), intent( IN) :: temp(FV_Atm(1)%bd%isc:FV_Atm(1)%bd%iec,FV_Atm(1)%bd%jsc:FV_Atm(1)%bd%jec,1:FV_Atm(1)%npz)