Adds the matsuno-gill heating pattern as an idealized forcing namelist option.

OkayHughes · OkayHughes · commit 2582627ae457 · 2025-01-14T14:03:00.000-05:00
It is enabled in standalone HOMME by using the namelist option `test_case = "matsuno_gill"`. This test case is described in, e.g., (this paper)[https://iopscience.iop.org/article/10.1016/j.fluiddyn.2004.03.003/pdf). Add the vertical coordinate files for the Matsuno-Gill heating test. The hybrid coefficients are analytically calculated to give constant ∆z between model levels for this test case, with a model top at 10km. Update CMakeLists.txt for theta-l_kokkos to include the new file for the matsuno gill heating profile.
diff --git a/components/homme/src/preqx_kokkos/CMakeLists.txt b/components/homme/src/preqx_kokkos/CMakeLists.txt
@@ -127,6 +127,7 @@ MACRO(PREQX_KOKKOS_SETUP)
     ${TEST_SRC_DIR}/dcmip2016-terminator.F90
     ${TEST_SRC_DIR}/dcmip2016-tropical-cyclone.F90
     ${TEST_SRC_DIR}/held_suarez_mod.F90
+    ${TEST_SRC_DIR}/matsuno_gill_mod.F90
     ${TEST_SRC_DIR}/dry_planar.F90
     ${TEST_SRC_DIR}/moist_planar.F90
     ${TEST_SRC_DIR}/mtests.F90
diff --git a/components/homme/src/share/namelist_mod.F90 b/components/homme/src/share/namelist_mod.F90
@@ -555,6 +555,7 @@ subroutine readnl(par)
 #ifdef _PRIM
        write(iulog,*) "reading physics namelist..."
        if (test_case      == "held_suarez0"   .or. &
+           test_case      == "matsuno_gill"   .or. &
            test_case(1:10)== "baroclinic"     .or. &
            test_case(1:13)== "jw_baroclinic"  .or. &
            test_case(1:5) == "dcmip"          .or. &
diff --git a/components/homme/src/test_mod.F90 b/components/homme/src/test_mod.F90
@@ -26,6 +26,7 @@ module test_mod
                                 dcmip2016_test1_forcing, dcmip2016_test2_forcing, dcmip2016_test3_forcing, &
                                 dcmip2016_pg_init, dcmip2016_test1_pg, dcmip2016_test1_pg_forcing, dcmip2016_init
 use held_suarez_mod,      only: hs0_init_state
+use matsuno_gill_mod,      only: mg_init_state
 
 use dry_planar_tests,     only: planar_hydro_gravity_wave_init, planar_nonhydro_gravity_wave_init
 use dry_planar_tests,     only: planar_hydro_mountain_wave_init, planar_nonhydro_mountain_wave_init, planar_schar_mountain_wave_init
@@ -84,6 +85,7 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n
     case('mtest1'); test_with_forcing = .true. ;
     case('mtest2'); test_with_forcing = .true. ;
     case('mtest3'); test_with_forcing = .true. ;
+    case('matsuno_gill'); test_with_forcing = .true. ;
     case('held_suarez0'); test_with_forcing = .true. ;
     case('jw_baroclinic');
     case('planar_hydro_gravity_wave');
@@ -139,6 +141,7 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n
       case('mtest2');             call mtest_init       (elem,hybrid,hvcoord,nets,nete,2)
       case('mtest3');             call mtest_init       (elem,hybrid,hvcoord,nets,nete,3)
       case('held_suarez0');       call hs0_init_state   (elem,hybrid,hvcoord,nets,nete,300.0_rl)
+      case('matsuno_gill');    call mg_init_state  (elem,hybrid,hvcoord,nets,nete)
       case('jw_baroclinic');      call jw_baroclinic    (elem,hybrid,hvcoord,nets,nete)
       case('planar_hydro_gravity_wave');            call planar_hydro_gravity_wave_init(elem,hybrid,hvcoord,nets,nete)
       case('planar_nonhydro_gravity_wave');         call planar_nonhydro_gravity_wave_init(elem,hybrid,hvcoord,nets,nete)
@@ -160,7 +163,6 @@ subroutine set_test_initial_conditions(elem, deriv, hybrid, hvcoord, tl, nets, n
 
     endselect
 !  endif
-
     if (midpoint_eta_dot_dpdn) then
        do ie = nets,nete
           elem(ie)%derived%eta_dot_dpdn = 0
@@ -211,6 +213,7 @@ subroutine compute_test_forcing(elem,hybrid,hvcoord,nt,ntQ,dt,nets,nete,tl)
 
   use dcmip12_wrapper, only:  dcmip2012_test2_x_forcing
   use held_suarez_mod, only: hs_forcing
+  use matsuno_gill_mod, only: mg_forcing
   use control_mod,     only: ftype
   implicit none
   type(element_t),  intent(inout) :: elem(:)                            ! element array
@@ -261,6 +264,11 @@ subroutine compute_test_forcing(elem,hybrid,hvcoord,nt,ntQ,dt,nets,nete,tl)
        do ie=nets,nete
           call hs_forcing(elem(ie),hvcoord,nt,ntQ,dt)
        enddo
+    case('matsuno_gill');
+       do ie=nets,nete
+          call mg_forcing(elem(ie),hvcoord,nt,ntQ,dt)
+       enddo
+
 
   endselect
 
diff --git a/components/homme/src/test_src/matsuno_gill_mod.F90 b/components/homme/src/test_src/matsuno_gill_mod.F90
@@ -0,0 +1,297 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+module matsuno_gill_mod
+
+
+  use coordinate_systems_mod, only: spherical_polar_t
+  use dimensions_mod,         only: nlev,np,qsize,nlevp
+  use element_mod,            only: element_t
+  use element_state,          only: timelevels
+  use element_ops,            only: set_thermostate, get_temperature
+  use hybrid_mod,             only: hybrid_t
+  use hybvcoord_mod,          only: hvcoord_t
+  use kinds,                  only: real_kind, iulog
+  use physical_constants,     only: p0, kappa,g, dd_pi, Rgas,Cp
+  use physics_mod,            only: prim_condense
+  use time_mod,               only: secpday
+#ifndef HOMME_WITHOUT_PIOLIBRARY
+  use common_io_mod,          only: infilenames
+#endif
+
+implicit none
+private
+
+  real (kind=real_kind), public, parameter :: sigma_b = 0.70D0
+  real (kind=real_kind), public, parameter :: k_Tv     = 1.0D0/(40.0D0*secpday)
+  real (kind=real_kind), public, parameter :: bv_sq   = 1.0D-4
+  real (kind=real_kind), public, parameter :: T0 = 300
+  real (kind=real_kind), public, parameter :: q0   = 5d0 /secpday
+  real (kind=real_kind), public, parameter :: dlam   = 1.0D0 * 30D0 * dd_pi/180D0 
+  real (kind=real_kind), public, parameter :: dphi   = 1.0D0 * 10D0 * dd_pi/180D0 
+
+  public :: mg_v_forcing
+  public :: mg_T_forcing
+  public :: mg_init_state
+  public :: mg_forcing
+
+contains
+
+  subroutine mg_forcing(elemin,hvcoord,nm1,nm1_Q,dt)
+
+    type (element_t)                 :: elemin
+    type (hvcoord_t)                 :: hvcoord
+    integer                          :: nm1,nm1_Q  ! timelevel to use
+    real (kind=real_kind)                   :: dt
+
+    ! local
+    real (kind=real_kind)                   :: pmid,r0,r1,dtf_q,dp,rdp,FQ
+    real (kind=real_kind), dimension(np,np) :: psfrc 
+    real (kind=real_kind)                   :: temperature(np,np,nlev)
+    real (kind=real_kind)                   :: v(np,np,3,nlev)
+    real (kind=real_kind)                   :: fv(np,np,3,nlev)
+    integer                                 :: i,j,k,q
+
+    dtf_q = dt
+    call get_temperature(elemin,temperature,hvcoord,nm1)
+        
+    do j=1,np
+       do i=1,np
+          psfrc(i,j) = (elemin%state%ps_v(i,j,nm1))
+       end do
+    end do
+
+    elemin%derived%FT(:,:,:) = elemin%derived%FT(:,:,:) + &
+         mg_T_forcing(hvcoord,psfrc(1,1),               &
+         temperature,elemin%spherep,np, nlev)
+
+    v(:,:,1:2,:) = elemin%state%v(:,:,1:2,:,nm1)
+#if ( defined MODEL_THETA_L ) 
+    v(:,:,3,:) = elemin%state%w_i(:,:,1:nlev,nm1)  ! dont apply at surface
+#else
+    v(:,:,3,:) = 0
+#endif
+
+    fv = mg_v_forcing(hvcoord,psfrc(1,1),v,np,nlev)
+
+#if ( defined MODEL_THETA_L ) 
+    elemin%derived%FM(:,:,1:3,:) = elemin%derived%FM(:,:,1:3,:) + fv(:,:,1:3,:)
+#else
+    elemin%derived%FM(:,:,1:2,:) = elemin%derived%FM(:,:,1:2,:) + fv(:,:,1:2,:)
+#endif
+
+    if (qsize>=1) then
+       ! HS with tracer  (Galewsky type forcing, with flux of  2.3e-5 kg/m^2/s
+       ! MASS in kg/m^2   = < Q dp_in_Pa / g >   
+       ! flux in kg/m^2/s = < FQ dp_in_Pa / g >   
+       ! We want < FQ dp_in_Pa / g > = 2.3e-5  so:  FQ = 2.3e-5*g/dp_in_Pa 
+
+       ! lowest layer thickness, in Pa
+       dp = ( hvcoord%hyai(nlev+1) - hvcoord%hyai(nlev) ) + &
+               ( hvcoord%hybi(nlev+1) - hvcoord%hybi(nlev) )*1000*100
+       rdp = 1./ dp
+       q=1
+       do j=1,np
+          do i=1,np
+             FQ = rdp * g * 2.3E-5 * COS(elemin%spherep(i,j)%lat)**2
+             elemin%derived%FQ(i,j,nlev,q) =elemin%derived%FQ(i,j,nlev,q)+FQ
+          enddo
+       enddo
+
+       do j=1,np
+          do i=1,np
+             do k=1,nlev
+                pmid = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*(elemin%state%ps_v(i,j,nm1))
+                r0=elemin%state%Q(i,j,k,q)
+                r1=r0
+                call Prim_Condense(r1,temperature(i,j,k),pmid)
+                elemin%derived%FQ(i,j,k,q) = elemin%derived%FQ(i,j,k,q) + &
+                     (r1-r0)/(dtf_q)
+             enddo
+          enddo
+       enddo
+    endif
+
+  end subroutine mg_forcing
+  function z_given_p(p) result(z)
+      real (kind=real_kind), intent(in) :: p
+      real (kind=real_kind):: z
+      z = (-g/bv_sq)*log( (T0/(g**2/(bv_sq*Cp)))*( (p/p0)**(Rgas/Cp) - 1.0  ) + 1.0 )
+      !z = (-g/Nsq)*np.log( (T0/bigG)*( (p/p0)**(Rd/Cp) - 1.0  ) + 1.0 )
+  end function 
+  function T_given_p(p) result(T)
+      real (kind=real_kind), intent(in) :: p
+      real (kind=real_kind) :: T
+      T = T0 * (p/p0)**kappa * exp(bv_sq/g * z_given_p(p))
+  end function
+ 
+  function mg_v_forcing(hvcoord,ps,v,npts,nlevels) result(mg_v_frc)
+
+    integer, intent(in)               :: npts
+    integer, intent(in)               :: nlevels
+    type (hvcoord_t), intent(in)       :: hvcoord
+    real (kind=real_kind), intent(in) :: ps(npts,npts)
+
+    real (kind=real_kind), intent(in) :: v(npts,npts,3,nlevels)
+    real (kind=real_kind)             :: mg_v_frc(npts,npts,3,nlevels)
+
+    ! Local variables
+
+    integer i,j,k
+    real (kind=real_kind) :: k_v
+    real (kind=real_kind) :: p,etam
+
+    do k=1,nlevels
+       do j=1,npts
+          do i=1,npts
+             p    = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*ps(i,j)
+             mg_v_frc(i,j,1,k) = -k_Tv*v(i,j,1,k)
+             mg_v_frc(i,j,2,k) = -k_Tv*v(i,j,2,k)
+
+             mg_v_frc(i,j,3,k) = -k_Tv*v(i,j,3,k)
+          end do
+       end do
+    end do
+
+  end function mg_v_forcing
+
+  function mg_T_forcing(hvcoord,ps,T,sphere,npts,nlevels) result(mg_T_frc)
+
+    integer, intent(in) :: npts
+    integer, intent(in) :: nlevels
+
+    type (hvcoord_t), intent(in)          :: hvcoord
+    real (kind=real_kind), intent(in)    :: ps(npts,npts)
+    real (kind=real_kind), intent(in)    :: T(npts,npts,nlevels)
+    type (spherical_polar_t), intent(in) :: sphere(npts,npts)
+
+    real (kind=real_kind)             :: mg_T_frc(npts,npts,nlevels)
+
+    ! Local variables
+
+    real (kind=real_kind) :: p,logprat,pratk,Teq,T_forcing
+    real (kind=real_kind) :: logps0,etam
+    real (kind=real_kind) :: lat,snlat
+
+    real (kind=real_kind) :: k_t(npts,npts)
+    real (kind=real_kind) :: snlatsq(npts,npts)
+    real (kind=real_kind) :: latlon_comp(npts,npts), z_top(npts,npts)
+    real (kind=real_kind) :: cslatsq(npts,npts)
+
+    real (kind=real_kind) :: rec_one_minus_sigma_b
+
+    integer i,j,k
+
+    logps0    = LOG(hvcoord%ps0 )
+
+    do j=1,npts
+       do i=1,npts
+         snlat        = SIN(sphere(i,j)%lat)
+         if (abs(sphere(i,j)%lon-dd_pi/2) < dlam .and. abs(sphere(i,j)%lat) < dphi ) then
+         latlon_comp(i,j) = cos((sphere(i,j)%lon-dd_pi/2)/(2 * dlam))**2 * cos(sphere(i,j)%lat/(2*dphi))**2
+         else
+           latlon_comp(i,j) = 0
+         end if
+         z_top(i,j) = z_given_p(hvcoord%hyai(1)*hvcoord%ps0 + hvcoord%hybi(1)*ps(i,j))
+       end do
+    end do
+
+    rec_one_minus_sigma_b = 1.0D0/(1.0D0 - sigma_b)
+
+    do k=1,nlevels
+       do j=1,npts
+          do i=1,npts
+             p         = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*ps(i,j)
+
+             T_forcing = q0 * latlon_comp(i,j) * sin(dd_pi * z_given_p(p)/z_top(i,j))
+#if 0
+             ! ======================================
+             ! This is a smooth forcing 
+             ! for debugging purposes only...
+             ! ======================================
+
+#endif
+             mg_T_frc(i,j,k)= -k_Tv*(T(i,j,k)-T_given_p(p)) + T_forcing
+          end do
+       end do
+    end do
+     
+  end function mg_T_forcing
+
+  subroutine mg_init_state(elem, hybrid, hvcoord,nets,nete)
+
+    type(element_t),        intent(inout) :: elem(:)
+    type(hybrid_t),         intent(in)    :: hybrid                   ! hybrid parallel structure
+    type (hvcoord_t),       intent(in)    :: hvcoord
+    integer,                intent(in)    :: nets
+    integer,                intent(in)    :: nete
+
+    ! Local variables
+    
+    integer ie,i,j,k,q,tl
+    integer :: nm1 
+    integer :: n0 
+    integer :: np1
+    real (kind=real_kind) :: lat_mtn,lon_mtn,r_mtn,h_mtn,rsq,lat,lon,p_tmp
+    real (kind=real_kind) :: temperature(np,np,nlev),p(np,np),exner(np,np),ps(np,np)
+
+    if (hybrid%masterthread) write(iulog,*) 'initializing Held-Suarez primitive equations test'
+
+    nm1= 1
+    n0 = 2
+    np1= 3
+
+    do ie=nets,nete
+
+       elem(ie)%state%ps_v(:,:,n0) =p0
+       elem(ie)%state%ps_v(:,:,nm1)=p0
+       elem(ie)%state%ps_v(:,:,np1)=p0
+
+       elem(ie)%state%v(:,:,:,:,n0) =0.0D0
+       elem(ie)%state%v(:,:,:,:,nm1)=elem(ie)%state%v(:,:,:,:,n0)
+       elem(ie)%state%v(:,:,:,:,np1)=elem(ie)%state%v(:,:,:,:,n0)
+
+#ifdef MODEL_THETA_L
+       elem(ie)%state%w_i = 0.0
+#endif
+
+
+       ! if topo file was given in the namelist, PHIS was initilized in prim_main
+       ! otherwise assume 0
+#ifndef HOMME_WITHOUT_PIOLIBRARY
+       if (infilenames(1)=='') then
+          elem(ie)%state%phis(:,:)=0.0D0
+       endif
+#endif
+
+       do k=1,nlev
+         p_tmp = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem(ie)%state%ps_v(1,1,n0)
+         temperature(1, 1, k) = T_given_p(p_tmp)
+       end do
+       do i=1,np
+         do j=1,np
+           temperature(i,j, :) = temperature(1,1, :)
+         end do
+       end do
+
+
+
+
+       !if (qsize>=1) then
+       !q=1
+       !elem(ie)%state%Q(:,:,:,q) =0  ! moist HS tracer IC=0
+       !do q=2,qsize
+       !   elem(ie)%state%Q(:,:,:,q) =temperature(:,:,:)/400
+       !enddo
+       !endif
+       ps=elem(ie)%state%ps_v(:,:,n0)
+       call set_thermostate(elem(ie),ps,temperature,hvcoord) ! sets z based on hydrostatic presumption
+    end do
+
+
+  end subroutine mg_init_state
+
+
+end module matsuno_gill_mod
+
diff --git a/components/homme/src/theta-l_kokkos/CMakeLists.txt b/components/homme/src/theta-l_kokkos/CMakeLists.txt
@@ -129,6 +129,7 @@ MACRO(THETAL_KOKKOS_SETUP)
     ${TEST_SRC_DIR}/dcmip2016-terminator.F90
     ${TEST_SRC_DIR}/dcmip2016-tropical-cyclone.F90
     ${TEST_SRC_DIR}/held_suarez_mod.F90
+    ${TEST_SRC_DIR}/matsuno_gill_mod.F90
     ${TEST_SRC_DIR}/dry_planar.F90
     ${TEST_SRC_DIR}/moist_planar.F90
     ${TEST_SRC_DIR}/mtests.F90
diff --git a/components/homme/test/vcoord/mg-20i.ascii b/components/homme/test/vcoord/mg-20i.ascii
@@ -0,0 +1,13 @@
+21 !hyai
+0.27381905404907 0.26590539035976 0.25752394878279
+0.24865429808716 0.23927536347865 0.22936541085572 0.21890203071294
+0.20786212168493 0.19622187372385 0.18395675090318 0.17104147384052
+0.15745000173179 0.14315551398919 0.12813039147516 0.11234619732416
+0.095773657344247 0.078382639990006 0.060142135898353
+0.041020236978492 0.020984115047143 0
+21 ! hybi
+0 0.028901070149364 0.059510487036309 0.091902866472543
+0.1261551745929 0.16234678535331 0.20055953931639 0.24087780374962
+0.28338853406205 0.32818133660555 0.37534853286773 0.42498522508382
+0.4771893632956 0.53206181388604 0.58970642961892 0.65023012121324
+0.71374293048299 0.78035810507338 0.85019217482527 0.92336502980036 1
diff --git a/components/homme/test/vcoord/mg-20m.ascii b/components/homme/test/vcoord/mg-20m.ascii
