Refactor nems field exchange; set default masks for mapping in med_in…

…ternalstate (#279)

Refactors esmFldsExchange_nems.F90 to use separate advertise and initialize phases and to check that a component is present before advertising a field to or from that component. Implements default src and dst mask values in place of the code currently in med_map_mod.F90. Fixes #63 and #64.

mediator/med.F90

@@ -25,7 +25,6 @@ module MED
   use med_constants_mod        , only : spval_init         => med_constants_spval_init
   use med_constants_mod        , only : spval              => med_constants_spval
   use med_constants_mod        , only : czero              => med_constants_czero
-  use med_constants_mod        , only : ispval_mask        => med_constants_ispval_mask
   use med_utils_mod            , only : chkerr             => med_utils_ChkErr
   use med_methods_mod          , only : Field_GeomPrint    => med_methods_Field_GeomPrint
   use med_methods_mod          , only : State_GeomPrint    => med_methods_State_GeomPrint
@@ -41,7 +40,7 @@ module MED
   use med_utils_mod            , only : memcheck           => med_memcheck
   use med_time_mod             , only : med_time_alarmInit
   use med_internalstate_mod    , only : InternalState, med_internalstate_init, med_internalstate_coupling
-  use med_internalstate_mod    , only : logunit, mastertask
+  use med_internalstate_mod    , only : med_internalstate_defaultmasks, logunit, mastertask
   use med_internalstate_mod    , only : ncomps, compname
   use med_internalstate_mod    , only : compmed, compatm, compocn, compice, complnd, comprof, compwav, compglc
   use med_internalstate_mod    , only : coupling_mode
@@ -654,13 +653,14 @@ subroutine InitializeIPDv03p1(gcomp, importState, exportState, clock, rc)
     ! TransferOfferGeomObject Attribute.
 
     use ESMF  , only : ESMF_GridComp, ESMF_State, ESMF_Clock, ESMF_SUCCESS, ESMF_LogFoundAllocError
-    use ESMF  , only : ESMF_StateIsCreated 
+    use ESMF  , only : ESMF_StateIsCreated
     use ESMF  , only : ESMF_LogMsg_Info, ESMF_LogWrite
     use ESMF  , only : ESMF_END_ABORT, ESMF_Finalize, ESMF_MAXSTR
     use NUOPC , only : NUOPC_AddNamespace, NUOPC_Advertise, NUOPC_AddNestedState
     use NUOPC , only : NUOPC_CompAttributeGet, NUOPC_CompAttributeSet, NUOPC_CompAttributeAdd
     use esmFlds, only : med_fldlist_init1
     use med_phases_history_mod, only : med_phases_history_init
+    use med_internalstate_mod , only : atm_name
 
     ! input/output variables
     type(ESMF_GridComp)  :: gcomp
@@ -783,8 +783,7 @@ subroutine InitializeIPDv03p1(gcomp, importState, exportState, clock, rc)
     if (trim(coupling_mode) == 'cesm') then
        call esmFldsExchange_cesm(gcomp, phase='advertise', rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    else if (trim(coupling_mode) == 'nems_orig' .or. trim(coupling_mode) == 'nems_frac' &
-       .or. trim(coupling_mode) == 'nems_orig_data') then
+    else if (trim(coupling_mode(1:4)) == 'nems') then
        call esmFldsExchange_nems(gcomp, phase='advertise', rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
     else if (trim(coupling_mode(1:4)) == 'hafs') then
@@ -795,6 +794,10 @@ subroutine InitializeIPDv03p1(gcomp, importState, exportState, clock, rc)
         call ESMF_Finalize(endflag=ESMF_END_ABORT)
     end if
 
+    ! Set default masking for mapping
+    call med_internalstate_defaultmasks(gcomp, rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
     !------------------
     ! Determine component present indices
     !------------------
@@ -1746,6 +1749,8 @@ subroutine DataInitialize(gcomp, rc)
       if (trim(coupling_mode) == 'cesm') then
          call esmFldsExchange_cesm(gcomp, phase='initialize', rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
+      else if (trim(coupling_mode(1:4)) == 'nems') then
+       call esmFldsExchange_nems(gcomp, phase='initialize', rc=rc)
       else if (trim(coupling_mode) == 'hafs') then
          call esmFldsExchange_hafs(gcomp, phase='initialize', rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return

mediator/med_internalstate_mod.F90

@@ -15,6 +15,7 @@ module med_internalstate_mod
   ! public routines
   public :: med_internalstate_init
   public :: med_internalstate_coupling
+  public :: med_internalstate_defaultmasks
 
   integer, public :: logunit            ! logunit for mediator log output
   integer, public :: diagunit           ! diagunit for budget output (med master only)
@@ -48,6 +49,9 @@ module med_internalstate_mod
   ! Coupling mode
   character(len=CS), public :: coupling_mode ! valid values are [cesm,nems_orig,nems_frac,nems_orig_data,hafs]
 
+  ! Default src and destination masks for mapping
+  integer, public, allocatable :: defaultMasks(:,:)
+
   ! Mapping
   integer , public, parameter :: mapunset          = 0
   integer , public, parameter :: mapbilnr          = 1
@@ -113,7 +117,7 @@ module med_internalstate_mod
     logical, pointer :: med_coupling_active(:,:)   ! computes the active coupling
     integer          :: num_icesheets              ! obtained from attribute
     logical          :: ocn2glc_coupling = .false. ! obtained from attribute
-    logical          :: lnd2glc_coupling = .false. 
+    logical          :: lnd2glc_coupling = .false.
     logical          :: accum_lnd2glc = .false.
 
     ! Mediator vm
@@ -187,8 +191,8 @@ module med_internalstate_mod
 
   subroutine med_internalstate_init(gcomp, rc)
 
-    use ESMF         , only : ESMF_LogFoundAllocError, ESMF_AttributeGet
-    use NUOPC_Comp   , only : NUOPC_CompAttributeGet
+    use ESMF              , only : ESMF_LogFoundAllocError, ESMF_AttributeGet
+    use NUOPC_Comp        , only : NUOPC_CompAttributeGet
 
     ! input/output variables
     type(ESMF_GridComp)            :: gcomp
@@ -205,7 +209,7 @@ subroutine med_internalstate_init(gcomp, rc)
     character(len=CL)          :: cname
     character(len=ESMF_MAXSTR) :: mesh_glc
     character(len=CX)          :: msgString
-    character(len=3)           :: name 
+    character(len=3)           :: name
     integer                    :: num_icesheets
     character(len=*),parameter :: subname=' (internalstate init) '
     !-----------------------------------------------------------
@@ -329,7 +333,7 @@ subroutine med_internalstate_init(gcomp, rc)
        ! Write out present flags
        write(logunit,*)
        do n1 = 1,ncomps
-          name = trim(compname(n1))  ! this trims the ice sheets index from the glc name 
+          name = trim(compname(n1))  ! this trims the ice sheets index from the glc name
           write(msgString,'(A,L4)') trim(subname)//' comp_present(comp'//name//') = ',&
                is_local%wrap%comp_present(n1)
           write(logunit,'(a)') trim(msgString)
@@ -353,7 +357,7 @@ subroutine med_internalstate_init(gcomp, rc)
     ! Obtain dststatus_print setting if present
     call NUOPC_CompAttributeGet(gcomp, name='dststatus_print', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    if (isPresent .and. isSet) dststatus_print=(trim(cvalue)=="true")
+    if (isPresent .and. isSet) dststatus_print=(trim(cvalue) == "true")
     write(msgString,*) trim(subname)//': Mediator dststatus_print is ',dststatus_print
     call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO)
 
@@ -551,4 +555,44 @@ subroutine med_internalstate_coupling(gcomp, rc)
 
   end subroutine med_internalstate_coupling
 
+  subroutine med_internalstate_defaultmasks(gcomp, rc)
+
+    use med_constants_mod        , only : ispval_mask        => med_constants_ispval_mask
+
+    ! input/output variables
+    type(ESMF_GridComp)            :: gcomp
+    integer          , intent(out) :: rc
+
+    ! local variables
+    type(InternalState)        :: is_local
+
+    !----------------------------------------------------------
+    ! Default masking: for each component, the first element is
+    ! when it is the src and the second element is when it is
+    ! the destination
+    !----------------------------------------------------------
+
+    nullify(is_local%wrap)
+    call ESMF_GridCompGetInternalState(gcomp, is_local, rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    allocate(defaultMasks(ncomps,2))
+    defaultMasks(:,:) = ispval_mask
+    if (is_local%wrap%comp_present(compocn)) defaultMasks(compocn,:) = 0
+    if (is_local%wrap%comp_present(compice)) defaultMasks(compice,:) = 0
+    if (is_local%wrap%comp_present(compwav)) defaultMasks(compwav,:) = 0
+    if ( trim(coupling_mode(1:4)) == 'nems') then
+       if (is_local%wrap%comp_present(compatm)) defaultMasks(compatm,:) = 1
+    endif
+    if ( trim(coupling_mode) == 'hafs') then
+       if (is_local%wrap%comp_present(compatm)) defaultMasks(compatm,1) = 1
+    endif
+    if ( trim(coupling_mode) /= 'cesm') then
+       if (is_local%wrap%comp_present(compatm) .and. trim(atm_name(1:4)) == 'datm') then
+          defaultMasks(compatm,1) = 0
+       end if
+    end if
+
+  end subroutine med_internalstate_defaultmasks
+
 end module med_internalstate_mod

mediator/med_map_mod.F90

@@ -342,7 +342,7 @@ subroutine med_map_routehandles_initfrom_field(n1, n2, fldsrc, flddst, mapindex,
     use med_internalstate_mod , only : mapfillv_bilnr, mapbilnr_nstod, mapconsf_aofrac
     use med_internalstate_mod , only : ncomps, compatm, compice, compocn, compwav, compname
     use med_internalstate_mod , only : coupling_mode, dststatus_print
-    use med_internalstate_mod , only : atm_name
+    use med_internalstate_mod , only : defaultMasks
     use med_constants_mod     , only : ispval_mask => med_constants_ispval_mask
 
     ! input/output variables
@@ -389,63 +389,33 @@ subroutine med_map_routehandles_initfrom_field(n1, n2, fldsrc, flddst, mapindex,
     ! set local flag to false
     ldstprint = .false.
 
-    polemethod=ESMF_POLEMETHOD_ALLAVG
+    ! set src and dst masking using defaults
+    srcMaskValue = defaultMasks(n1,1)
+    dstMaskValue = defaultMasks(n2,2)
+
+    ! override defaults for specific cases
     if (trim(coupling_mode) == 'cesm') then
-       dstMaskValue = ispval_mask
-       srcMaskValue = ispval_mask
-       if (n1 == compocn .or. n1 == compice) srcMaskValue = 0
-       if (n2 == compocn .or. n2 == compice) dstMaskValue = 0
        if (n1 == compwav .and. n2 == compocn) then
          srcMaskValue = 0
          dstMaskValue = ispval_mask
       endif
-      if (n1 == compwav .or. n2 == compwav) then
-        polemethod = ESMF_POLEMETHOD_NONE ! todo: remove this when ESMF tripolar mapping fix is in place.
-      endif
-    else if (coupling_mode(1:4) == 'nems') then
-       if ( (n1 == compocn .or. n1 == compice .or. n1 == compwav) .and. &
-            (n2 == compocn .or. n2 == compice .or. n2 == compwav) ) then
-          srcMaskValue = 0
-          dstMaskValue = 0
-       else if (n1 == compatm .and. (n2 == compocn .or. n2 == compice .or. n2 == compwav)) then
-          srcMaskValue = 1
-          dstMaskValue = 0
-          if (atm_name(1:4).eq.'datm') then
-             srcMaskValue = 0
-          endif
-       else if (n2 == compatm .and. (n1 == compocn .or. n1 == compice .or. n1 == compwav)) then
-          srcMaskValue = 0
-          dstMaskValue = 1
-       else
-          ! TODO: what should the condition be here?
-          dstMaskValue = ispval_mask
+    end if
+    if (trim(coupling_mode) == 'hafs') then
+       if (n1 == compatm .and. n2 == compwav) then
           srcMaskValue = ispval_mask
        end if
-    else if (trim(coupling_mode) == 'hafs') then
-       dstMaskValue = ispval_mask
-       srcMaskValue = ispval_mask
-       if (n1 == compocn .or. n1 == compice) srcMaskValue = 0
-       if (n2 == compocn .or. n2 == compice) dstMaskValue = 0
-       if (n1 == compatm .and. n2 == compocn) then
-          if (trim(atm_name).ne.'datm') then
-             srcMaskValue = 1
-          endif
-          dstMaskValue = 0
-       elseif (n1 == compocn .and. n2 == compatm) then
-          srcMaskValue = 0
-          dstMaskValue = ispval_mask
-       elseif (n1 == compatm .and. n2 == compwav) then
-          dstMaskValue = 0
-       elseif (n1 == compwav .and. n2 == compatm) then
-          srcMaskValue = 0
-          dstMaskValue = ispval_mask
-       endif
     end if
-
     write(string,'(a,i10,a,i10)') trim(compname(n1))//' to '//trim(compname(n2))//' srcMask = ', &
                srcMaskValue,' dstMask = ',dstMaskValue
     call ESMF_LogWrite(trim(string), ESMF_LOGMSG_INFO)
 
+    polemethod=ESMF_POLEMETHOD_ALLAVG
+    if (trim(coupling_mode) == 'cesm') then
+      if (n1 == compwav .or. n2 == compwav) then
+        polemethod = ESMF_POLEMETHOD_NONE ! todo: remove this when ESMF tripolar mapping fix is in place.
+      endif
+    end if
+
     ! Create route handle
     if (mapindex == mapfcopy) then
        if (mastertask) then

mediator/med_phases_post_lnd_mod.F90

@@ -80,7 +80,7 @@ subroutine med_phases_post_lnd(gcomp, rc)
        if (is_local%wrap%lnd2glc_coupling) then
           call med_phases_prep_glc_accum_lnd(gcomp, rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
-          ! Note that in this case med_phases_prep_glc_avg is called 
+          ! Note that in this case med_phases_prep_glc_avg is called
           ! from med_phases_prep_glc in the run sequence
        else if (is_local%wrap%accum_lnd2glc) then
           call med_phases_prep_glc_accum_lnd(gcomp, rc)

mediator/med_phases_prep_atm_mod.F90

@@ -242,17 +242,17 @@ subroutine med_phases_prep_atm(gcomp, rc)
   end subroutine med_phases_prep_atm
 
   !-----------------------------------------------------------------------------
-  subroutine med_phases_prep_atm_enthalpy_correction (gcomp, hcorr, rc) 
+  subroutine med_phases_prep_atm_enthalpy_correction (gcomp, hcorr, rc)
 
-    ! Enthalpy correction term calculation called by med_phases_prep_ocn_accum in 
+    ! Enthalpy correction term calculation called by med_phases_prep_ocn_accum in
     ! med_phases_prep_ocn_mod
     ! Note that this is only called if the following fields are in FBExp(compocn)
     ! 'Faxa_rain','Foxx_hrain','Faxa_snow' ,'Foxx_hsnow',
-    ! 'Foxx_evap','Foxx_hevap','Foxx_hcond','Foxx_rofl', 
+    ! 'Foxx_evap','Foxx_hevap','Foxx_hcond','Foxx_rofl',
     ! 'Foxx_hrofl','Foxx_rofi','Foxx_hrofi'
 
     use ESMF            , only : ESMF_VMAllreduce, ESMF_GridCompGet, ESMF_REDUCE_SUM
-    use ESMF            , only : ESMF_VM 
+    use ESMF            , only : ESMF_VM
 
     ! input/output variables
     type(ESMF_GridComp) , intent(in)  :: gcomp

mediator/med_time_mod.F90

@@ -73,7 +73,7 @@ subroutine med_time_alarmInit( clock, alarm, option, &
     integer          , optional , intent(in)    :: opt_tod       ! alarm tod (sec)
     type(ESMF_Time)  , optional , intent(in)    :: reftime       ! reference time
     character(len=*) , optional , intent(in)    :: alarmname     ! alarm name
-    logical          , optional , intent(in)    :: advance_clock ! advance clock to trigger alarm 
+    logical          , optional , intent(in)    :: advance_clock ! advance clock to trigger alarm
     integer                     , intent(out)   :: rc            ! Return code
 
     ! local variables
@@ -264,7 +264,7 @@ subroutine med_time_alarmInit( clock, alarm, option, &
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     ! Advance model clock to trigger alarm then reset model clock back to currtime
-    if (present(advance_clock)) then 
+    if (present(advance_clock)) then
        if (advance_clock) then
           call ESMF_AlarmSet(alarm, clock=clock, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return