Fixes for incorrect coupling field data (ESCOMP#10)

Commit summary:
* 71baab0 Fixed wrong domain decomposition by switching to ORANGE partitioning scheme
* 6a0fce3 Fixed syntax for checking missing dict keys
* 1bea6fd Changed to bundled coupling fields (w/ extra dim = # of soil levels)
* fb6b252 Always compile with OpenMP.
* 906ff7b Fixed mismatched coupling timestep with ParFlow

cmake/SetBuildOptions.cmake

@@ -23,18 +23,18 @@ endif()
 # Set compiler specific flags.
 if(COMPILER STREQUAL "GNU")
     add_compile_definitions(CPRGNU)
-    set(CMAKE_C_FLAGS "-std=gnu99")
+    set(CMAKE_C_FLAGS "-std=gnu99 -fopenmp")
     set(CMAKE_C_FLAGS_DEBUG "-fcheck=bounds")
     set(CMAKE_C_FLAGS_RELEASE "-O")
-    set(CMAKE_Fortran_FLAGS "-fconvert=big-endian -ffree-line-length-none -ffixed-line-length-none -ffree-form")
+    set(CMAKE_Fortran_FLAGS "-fconvert=big-endian -ffree-line-length-none -ffixed-line-length-none -ffree-form -fopenmp")
     set(CMAKE_Fortran_FLAGS_DEBUG "-g -Wall -Og -fbacktrace -ffpe-trap=zero,overflow -fcheck=bounds")
     set(CMAKE_Fortran_FLAGS_RELEASE "-O") 
 elseif(COMPILER STREQUAL "Intel")
     add_compile_definitions(CPRINTEL)
-    set(CMAKE_C_FLAGS "-qno-opt-dynamic-align -std=gnu99 -fp-model precise")
+    set(CMAKE_C_FLAGS "-qno-opt-dynamic-align -std=gnu99 -fp-model precise -qopenmp")
     set(CMAKE_C_FLAGS_DEBUG "-O0 -g")
     set(CMAKE_C_FLAGS_RELEASE "-O2 -debug minimal")
-    set(CMAKE_Fortran_FLAGS "-free -qno-opt-dynamic-align -ftz -traceback -convert big_endian -assume byterecl -assume realloc_lhs -fp-model source")
+    set(CMAKE_Fortran_FLAGS "-free -qno-opt-dynamic-align -ftz -traceback -convert big_endian -assume byterecl -assume realloc_lhs -fp-model source -qopenmp")
     set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created")
     set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -debug minimal") 
 else()

namelist_generator/clm5nl/generators/gen_lnd_in.py

@@ -391,7 +391,7 @@ def setup_logic_lnd_frac():
 def setup_logic_co2_type():
     _nl.clm_inparm.co2_type = _opts["co2_type"]
     if _opts["co2_type"] == "constant":     
-        if _opts["co2_ppmv"] is None:
+        if "co2_ppmv" not in _opts:
             if _opts["sim_year"] == 2100:
                 ssp_co2_defaults = {"SSP5-8.5":1135.2, "SSP5-3.4":496.6, "SSP1-2.6":445.6}
                 _nl.clm_inparm.co2_ppmv = ssp_co2_defaults.get(_opts["ssp_rcp"], "Invalid ssp_rcp value")
@@ -420,7 +420,7 @@ def setup_logic_start_type():
         _nl.clm_inparm.override_nsrest = _opts["override_nsrest"]
 
     if my_start_type == "branch":
-        if _opts["nrevsn"] is None:
+        if "nrevsn" not in _opts:
             error("nrevsn is required for a branch type.")
         else:
             _nl.clm_inparm.nrevsn = _opts["nrevsn"]
@@ -433,12 +433,12 @@ def setup_logic_start_type():
         _nl.clm_inparm.use_init_interp = True
 
 def setup_logic_delta_time():
-    if _opts["l_ncpl"] is None:
+    if "l_ncpl" not in _opts:
         _nl.clm_inparm.dtime = 1800
     else:
         if _opts["l_ncpl"] <= 0: 
             error("bad value for -l_ncpl option")
-        if _opts["dtime"] is None:
+        if "dtime" not in _opts:
             _nl.clm_inparm.dtime = int((3600 * 24) / _opts["l_ncpl"])
         else:
             error("can NOT set both -l_ncpl option (via LND_NCPL env variable) AND dtime namelist variable.")
@@ -601,20 +601,20 @@ def setup_logic_surface_dataset():
 
 def setup_logic_initial_conditions():
     if _opts["clm_start_type"] == "cold":
-        if not _user_nl["finidat"] is None:
+        if "finidat" in _user_nl:
             print("""
             WARNING: setting finidat (either explicitly in your user_nl_clm or by doing a hybrid or branch RUN_TYPE) is 
             incomptable with using a cold start (by setting CLM_FORCE_COLDSTART=on)
             Overridding input finidat file with one specifying that this is a cold start from arbitrary initial conditions.""")
         _opts["finidat"] = "' '"
-    elif not _user_nl["finidat"] is None and _user_nl["finidat"] == "' '":
+    elif "finidat" in _user_nl and _user_nl["finidat"] == "' '":
         error("""You are setting finidat to blank which signals arbitrary initial conditions.
         But, CLM_FORCE_COLDSTART is off which is a contradiction. For arbitrary initial conditions just use the CLM_FORCE_COLDSTART option
         To do a cold-start set ./xmlchange CLM_FORCE_COLDSTART=on, and remove the setting of finidat in the user_nl_clm file""")
 
-    if _user_nl["finidat"] is None:
+    if "finidat" not in _user_nl:
         #TODO
-        pass
+        assert False, "No implementation yet for  \"if \"finidat\" not in _user_nl\" (TODO)"
     else:
         _nl.clm_inparm.finidat = _user_nl["finidat"]
 

src/clm5/cpl/lnd_comp_mct.F90

@@ -12,6 +12,10 @@ module lnd_comp_mct
   use mct_mod          , only : mct_avect, mct_gsmap, mct_gGrid
   use decompmod        , only : bounds_type, ldecomp
   use lnd_import_export, only : lnd_import, lnd_export
+#if defined(USE_OASIS)
+  use oas_defineMod     , only : oas_definitions_init
+  use oas_sendReceiveMod, only : oas_receive, oas_send
+#endif
   !
   ! !public member functions:
   implicit none
@@ -64,9 +68,6 @@ subroutine lnd_init_mct( EClock, cdata_l, x2l_l, l2x_l, NLFilename )
     use clm_varctl       , only : nsrStartup, nsrContinue, nsrBranch
     use clm_cpl_indices  , only : clm_cpl_indices_set
     use mct_mod          , only : mct_aVect_init, mct_aVect_zero, mct_gsMap_lsize
-#if defined(USE_OASIS)
-    use oas_defineMod    , only : oas_definitions_init
-#endif
     use ESMF
     !
     ! !ARGUMENTS:
@@ -331,6 +332,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
     integer      :: tod                  ! CLM current time of day (sec)
     integer      :: dtime                ! time step increment (sec)
     integer      :: nstep                ! time step index
+    integer      :: time_elapsed         ! elapsed time (dtime * nstep)
     logical      :: rstwr_sync           ! .true. ==> write restart file before returning
     logical      :: rstwr                ! .true. ==> write restart file before returning
     logical      :: nlend_sync           ! Flag signaling last time-step
@@ -411,7 +413,6 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
          atm2lnd_inst = atm2lnd_inst, &
          glc2lnd_inst = glc2lnd_inst)
     call t_stopf ('lc_lnd_import')
-
     ! Use infodata to set orbital values if updated mid-run
 
     call seq_infodata_GetData( infodata, orb_eccen=eccen, orb_mvelpp=mvelpp, &
@@ -434,9 +435,10 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
        ! Determine doalb based on nextsw_cday sent from atm model
 
        nstep = get_nstep()
+       time_elapsed = dtime * nstep
        caldayp1 = get_curr_calday(offset=dtime)
        if (nstep == 0) then
-	  doalb = .false. 	
+          doalb = .false.
        else if (nstep == 1) then 
           doalb = (abs(nextsw_cday- caldayp1) < 1.e-10_r8) 
        else
@@ -451,8 +453,10 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
        nlend = .false.
        if (nlend_sync .and. dosend) nlend = .true.
 
+#if defined(USE_OASIS)
+       call oas_receive(bounds, time_elapsed, atm2lnd_inst)
+#endif
        ! Run clm 
-
        call t_barrierf('sync_clm_run1', mpicom)
        call t_startf ('clm_run')
        call t_startf ('shr_orb_decl')
@@ -463,8 +467,10 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
        call clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, rof_prognostic)
        call t_stopf ('clm_run')
 
+#if defined(USE_OASIS)
+       call oas_send(bounds, time_elapsed, lnd2atm_inst)
+#endif
        ! Create l2x_l export state - add river runoff input to l2x_l if appropriate
-
        call t_startf ('lc_lnd_export')
        call lnd_export(bounds, lnd2atm_inst, lnd2glc_inst, l2x_l%rattr)
        call t_stopf ('lc_lnd_export')
@@ -474,7 +480,6 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l)
        call t_startf ('lc_clm2_adv_timestep')
        call advance_timestep()
        call t_stopf ('lc_clm2_adv_timestep')
-
     end do
 
     ! Check that internal clock is in sync with master clock

src/clm5/cpl/lnd_import_export.F90

@@ -8,9 +8,6 @@ module lnd_import_export
   use atm2lndType  , only: atm2lnd_type
   use glc2lndMod   , only: glc2lnd_type 
   use clm_cpl_indices
-#if defined(USE_OASIS)
-  use oas_sendReceiveMod
-#endif
   !
   implicit none
   !===============================================================================
@@ -279,9 +276,7 @@ subroutine lnd_import( bounds, x2l, glc_present, atm2lnd_inst, glc2lnd_inst)
          index_x2l_Flgg_hflx = index_x2l_Flgg_hflx, &
          index_x2l_Sg_icemask = index_x2l_Sg_icemask, &
          index_x2l_Sg_icemask_coupled_fluxes = index_x2l_Sg_icemask_coupled_fluxes)
-#if defined(USE_OASIS)
-    call oas_receive(bounds, atm2lnd_inst)
-#endif
+
   end subroutine lnd_import
 
   !===============================================================================
@@ -428,9 +423,7 @@ subroutine lnd_export( bounds, lnd2atm_inst, lnd2glc_inst, l2x)
        end if
 
     end do
-#if defined(USE_OASIS)
-    call oas_send(bounds, lnd2atm_inst)
-#endif
+
   end subroutine lnd_export
 
 end module lnd_import_export

src/clm5/oasis3/oas_defineMod.F90

@@ -10,17 +10,23 @@ module oas_defineMod
 
   subroutine oas_definitions_init(bounds)
     use spmdMod      , only : masterproc
-    use clm_varpar   , only : nlevsoi
-    use decompMod    , only : bounds_type
+    use clm_varpar   , only : nlevsoi, nlevgrnd
+    use decompMod    , only : bounds_type, ldecomp
+    use domainMod    , only : ldomain
     use oas_vardefMod
 
-    type(bounds_type) , intent(in)  :: bounds
-    integer          :: partition(3)
-    integer          :: grid_id       ! id returned by oasis_def_partition
-    integer          :: var_nodims(2)
-    integer          :: var_shape(1)  ! not used by oasis_def_var
-    character(len=2) :: soil_layer         
-    integer          :: i, n_grid_points
+    type(bounds_type) , intent(in)  :: bounds ! start and end gridcell indices for this MPI task
+
+    ! oasis_def_partition
+    integer, allocatable :: partition(:)      ! partition descriptor; input to oasis_def_partition
+    integer              :: gcell_start       ! starting gridcell index
+    integer              :: gcell_previous    ! gridcell index from previous loop iteration
+    integer              :: k, g              ! array/loop indices
+    integer              :: grid_id           ! id returned after call to oasis_def_partition
+
+    ! oasis_def_var
+    integer              :: var_nodims(2)     ! var dimension parameters
+    integer              :: var_shape(1)      ! unused dummy input to oasis_def_var
 
     if (masterproc) then
       call define_grid()
@@ -29,33 +35,51 @@ subroutine oas_definitions_init(bounds)
     ! -----------------------------------------------------------------
     ! ... Define partition
     ! -----------------------------------------------------------------
-    n_grid_points = (bounds%endg - bounds%begg) + 1
-    partition(1) = 1                ! Apple style partition
-    partition(2) = bounds%begg - 1  ! Global offset
-    partition(3) = n_grid_points    ! # of grid cells allocated to this MPI task
+    allocate(partition(200))
+    partition(:) = 0; k = 0
+
+    ! Use ORANGE partitioning scheme. This scheme defines an ensemble
+    ! of gridcell segments. See OASIS3-MCT User's guide for more info.
+    partition(1) = 3
+
+    ! Mark 1st segment
+    gcell_start = ldecomp%gdc2glo(bounds%begg)
+    partition(2) = 1
+    gcell_previous = gcell_start
+
+    ! Capture segments by detecting segment boundaries. A boundary is 
+    ! detected when the current and previous gridcells are not consecutive.
+    do g = bounds%begg+1, bounds%endg
+      if (ldecomp%gdc2glo(g) - gcell_previous /= 1) then
+        ! Previous segment complete; its partition params could now be defined
+        partition(3+k) = gcell_start - 1                  ! segment global offset (0-based)
+        partition(4+k) = gcell_previous - gcell_start + 1 ! segment length
+        k = k + 2
+
+        gcell_start  = ldecomp%gdc2glo(g) ! current gridcell marks the start of a new segment 
+        partition(2) = partition(2) + 1   ! increment number of segments
+      end if
+      gcell_previous = ldecomp%gdc2glo(g)
+    enddo
+
+    ! Define partition params for last segment
+    partition(3+k) = gcell_start - 1
+    partition(4+k) = gcell_previous - gcell_start + 1
+
     call oasis_def_partition(grid_id, partition, ierror)
+    deallocate(partition)
 
     ! -----------------------------------------------------------------
     ! ... Define coupling fields
     ! -----------------------------------------------------------------
-    var_nodims(1) = 1               ! var_nodims(1) is not used anymore in OASIS
-    var_nodims(2) = 1               ! number of fields in a bundle
+    var_nodims(1) = 1         ! unused
+    var_nodims(2) = nlevsoi   ! number of fields in a bundle
 
-    allocate(et_loss(nlevsoi))
-    allocate(watsat(nlevsoi))
-    allocate(psi(nlevsoi))
-
-    do i = 1, nlevsoi
-      write (soil_layer, '(I2.2)') i ! soil layer index (01-10)
-
-      et_loss(i)%name = 'CLMFLX'//soil_layer ! Evapotranspiration fluxes sent to Parflow
-      watsat(i)%name  = 'CLMSAT'//soil_layer ! Water saturation received from Parflow
-      psi(i)%name     = 'CLMPSI'//soil_layer ! Pressure head received from Parflow
-
-      call oasis_def_var(et_loss(i)%id, et_loss(i)%name, grid_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror) 
-      call oasis_def_var(watsat(i)%id, watsat(i)%name, grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
-      call oasis_def_var(psi(i)%id, psi(i)%name, grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
-    end do
+    call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror) 
+
+    var_nodims(2) = nlevgrnd  ! number of fields in a bundle
+    call oasis_def_var(oas_sat_id, "ECLM_SAT", grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
+    call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror)
 
     ! End definition phase
     call oasis_enddef(ierror)

src/clm5/oasis3/oas_sendReceiveMod.F90

@@ -1,8 +1,9 @@
 module oas_sendReceiveMod
   use shr_kind_mod     , only: r8 => shr_kind_r8
-  use clm_time_manager , only: get_curr_time, get_prev_time
+  use clm_time_manager , only: get_nstep, get_step_size
   use decompMod        , only: bounds_type
   use clm_varpar       , only: nlevsoi
+  use clm_varctl       , only: iulog
   use oas_vardefMod
   use mod_oasis
   implicit none
@@ -11,40 +12,40 @@ module oas_sendReceiveMod
 
   public  :: oas_send
   public  :: oas_receive
-  integer :: days_elapsed, seconds_elapsed
-  integer :: n_grid_points
-  integer :: i, ierror
 
 contains
 
-  subroutine oas_receive(bounds, atm2lnd_inst)
+  subroutine oas_receive(bounds, seconds_elapsed, atm2lnd_inst)
     use atm2lndType, only: atm2lnd_type
+
     type(bounds_type),  intent(in)    :: bounds
+    integer          ,  intent(in)    :: seconds_elapsed
     type(atm2lnd_type), intent(inout) :: atm2lnd_inst
-    real(kind=r8),      allocatable   :: buffer(:)
-
-    n_grid_points = (bounds%endg - bounds%begg) + 1
-    allocate(buffer(n_grid_points))
-    call get_curr_time(days_elapsed, seconds_elapsed)    
-    do i = 1, nlevsoi
-      call oasis_get(psi(i)%id, seconds_elapsed, buffer, ierror)
-      atm2lnd_inst%parflow_psi_grc(:,i) = buffer
-    end do
+    real(kind=r8),      allocatable   :: buffer(:,:)
+    integer                           :: num_grid_points
+    integer                           :: info
+
+    num_grid_points = (bounds%endg - bounds%begg) + 1
+    allocate(buffer(num_grid_points, nlevsoi))
+
+    call oasis_get(oas_psi_id, seconds_elapsed, atm2lnd_inst%parflow_psi_grc, info)
+    call oasis_get(oas_sat_id, seconds_elapsed, buffer, info)
+
   end subroutine oas_receive
 
-  subroutine oas_send(bounds, lnd2atm_inst)
+  subroutine oas_send(bounds, seconds_elapsed, lnd2atm_inst)
     use lnd2atmType, only : lnd2atm_type
+    use spmdMod,     only : mpicom
+    use shr_mpi_mod, only: shr_mpi_barrier
+
     type(bounds_type),  intent(in)    :: bounds
+    integer          ,  intent(in)    :: seconds_elapsed
     type(lnd2atm_type), intent(inout) :: lnd2atm_inst
-    real(kind=r8),      allocatable   :: buffer(:)
-
-    n_grid_points = (bounds%endg - bounds%begg) + 1
-    allocate(buffer(n_grid_points))
-    call get_prev_time(days_elapsed, seconds_elapsed)
-    do i = 1, nlevsoi
-      buffer = lnd2atm_inst%qflx_parflow_grc(:,i)
-      call oasis_put(et_loss(i)%id, seconds_elapsed, buffer, ierror)
-    end do
+
+    integer                           :: info
+
+    call oasis_put(oas_et_loss_id, seconds_elapsed, lnd2atm_inst%qflx_parflow_grc, info)
+
   end subroutine oas_send
 
 end module oas_sendReceiveMod

src/clm5/oasis3/oas_vardefMod.F90

@@ -2,17 +2,6 @@ module oas_vardefMod
   implicit none
   save
 
-  ! Type for coupling field information
-  type :: oas_var
-    character(len=12) :: name  ! Name of the coupling field
-    integer           :: id    ! Id of the field
-  end type oas_var
-
-  ! Sent fields
-  type(oas_var), allocatable :: et_loss(:) ! soil ET loss
-
-  ! Received fields from Parflow
-  type(oas_var), allocatable :: watsat(:)  ! water saturation
-  type(oas_var), allocatable :: psi(:)     ! pressure head
+  integer :: oas_psi_id, oas_et_loss_id, oas_sat_id
 
 end module oas_vardefMod