Merge branch 'ambrad:ambrad/eam/physgrid-ho-pp-remap' into master (PR #…

…3141)

Adds high-order, property-preserving remap between physics and GLL grids.

Provide a module, gllfvremap_mod, whose use in dynamics/se is optionally enabled
with atmosphere namelist option se_fv_phys_remap_alg=1 (default 0 uses the
original remap algorithm), to remap between physics and GLL grids. The remap

-- is mass conserving;
-- produces tracers that do not violate extrema established before the remap;
-- remaps tendencies only, so that the background state is not affected by remap;
-- has order of accuracy N, where N is the parameter in the physgrid specification pgN.

[BFB]

cime/config/e3sm/tests.py

@@ -43,6 +43,7 @@
             "SMS_R_Ld5.ne4_ne4.FSCM5A97",
             "SMS_D_Ln5.ne4_ne4.FC5AV1C-L",
             "SMS_Ln5.ne4pg2_ne4pg2.FC5AV1C-L",
+            "SMS_Ln5.ne4pg2_ne4pg2.FC5AV1C-L.cam-thetahy_pg2",
             )
         },
 

components/cam/bld/namelist_files/namelist_defaults_cam.xml

@@ -1553,6 +1553,8 @@
 <semi_lagrange_cdr_check> .false. </semi_lagrange_cdr_check>
 <semi_lagrange_nearest_point_lev> 0 </semi_lagrange_nearest_point_lev>
 
+<se_fv_phys_remap_alg> 0 </se_fv_phys_remap_alg>
+
 <!-- ================================================================== -->
 <!-- Defaults for driver namelist seq_infodata_inparm                   -->
 <!-- ================================================================== -->

components/cam/bld/namelist_files/namelist_definition.xml

@@ -6065,8 +6065,22 @@ nearest point in the halo to the departure point, when necessary.
 Default: 0
 </entry>
 
+<entry id="semi_lagrange_hv_q" type="integer" category="se"
+       group="ctl_nl" valid_values="">
+Number of tracers, starting from 1, to which to apply hyperviscosity. For
+example, to apply hyperviscosity to moisture, the first tracer, set the value to
+1. Hyperviscosity is applied using nu_q as the coefficient.
+</entry>
+
 <!-- Physics grid -->
 
+<entry id="se_fv_phys_remap_alg" type="integer" category="se"
+       group="ctl_nl" valid_values="0,1">
+Remap algorithm to use to map between FV physics grid and GLL dynamics
+grid. 0, default, provides the original algorithms; 1 provides
+high-order algorithms implemented in gllfvremap_mod.
+</entry>
+
 <entry id="se_fv_nphys" type="integer" category="se"
        group="dyn_se_inparm" valid_values="" >
 Number of equally-spaced horizontal physics points per spectral
@@ -6076,13 +6090,6 @@ result in 4 equally-spaced physics points per element).
 Default: 0 = feature disabled, use dynamics GLL points.
 </entry>
 
-<entry id="semi_lagrange_hv_q" type="integer" category="se"
-       group="ctl_nl" valid_values="">
-Number of tracers, starting from 1, to which to apply hyperviscosity. For
-example, to apply hyperviscosity to moisture, the first tracer, set the value to
-1. Hyperviscosity is applied using nu_q as the coefficient.
-</entry>
-
 <!-- CAM I/O  -->
 
 <entry id="pio_stride" type="integer" category="pio"

components/cam/cime_config/testdefs/testmods_dirs/cam/thetahy_pg2/shell_commands

@@ -0,0 +1,2 @@
+#!/bin/bash
+./xmlchange CAM_TARGET=theta-l

components/cam/cime_config/testdefs/testmods_dirs/cam/thetahy_pg2/user_nl_cam

@@ -0,0 +1,2 @@
+theta_hydrostatic_mode=.true.
+se_fv_phys_remap_alg=1

components/cam/src/dynamics/se/dp_coupling.F90

@@ -40,6 +40,9 @@ subroutine d_p_coupling(phys_state, phys_tend,  pbuf2d, dyn_out)
     use dyn_comp,                only: frontgf_idx, frontga_idx, hvcoord
     use phys_control,            only: use_gw_front
     use fv_physics_coupling_mod, only: dyn_to_fv_phys
+    use control_mod,             only: se_fv_phys_remap_alg
+    use dyn_comp,                only: dom_mt
+    use gllfvremap_mod,          only: gfr_dyn_to_fv_phys
 
     implicit none
     !---------------------------------------------------------------------------
@@ -113,9 +116,14 @@ subroutine d_p_coupling(phys_state, phys_tend,  pbuf2d, dyn_out)
         ! Map dynamics state to FV physics grid
         !-----------------------------------------------------------------------
         call t_startf('dyn_to_fv_phys')
-        call dyn_to_fv_phys(elem,ps_tmp(1:nphys_sq,:),zs_tmp(1:nphys_sq,:),     &
-                                 T_tmp(1:nphys_sq,:,:),uv_tmp(1:nphys_sq,:,:,:),&
-                                 om_tmp(1:nphys_sq,:,:),q_tmp(1:nphys_sq,:,:,:))
+        if (se_fv_phys_remap_alg == 0) then
+          call dyn_to_fv_phys(elem,ps_tmp(1:nphys_sq,:),zs_tmp(1:nphys_sq,:),     &
+               T_tmp(1:nphys_sq,:,:),uv_tmp(1:nphys_sq,:,:,:),&
+               om_tmp(1:nphys_sq,:,:),q_tmp(1:nphys_sq,:,:,:))
+        else
+          call gfr_dyn_to_fv_phys(par, dom_mt, tl_f, hvcoord, elem, ps_tmp, zs_tmp, &
+               T_tmp, uv_tmp, om_tmp, q_tmp)
+        end if
         call t_stopf('dyn_to_fv_phys')
 
         !-----------------------------------------------------------------------
@@ -334,7 +342,10 @@ subroutine p_d_coupling(phys_state, phys_tend,  dyn_in)
     use shr_vmath_mod,           only: shr_vmath_log
     use cam_control_mod,         only: adiabatic
     use fv_physics_coupling_mod, only: fv_phys_to_dyn
-    use control_mod,             only: ftype
+    use control_mod,             only: ftype, se_fv_phys_remap_alg
+    use dyn_comp,                only: dom_mt, hvcoord
+    use gllfvremap_mod,          only: gfr_fv_phys_to_dyn
+    use time_manager,            only: get_step_size
     implicit none
     ! INPUT PARAMETERS:
     type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
@@ -349,6 +360,7 @@ subroutine p_d_coupling(phys_state, phys_tend,  dyn_in)
     real (kind=real_kind), dimension(npsq,pver,nelemd)       :: T_tmp  ! temp array to hold T
     real (kind=real_kind), dimension(npsq,2,pver,nelemd)     :: uv_tmp ! temp array to hold u and v
     real (kind=real_kind), dimension(npsq,pver,pcnst,nelemd) :: q_tmp  ! temp to hold advected constituents
+    real (kind=real_kind)    :: dtime
     integer(kind=int_kind)   :: m, i, j, k                 ! loop iterators
     integer(kind=int_kind)   :: gi(2), gj(2)               ! index list used to simplify pg2 case
     integer(kind=int_kind)   :: di, dj
@@ -466,11 +478,18 @@ subroutine p_d_coupling(phys_state, phys_tend,  dyn_in)
 
     if (par%dynproc) then
       if (fv_nphys > 0) then
-
+        call t_startf('fv_phys_to_dyn')
         ! Map FV physics state to dynamics grid
-        call fv_phys_to_dyn(elem,T_tmp(1:nphys_sq,:,:),   &
-                                 uv_tmp(1:nphys_sq,:,:,:),&
-                                 q_tmp(1:nphys_sq,:,:,:))
+        if (se_fv_phys_remap_alg == 0) then
+          call fv_phys_to_dyn(elem,T_tmp(1:nphys_sq,:,:),   &
+               uv_tmp(1:nphys_sq,:,:,:),&
+               q_tmp(1:nphys_sq,:,:,:))
+        else
+          dtime = get_step_size()
+          call gfr_fv_phys_to_dyn(par, dom_mt, TimeLevel%n0, dtime, hvcoord, elem, T_tmp, &
+               uv_tmp, q_tmp)
+        end if
+        call t_stopf('fv_phys_to_dyn')
 
       else ! physics is on GLL nodes
 

components/cam/src/dynamics/se/dyn_grid.F90

@@ -1506,12 +1506,13 @@ end subroutine dyn_grid_get_elem_coords
   !=================================================================================================
   !
   subroutine fv_physgrid_init()
-    use control_mod,            only: cubed_sphere_map
+    use control_mod,            only: cubed_sphere_map, se_fv_phys_remap_alg
     use kinds,                  only: lng_dbl => longdouble_kind
     use cube_mod,               only: ref2sphere
     use coordinate_systems_mod, only: spherical_polar_t, cartesian3D_t
     use coordinate_systems_mod, only: sphere_tri_area, change_coordinates 
     use derivative_mod,         only: allocate_subcell_integration_matrix
+    use gllfvremap_mod,         only: gfr_init
     !------------------------------Arguments------------------------------------
     ! type(element_t)         , intent(in   ) :: elem(:)
     ! type(fv_physgrid_struct), intent(inout) :: fv_physgrid(nelemd)
@@ -1625,11 +1626,15 @@ subroutine fv_physgrid_init()
       end do ! j
     end do ! ie
 
+    if (se_fv_phys_remap_alg == 1) call gfr_init(par, elem, fv_nphys)
   end subroutine fv_physgrid_init
   !
   !=================================================================================================
   !
   subroutine fv_physgrid_final()
+    use control_mod,      only: se_fv_phys_remap_alg
+    use gllfvremap_mod,   only: gfr_finish
+
     !----------------------------Local-Variables--------------------------------
     integer :: ie
     !---------------------------------------------------------------------------  
@@ -1640,6 +1645,8 @@ subroutine fv_physgrid_final()
       deallocate( fv_physgrid(ie)%corner_lat )
       deallocate( fv_physgrid(ie)%corner_lon )
     end do ! ie
+
+    if (se_fv_phys_remap_alg == 1) call gfr_finish()
   end subroutine fv_physgrid_final
   !
   !=================================================================================================

components/cam/src/dynamics/se/inidat.F90

@@ -28,7 +28,7 @@ module inidat
 contains
 
   subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)
-    use dyn_comp,                only: dyn_import_t, hvcoord
+    use dyn_comp,                only: dyn_import_t, hvcoord, dom_mt
     use parallel_mod,            only: par
     use bndry_mod,               only: bndry_exchangev
     use constituents,            only: cnst_name, cnst_read_iv, qmin
@@ -60,6 +60,8 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)
     use se_single_column_mod,    only: scm_setinitial
     use element_ops,             only: set_thermostate
     use fv_physics_coupling_mod, only: fv_phys_to_dyn_topo
+    use control_mod,             only: se_fv_phys_remap_alg
+    use gllfvremap_mod,          only: gfr_fv_phys_to_dyn_topo
     implicit none
     type(file_desc_t),intent(inout) :: ncid_ini, ncid_topo
     type (dyn_import_t), target, intent(inout) :: dyn_in   ! dynamics import
@@ -455,10 +457,14 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)
       fieldname = 'PHIS'
       tmp(:,1,:) = 0.0_r8
       if (fv_nphys > 0) then
-         call infld(fieldname, ncid_topo, 'ncol', 1, nphys_sq, &
-                    1, nelemd, phys_tmp, found, gridname='physgrid_d')
          ! Copy phis field to GLL grid
-         call fv_phys_to_dyn_topo(elem,phys_tmp)
+         call infld(fieldname, ncid_topo, 'ncol', 1, nphys_sq, &
+              1, nelemd, phys_tmp, found, gridname='physgrid_d')
+         if (se_fv_phys_remap_alg == 0) then
+            call fv_phys_to_dyn_topo(elem,phys_tmp)
+         else
+            call gfr_fv_phys_to_dyn_topo(par, dom_mt, elem, phys_tmp)
+         end if
       else
          call infld(fieldname, ncid_topo, ncol_name,      &
             1, npsq, 1, nelemd, tmp(:,1,:), found, gridname=grid_name)

components/homme/src/share/compose_mod.F90

@@ -182,6 +182,7 @@ subroutine compose_init(par, elem, GridVertex)
     logical(kind=c_bool) :: use_sgi, owned
     integer, allocatable :: owned_ids(:)
     integer, pointer :: rank2sfc(:) => null()
+    integer, target :: null_target(1)
 
 #ifdef HOMME_ENABLE_COMPOSE
     call t_startf('compose_init')
@@ -207,6 +208,12 @@ subroutine compose_init(par, elem, GridVertex)
              sc2rank(sc) = GridVertex(i)%processor_number - 1
           end do
        end if
+    else
+       ! These lines fix ifort -check catches. The data are not used,
+       ! but the function call cedr_init_impl makes -check think they
+       ! are used.
+       allocate(owned_ids(1))
+       rank2sfc => null_target
     end if
     if (use_sgi) then
        call cedr_init_impl(par%comm, semi_lagrange_cdr_alg, &

components/homme/src/share/compose_test_mod.F90

@@ -115,7 +115,7 @@ subroutine compose_test(par, hvcoord, dom_mt, elem)
 
     ! 2. Standalone tracer advection test, useful as a basic but comprehensive
     ! correctness test and also as part of a convergence test.
-    call compose_stt(hybrid, nets, nete, hvcoord, deriv, elem)
+    call compose_stt(hybrid, dom_mt, nets, nete, hvcoord, deriv, elem)
 #if (defined HORIZ_OPENMP)
     !$omp end parallel
 #endif
@@ -153,9 +153,10 @@ subroutine print_software_statistics(hybrid, nets, nete)
     end if
   end subroutine print_software_statistics
 
-  subroutine compose_stt(hybrid, nets, nete, hvcoord, deriv, elem)
+  subroutine compose_stt(hybrid, dom_mt, nets, nete, hvcoord, deriv, elem)
     use iso_c_binding, only: c_loc
     use parallel_mod, only: parallel_t
+    use domain_mod, only: domain1d_t
     use hybrid_mod, only: hybrid_t
     use element_mod, only: element_t
     use time_mod, only: timelevel_t, timelevel_init_default, timelevel_qdp
@@ -164,13 +165,16 @@ subroutine compose_stt(hybrid, nets, nete, hvcoord, deriv, elem)
     use derivative_mod, only: derivative_t, derivinit
     use dimensions_mod, only: ne, np, nlev, qsize, qsize_d, nelemd
     use coordinate_systems_mod, only: spherical_polar_t
-    use control_mod, only: qsplit, statefreq
+    use control_mod, only: qsplit, statefreq, se_fv_phys_remap_alg
     use time_mod, only: nmax
     use hybvcoord_mod, only: hvcoord_t
     use perf_mod
     use sl_advection
+    use gllfvremap_mod
+    use gllfvremap_test_mod
 
     type (hybrid_t), intent(in) :: hybrid
+    type (domain1d_t), pointer, intent(in) :: dom_mt(:)
     type (derivative_t), intent(in) :: deriv
     type (element_t), intent(inout) :: elem(:)
     type (hvcoord_t) , intent(in) :: hvcoord
@@ -182,6 +186,12 @@ subroutine compose_stt(hybrid, nets, nete, hvcoord, deriv, elem)
     integer :: nsteps, n0_qdp, np1_qdp, ie, i, j
     real (kind=real_kind) :: dt, tprev, t
 
+    if (se_fv_phys_remap_alg == 1) then
+       call gfr_test(hybrid, dom_mt, hvcoord, deriv, elem)
+       call gfr_check_api(hybrid, nets, nete, hvcoord, elem)
+       return
+    end if
+
 #ifdef HOMME_ENABLE_COMPOSE  
     call t_startf('compose_stt')
     ! Set up time stepping and initialize q and density.

components/homme/src/share/control_mod.F90

@@ -13,7 +13,6 @@ module control_mod
   integer, public, parameter :: MAX_FILE_LEN=240
   character(len=MAX_STRING_LEN)    , public :: integration    ! time integration (explicit, or full imp)
 
-! experimental option for preqx model:
   ! Tracer transport algorithm type:
   !     0  spectral-element Eulerian
   !    12 interpolation semi-Lagrangian
@@ -182,6 +181,9 @@ module control_mod
 
   logical, public :: disable_diagnostics  = .FALSE.
 
+  ! Physgrid parameters
+  integer, public :: se_fv_phys_remap_alg = 0
+
 
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!