wmesmfmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/
!/ ------------------------------------------------------------------- /
!/ Macros for ESMF logging
!/
#define FILENAME "wmesmfmd.ftn"
#define CONTEXT  line=__LINE__,file=FILENAME,method=METHOD
#define PASSTHRU msg=ESMF_LOGERR_PASSTHRU,CONTEXT
!/
!/ ------------------------------------------------------------------- /
!/ Define real kind for data passed through ESMF interface
!/
#if defined(ESMF_R8)
#define ESMF_KIND_RX ESMF_KIND_R8
#define ESMF_TYPEKIND_RX ESMF_TYPEKIND_R8
#else
#define ESMF_KIND_RX ESMF_KIND_R4
#define ESMF_TYPEKIND_RX ESMF_TYPEKIND_R4
#endif
!/
!/ ------------------------------------------------------------------- /
!/ Macro for enabling using W3OUTG to calculate export fields
!/
#define USE_W3OUTG_FOR_EXPORT___disabled
!/
!/ ------------------------------------------------------------------- /
!/ Macros for enabling test output
!/
#define TEST_WMESMFMD___disabled
#define TEST_WMESMFMD_GETIMPORT___disabled
#define TEST_WMESMFMD_SETEXPORT___disabled
#define TEST_WMESMFMD_CREATEIMPGRID___disabled
#define TEST_WMESMFMD_CREATEEXPGRID___disabled
#define TEST_WMESMFMD_SETUPIMPBMSK___disabled
#define TEST_WMESMFMD_SETUPIMPMMSK___disabled
#define TEST_WMESMFMD_CHARNK___disabled
#define TEST_WMESMFMD_ROUGHL___disabled
#define TEST_WMESMFMD_BOTCUR___disabled
#define TEST_WMESMFMD_RADSTR2D___disabled
#define TEST_WMESMFMD_STOKES3D___disabled
#define TEST_WMESMFMD_PSTOKES___disabled
#define TEST_WMESMFMD_READFROMFILE___disabled
 
!/ ------------------------------------------------------------------- /
module wmesmfmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |        T. J. Campbell, NRL        |
  !/                  |         J. Meixner, NCEP          |
  !/                  |     A. J. van der Westhuysen      |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         09-Aug-2017 |
  !/                  +-----------------------------------+
  !/
  !/    20-Jan-2017 : Origination.                        ( version 6.02 )
  !/    09-Aug-2017 : Add ocean forcing export fields     ( version 6.03 )
  !/    28-Feb-2018 : Modifications for unstruc meshes    ( version 6.06 )
  !/
  !/    Copyright 2009-2014 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     National Unified Prediction Capability (NUOPC) based
  !     Earth System Modeling Framework (ESMF) interface module for
  !     multi-grid wave model.
  !
  !  2. Variables and types :
  !
  !     All module variables and types are scoped private by default.
  !     The private module variables and types are not listed in this section.
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !     All module subroutines and functions are scoped private by default.
  !
  !      Name             Type   Scope    Description
  !     -----------------------------------------------------------------
  !      SetServices      Subr.  Public   Wave model ESMF Set Services
  !     -----------------------------------------------------------------
  !      InitializeP0     Subr.  Private  NUOPC/ESMF Initialize phase 0
  !      InitializeP1     Subr.  Private  NUOPC/ESMF Initialize phase 1
  !      InitializeP3     Subr.  Private  NUOPC/ESMF Initialize phase 3
  !      Finalize         Subr.  Private  NUOPC/ESMF Finalize
  !      DataInitialize   Subr.  Private  NUOPC/ESMF Data Initialize
  !      ModelAdvance     Subr.  Private  NUOPC/ESMF Model Advance
  !      GetImport        Subr.  Private  Get fields from import state
  !      SetExport        Subr.  Private  Set fields from export state
  !      CreateImpGrid    Subr.  Private  Create ESMF grid for import
  !      CreateExpGrid    Subr.  Private  Create ESMF grid for export
  !      CreateImpMesh    Subr.  Private  Create ESMF mesh for import
  !      CreateExpMesh    Subr.  Private  Create ESMF mesh for export
  !      SetupImpBmsk     Subr.  Private  Setup background blending mask
  !      BlendImpField    Subr.  Private  Blend import field with background
  !      SetupImpMmsk     Subr.  Private  Setup merging mask
  !      FieldFill        Subr.  Private  Fill ESMF field
  !      FieldGather      Subr.  Private  Gather ESMF field
  !      FieldIndex       Func.  Private  Return field index
  !      PrintTimers      Subr.  Private  Print wallclock timers
  !      CalcDecomp       Subr.  Private  Calculate a 2D processor layout
  !      GetEnvValue      Subr.  Private  Get value of env. variable
  !      GetZlevels       Subr.  Private  Get z-levels from file for SDC
  !      CalcCharnk       Subr.  Private  Calculate Charnock for export
  !      CalcRoughl       Subr.  Private  Calculate roughness length for export
  !      CalcBotcur       Subr.  Private  Calculate wave-bottom currents for export
  !      CalcRadstr2D     Subr.  Private  Calculate 2D radiation stresses for export
  !      CalcStokes3D     Subr.  Private  Calculate 3D Stokes drift current for export
  !      CalcPStokes      Subr.  Private  Calculate partitioned Stokes drift for export
  !      ReadFromFile     Subr.  Private  Read input file
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !     See subroutine documentation.
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !     See subroutine documentation.
  !
  !     !/MPI   Switch for enabling Message Passing Interface API
  !     !/SHRD  Switch for shared memory architecture
  !     !/DIST  Switch for distributed memory architecture
  !     !/ST3   WAM 4+ input and dissipation.
  !     !/ST4   Ardhuin et al. (2009, 2010)
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Use associated modules
  !/
  ! --- ESMF Module
  use esmf
 
  ! --- NUOPC modules
  use nuopc
  use nuopc_model, parent_setservices => setservices
 
  ! --- WW3 modules
  use constants
  use wminitmd, only: wminit, wminitnml
  use wmwavemd, only: wmwave
  use wmfinlmd, only: wmfinl
  use wmmdatmd
  use w3gdatmd
  use w3idatmd
  use w3odatmd
  use w3wdatmd
  use w3adatmd
  use w3timemd
  use wmupdtmd, only: wmupd2
  use w3updtmd, only: w3uini
#ifdef W3_ST3
  use w3src3md, only: w3spr3
#endif
#ifdef W3_ST4
  use w3src4md, only: w3spr4
#endif
  use w3iogomd, only: w3outg
#ifdef W3_SCRIP
  use wmscrpmd, only: get_scrip_info_structured
#endif
  !/
  !/ Specify default data typing
  !/
  implicit none
  !/
  !/ Include MPI definitions
  !/
#ifdef W3_MPI
  include "mpif.h"
#endif
  !/
  !/ Specify default accessibility
  !/
  private
  save
  !/
  !/ Public module methods
  !/
  public setservices, setvm
  !/
  !/ Private module parameters
  !/
  ! --- Default Mask Convention for import/export fields
  INTEGER, PARAMETER :: DEFAULT_MASK_WATER =  0
  INTEGER, PARAMETER :: DEFAULT_MASK_LAND  =  1
 
  ! --- Miscellaneous
  integer, parameter :: stdo = 6
  type(ESMF_VM) :: vm
  integer :: lpet
  integer :: npet
  integer :: verbosity
  logical :: realizeAllExport = .false.           
  integer :: maskValueWater = default_mask_water  
  integer :: maskValueLand  = default_mask_land   
  integer              :: nz
  real(4), allocatable :: zl(:)
  character(256)       :: zlfile = 'none'
  character(ESMF_MAXSTR) :: msg
  real(ESMF_KIND_RX) :: zeroValue
  real(ESMF_KIND_RX) :: missingValue
  real(ESMF_KIND_RX) :: fillValue
  !
  ! --- Timing
  integer, parameter :: numwt=10
  character(32) :: wtnam(numwt)
  integer       :: wtcnt(numwt)
  real(8)       :: wtime(numwt)
  !
  ! --- Import fields
  type(ESMF_ArraySpec)          :: impArraySpec2D
  type(ESMF_StaggerLoc)         :: impStaggerLoc
  type(ESMF_Index_Flag)         :: impIndexFlag
  type(ESMF_Grid)               :: impGrid
  integer                       :: impGridID
  logical                       :: impGridIsLocal
  integer, parameter            :: impHaloWidth = 3
  integer                       :: impHaloLWidth(2)
  integer                       :: impHaloUWidth(2)
  type(ESMF_RouteHandle)        :: impHaloRH
  type(ESMF_Field)              :: impMask
  logical                       :: noActiveImpFields
  integer                       :: numImpFields
  character(64), allocatable    :: impFieldName(:)
  character(128), allocatable   :: impFieldStdName(:)
  logical, allocatable          :: impFieldInitRqrd(:)
  logical, allocatable          :: impFieldActive(:)
  type(ESMF_Field), allocatable :: impField(:)
  !
  ! --- Background import fields
  character(10), allocatable    :: mbgFieldName(:)
  character(128), allocatable   :: mbgFieldStdName(:)
  logical, allocatable          :: mbgFieldActive(:)
  type(ESMF_Field), allocatable :: mbgField(:)
  type(ESMF_Field), allocatable :: bmskField(:)
  !
  ! --- Unstructured import meshes
  type(ESMF_Mesh)               :: impMesh
  !      integer                       :: impMeshID         !< impMeshID
  !      logical                       :: impMeshIsLocal         !< impMeshIsLocal
  !
  ! --- Export fields
  type(ESMF_ArraySpec)          :: expArraySpec2D
  type(ESMF_ArraySpec)          :: expArraySpec3D
  type(ESMF_StaggerLoc)         :: expStaggerLoc
  type(ESMF_Index_Flag)         :: expIndexFlag
  type(ESMF_Grid)               :: expGrid
  integer                       :: expGridID = 1
  logical                       :: expGridIsLocal
  integer, parameter            :: expHaloWidth = 3
  integer                       :: expHaloLWidth(2)
  integer                       :: expHaloUWidth(2)
  type(ESMF_RouteHandle)        :: expHaloRH
  type(ESMF_Field)              :: expMask
  logical                       :: noActiveExpFields
  integer                       :: numExpFields
  character(64), allocatable    :: expFieldName(:)
  character(128), allocatable   :: expFieldStdName(:)
  integer, allocatable          :: expFieldDim(:)
  logical, allocatable          :: expFieldActive(:)
  type(ESMF_Field), allocatable :: expField(:)
  !
  ! --- Unstructured export meshes
  type(ESMF_Mesh)               :: expMesh
  integer                       :: expMeshID
  logical                       :: expMeshIsLocal
  !
  ! --- Native field stuff
  type(ESMF_ArraySpec)  :: natArraySpec1D
  type(ESMF_ArraySpec)  :: natArraySpec2D
  type(ESMF_ArraySpec)  :: natArraySpec3D
  type(ESMF_StaggerLoc) :: natStaggerLoc
  type(ESMF_Index_Flag) :: natIndexFlag
  type(ESMF_Grid)       :: natGrid
  integer               :: natGridID
  logical               :: natGridIsLocal
  type(ESMF_RouteHandle):: n2eRH
  !
  ! --- Mediator
  logical        :: med_present = .false.         
  character(256) :: flds_scalar_name = ''
  integer        :: flds_scalar_num = 0
  ! flds_scalar_index_nx and flds_scalar_index_nx are domain
  ! metadata that allows CMEPS to convert a mesh back to 2d
  ! space for mediator restart and history outputs
  integer        :: flds_scalar_index_nx = 0
  integer        :: flds_scalar_index_ny = 0
  ! --- Memory Profiling
  logical        :: profile_memory = .false.         
  !
  ! --- Coupling stuff for non completely overlapped domains
  logical                       :: merge_import = .false.         
  logical, allocatable          :: mmskCreated(:)
  type(ESMF_Field), allocatable :: mmskField(:)
  type(ESMF_Field), allocatable :: mdtField(:)
  !/
  !/ ------------------------------------------------------------------- /
 
contains
 
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "SetServices"
 
  subroutine setservices ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Wave model ESMF set services.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp   Type   I/O Gridded component
    !       rc      Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name            Type   Module    Description
    !     ----------------------------------------------------------------
    !      InitializeP0    Subr.  WMESMFMD  Wave model NUOPC/ESMF Initialize phase 0
    !      InitializeP1    Subr.  WMESMFMD  Wave model NUOPC/ESMF Initialize phase 1
    !      InitializeP3    Subr.  WMESMFMD  Wave model NUOPC/ESMF Initialize phase 3
    !      Finalize        Subr.  WMESMFMD  Wave model NUOPC/ESMF Finalize
    !      DataInitialize  Subr.  WMESMFMD  Wave model NUOPC/ESMF Data Initialize
    !      ModelAdvance    Subr.  WMESMFMD  Wave model NUOPC/ESMF Model Advance
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_gridcomp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    !NONE
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = esmf_success
 
    ! --- Initialize wallclock timers
 
    wtnam( 1) = 'InitializeP0'
    wtnam( 2) = 'InitializeP1'
    wtnam( 3) = 'InitializeP3'
    wtnam( 4) = 'DataInitialize'
    wtnam( 5) = 'ModelAdvance'
    wtnam( 6) = 'Finalize'
    wtnam( 7) = 'GetImport'
    wtnam( 8) = 'SetExport'
    wtnam( 9) = 'FieldGather'
    wtnam(10) = 'FieldFill'
    wtcnt( :) = 0
    wtime( :) = 0d0
    !
    ! -------------------------------------------------------------------- /
    ! 1.  NUOPC model component will register the generic methods
    !
    call nuopc_compderive(gcomp, parent_setservices, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Set model entry points
    !
    ! --- Initialize - phase 0 (requires use of ESMF method)
 
    call esmf_gridcompsetentrypoint(gcomp, esmf_method_initialize, &
         userroutine=initializep0, phase=0, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    ! --- Set entry points for initialize methods
 
    ! >= IPDv03 supports satisfying inter-model data dependencies and
    ! the transfer of ESMF Grid & Mesh objects between Model and/or
    ! Mediator components during initialization
    ! IPDv03p1: advertise import & export fields
    call nuopc_compsetentrypoint(gcomp, esmf_method_initialize, &
         phaselabellist=(/"IPDv03p1"/), userroutine=initializep1, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    ! IPDv03p2: unspecified by NUOPC -- not required
    ! IPDv03p3: realize import & export fields
    call nuopc_compsetentrypoint(gcomp, esmf_method_initialize, &
         phaselabellist=(/"IPDv03p3"/), userroutine=initializep3, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    ! IPDv03p4: relevant for TransferActionGeomObject=="accept"
    ! IPDv03p5: relevant for TransferActionGeomObject=="accept"
    ! IPDv03p6: check compatibility of fields connected status
    ! IPDv03p7: handle field data initialization
 
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Register specializing methods
    !
    ! --- Model initialize export data method
 
    call nuopc_compspecialize(gcomp, speclabel=label_datainitialize, &
         specroutine=datainitialize, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    ! --- Model checkImport method (overriding default)
 
    call esmf_methodremove(gcomp, label_checkimport, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    call nuopc_compspecialize(gcomp, speclabel=label_checkimport, &
         specroutine=nuopc_noop, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    ! --- Model advance method
 
    call nuopc_compspecialize(gcomp, speclabel=label_advance, &
         specroutine=modeladvance, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    ! --- Model finalize method
 
    call nuopc_compspecialize(gcomp, speclabel=label_finalize, &
         specroutine=finalize, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    rc = esmf_success
    !/
    !/ End of SetServices ------------------------------------------------ /
    !/
  end subroutine setservices
  !/ ------------------------------------------------------------------- /
 
#undef METHOD
#define METHOD "InitializeP0"
 
  subroutine initializep0 ( gcomp, impState, expState, extClock, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Initialize wave model (phase 0)
    !     * Define the NUOPC Initialize Phase Mapping
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       impState  Type   I/O Import state
    !       expState  Type   I/O Export state
    !       extClock  Type   I   External clock
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_gridcomp) :: gcomp
    type(esmf_state)    :: impstate
    type(esmf_state)    :: expstate
    type(esmf_clock)    :: extclock
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    character(ESMF_MAXSTR) :: valuestring
    integer, parameter :: iwt=1
    real(8) :: wstime, wftime
    logical :: ispresent, isset
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = esmf_success
    call esmf_vmwtime(wstime)
    call esmf_gridcompget(gcomp, name=cname, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    !
    ! -------------------------------------------------------------------- /
    ! Determine verbosity
    !
    call nuopc_compattributeget(gcomp, name='Verbosity', &
         value=valuestring, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    verbosity = esmf_utilstring2int( valuestring, &
         specialstringlist=(/'high','max '/), &
         specialvaluelist=(/    255,   255/), rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': entered InitializeP0', esmf_logmsg_info)
    !
    ! -------------------------------------------------------------------- /
    ! Define initialization phases
    ! * switch to IPDv03 by filtering all other phaseMap entries
    !
    call nuopc_compfilterphasemap(gcomp, esmf_method_initialize, &
         acceptstringlist=(/"IPDv03p"/), rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    !
    ! -------------------------------------------------------------------- /
    ! Check if coupled with CMEPS mediator or not
    !
    call nuopc_compattributeget(gcomp, name="mediator_present", &
         value=valuestring, ispresent=ispresent, isset=isset, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    if (ispresent .and. isset) then
      read(valuestring,*) med_present
      call esmf_logwrite(trim(cname)//': mediator_present = '// &
           trim(valuestring), esmf_logmsg_info, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
    end if
    !
    ! -------------------------------------------------------------------- /
    ! Set memory profiling
    !
    call nuopc_compattributeget(gcomp, name="ProfileMemory", &
         value=valuestring, ispresent=ispresent, isset=isset, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
 
    if (ispresent .and. isset) then
      read(valuestring,*) profile_memory
      call esmf_logwrite(trim(cname)//': profile_memory = '// &
           trim(valuestring), esmf_logmsg_info, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
    end if
 
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    rc = esmf_success
    call esmf_vmwtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': leaving InitializeP0', esmf_logmsg_info)
    !/
    !/ End of InitializeP0 ----------------------------------------------- /
    !/
  end subroutine initializep0
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "InitializeP1"
 
  subroutine initializep1 ( gcomp, impState, expState, extClock, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    09-Aug-2017 : Add ocean forcing export fields     ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Initialize wave model (phase 1)
    !     * Advertise fields in import and export states.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       impState  Type   I/O Import state
    !       expState  Type   I/O Export state
    !       extClock  Type   I   External clock
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name            Type  Module   Description
    !     ----------------------------------------------------------------
    !      WMINIT          Subr. WMINITMD Wave model initialization
    !      WMINITNML       Subr. WMINITMD Wave model initialization
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     ----------------------------------------------------------------
    !      1.  Initialization necessary for driver
    !        a General I/O: (implicit in WMMDATMD)
    !        b MPI environment
    !        c Identifying output to "screen" unit
    !      2.  Initialization of all wave models / grids
    !      3.  Advertise import fields
    !      4.  Advertise export fields
    !     ----------------------------------------------------------------
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_gridcomp) :: gcomp
    type(esmf_state)    :: impstate
    type(esmf_state)    :: expstate
    type(esmf_clock)    :: extclock
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=2
    real(8) :: wstime, wftime
    integer :: idsi, idso, idss, idst, idse
    integer :: mpicomm = -99
    logical :: configispresent
    type(esmf_config) :: config
    character(ESMF_MAXSTR) :: wrkdir = '.'
    character(ESMF_MAXSTR) :: preamb = '.'
    character(ESMF_MAXSTR) :: ifname = 'ww3_multi.inp'
    logical :: lsep_ss = .true.
    logical :: lsep_st = .true.
    logical :: lsep_se = .true.
    character(ESMF_MAXSTR) :: attstr
    integer(ESMF_KIND_I4) :: yy,mm,dd,h,m,s
    type(esmf_time) :: ttmp, cttmp
    type(esmf_timeinterval) :: tstep, etstep
    integer :: i, j, n, istep, imod, jmod
    integer, allocatable :: cplmap(:,:)
    logical :: includeobg, includeabg, includeibg
    character(256) :: cvalue, logmsg
    logical :: ispresent, isset
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = esmf_success
    call esmf_vmwtime(wstime)
    call esmf_gridcompget(gcomp, name=cname, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': entered InitializeP1', esmf_logmsg_info)
    !
    ! -------------------------------------------------------------------- /
    ! Query mediator specific attributes
    !
    if (med_present) then
      call nuopc_compattributeget(gcomp, name="ScalarFieldName", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        flds_scalar_name = trim(cvalue)
        call esmf_logwrite(trim(cname)//': flds_scalar_name = '// &
             trim(flds_scalar_name), esmf_logmsg_info, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
      end if
 
      call nuopc_compattributeget(gcomp, name="ScalarFieldCount", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        flds_scalar_num = esmf_utilstring2int(cvalue, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
        if (verbosity.gt.0) then
          write(logmsg,*) flds_scalar_num
          call esmf_logwrite(trim(cname)//': flds_scalar_num = '// &
               trim(logmsg), esmf_logmsg_info, rc=rc)
          if (esmf_logfounderror(rc, passthru)) return
        end if
      end if
 
      call nuopc_compattributeget(gcomp, name="ScalarFieldIdxGridNX", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        flds_scalar_index_nx = esmf_utilstring2int(cvalue, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
        if (verbosity.gt.0) then
          write(logmsg,*) flds_scalar_index_nx
          call esmf_logwrite(trim(cname)//': flds_scalar_index_nx = '// &
               trim(logmsg), esmf_logmsg_info, rc=rc)
          if (esmf_logfounderror(rc, passthru)) return
        end if
      end if
 
      call nuopc_compattributeget(gcomp, name="ScalarFieldIdxGridNY", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        flds_scalar_index_ny = esmf_utilstring2int(cvalue, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
        if (verbosity.gt.0) then
          write(logmsg,*) flds_scalar_index_ny
          call esmf_logwrite(trim(cname)//': flds_scalar_index_ny = '// &
               trim(logmsg), esmf_logmsg_info, rc=rc)
          if (esmf_logfounderror(rc, passthru)) return
        end if
      end if
 
      call nuopc_compattributeget(gcomp, name="mask_value_water", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        maskvaluewater = esmf_utilstring2int(cvalue, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
        if (verbosity.gt.0) then
          write(logmsg,*) maskvaluewater
          call esmf_logwrite(trim(cname)//': mask_value_water = '// &
               trim(logmsg), esmf_logmsg_info, rc=rc)
          if (esmf_logfounderror(rc, passthru)) return
        end if
      end if
 
      call nuopc_compattributeget(gcomp, name="mask_value_land", &
           value=cvalue, ispresent=ispresent, isset=isset, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      if (ispresent .and. isset) then
        maskvalueland = esmf_utilstring2int(cvalue, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return
        if (verbosity.gt.0) then
          write(logmsg,*) maskvalueland
          call esmf_logwrite(trim(cname)//': mask_value_land = '// &
               trim(logmsg), esmf_logmsg_info, rc=rc)
          if (esmf_logfounderror(rc, passthru)) return
        end if
      end if
    end if
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Initialization necessary for driver
    !
    ! 1.a Set global flag indicating that model is an ESMF Component
    !
    is_esmf_component = .true.
    zerovalue = real(0,esmf_kind_rx)
    missingvalue = real(0,esmf_kind_rx)
    fillvalue = real(9.99e20,esmf_kind_rx)
    !
    !
    ! 1.b Get MPI environment from ESMF VM and set WW3 MPI related variables
    !
    call esmf_gridcompget(gcomp, vm=vm, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_vmget(vm, petcount=npet, localpet=lpet, &
         mpicommunicator=mpicomm, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    nmproc = npet
    improc = lpet + 1
    nmpscr = 1
    if ( improc .eq. nmpscr ) write (*,900)
    !
    ! 1.c Get background model info
    !
#if defined(COAMPS)
    call esmf_attributeget(gcomp, name="OcnBackground", &
         value=attstr, defaultvalue="none", &
         convention="COAMPS", purpose="General", rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    includeobg = trim(attstr).eq."model"
    call esmf_attributeget(gcomp, name="AtmBackground", &
         value=attstr, defaultvalue="none", &
         convention="COAMPS", purpose="General", rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    includeabg = trim(attstr).eq."model"
    call esmf_attributeget(gcomp, name="IceBackground", &
         value=attstr, defaultvalue="none", &
         convention="COAMPS", purpose="General", rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    includeibg = trim(attstr).eq."model"
    call esmf_attributeget(gcomp, name="MissingValue", &
         value=missingvalue, defaultvalue=real(0,esmf_kind_rx), &
         convention="COAMPS", purpose="General", rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
#else
    includeobg = .false.
    includeabg = .false.
    includeibg = .false.
#endif
    !
    ! 1.d Config input
    !
    call esmf_gridcompget(gcomp, configispresent=configispresent, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    if (configispresent) then
      call esmf_gridcompget(gcomp, config=config, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! working directory
      call esmf_configgetattribute(config, wrkdir, &
           label=trim(cname)//'_work_dir:', default='.', rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! I/O options
      call esmf_configgetattribute(config, ifname, &
           label=trim(cname)//'_input_file_name:', &
           default='ww3_multi.inp', rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_configgetattribute(config, lsep_ss, &
           label=trim(cname)//'_stdo_output_to_file:', default=.false., rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_configgetattribute(config, lsep_st, &
           label=trim(cname)//'_test_output_to_file:', default=.false., rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_configgetattribute(config, lsep_se, &
           label=trim(cname)//'_error_output_to_file:', default=.false., rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! export grid id
      call esmf_configgetattribute(config, expgridid, &
           label=trim(cname)//'_export_grid_id:', default=1, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! realize all export flag
      call esmf_configgetattribute(config, realizeallexport, &
           label=trim(cname)//'_realize_all_export:', default=.false., rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! grid mask convention
      call esmf_configgetattribute(config, maskvaluewater, &
           label='mask_value_water:', default=default_mask_water, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_configgetattribute(config, maskvalueland, &
           label='mask_value_land:', default=default_mask_land, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      ! z-level file
      call esmf_configgetattribute(config, zlfile, &
           label=trim(cname)//'_zlevel_exp_file:', default='none', rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
    endif
    !     preamb = trim(wrkdir)//'/'
    preamb = trim(preamb)//'/' !TODO: have separate paths for .inp, logs and data?
    !
    ! 1.e Set internal start/stop time from external start/stop time
    !
 
    call esmf_clockget(extclock, starttime=ttmp, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_clockget(extclock, currtime=cttmp, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    !
    ! Adjust internal start time to currTime in case of delayed start
    !
    if ( cttmp.gt.ttmp ) ttmp=cttmp
    call esmf_timeget(ttmp, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    stime(1) = 10000*yy + 100*mm + dd
    stime(2) = 10000*h  + 100*m  + s
 
    call esmf_clockget(extclock, stoptime=ttmp, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_timeget(ttmp, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    etime(1) = 10000*yy + 100*mm + dd
    etime(2) = 10000*h  + 100*m  + s
    !
    ! 1.f Identify available unit numbers
    ! Each ESMF_UtilIOUnitGet is followed by an OPEN statement for that
    ! unit so that subsequent ESMF_UtilIOUnitGet calls do not return the
    ! the same unit.  After getting all the available unit numbers, close
    ! the units since they will be opened within WMINIT.
    !
    call esmf_utiliounitget(idsi); open(unit=idsi, status='scratch');
    call esmf_utiliounitget(idso); open(unit=idso, status='scratch');
    call esmf_utiliounitget(idss); open(unit=idss, status='scratch');
    call esmf_utiliounitget(idst); open(unit=idst, status='scratch');
    call esmf_utiliounitget(idse); open(unit=idse, status='scratch');
    close(idsi); close(idso); close(idss); close(idst); close(idse);
    !
    ! 1.g Get merging option for regional coupling that domians does not
    ! overlap complately. This will blend the data coming from forcing with
    ! the data coming from coupling.
    !
    call nuopc_compattributeget(gcomp, name="merge_import", &
         value=attstr, ispresent=ispresent, isset=isset, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    if (ispresent .and. isset) then
      if (trim(attstr) .eq. '.true.') then
        merge_import = .true.
      end if
    end if
    if (verbosity.gt.0) then
      write(logmsg,'(l)') merge_import
      call esmf_logwrite(trim(cname)//import = '// &
           trim(logmsg), ESMF_LOGMSG_INFO, rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    end if
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Initialization of all wave models / grids
    !
    ! 2.a Call into WMINIT
    !
.not.    if ( lsep_ss ) idss = stdo
.not.    if ( lsep_st ) idst = stdo
.not.    if ( lsep_se ) idse = stdo
.eq.    if ( trim(ifname)'ww3_multi.nml' ) then
      call wminitnml ( idsi, idso, idss, idst, idse, trim(ifname), &
           mpicomm, preamb=preamb )
    else
      call wminit ( idsi, idso, idss, idst, idse, trim(ifname), &
           mpicomm, preamb=preamb )
    endif
    !
    ! 2.b Check consistency between internal timestep and external
    !     timestep (coupling interval)
    !
    call ESMF_ClockGet(extClock, timeStep=etstep, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 2.c Trap unsupported CPL input forcing settings
    !
.lt.    if ( any(inpmap0) ) then
.gt.      if ( nrgrd1 ) then
.eq.        if ( any(inpmap-999) ) then
          write (msg,'(a)') 'cpl input forcing defined on a '// &
               'native grid is not supported with multiple model grids'
.eq.          if ( improc  nmpscr ) write (idse,'(a)') trim(msg)
          call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
          rc = ESMF_FAILURE
          return
        endif
      endif
      allocate (cplmap(nrgrd,jfirst:8), stat=rc)
      if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
      jmod = minval(inpmap)
      cplmap = inpmap
.lt.      where ( inpmap0 ) cplmap = jmod
.ne.      if ( any(inpmapcplmap) ) then
        write (msg,'(a)') 'all cpl input forcing must be '// &
             'defined on the same grid'
.eq.        if ( improc  nmpscr ) write (idse,'(a)') trim(msg)
        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
        rc = ESMF_FAILURE
        return
      endif
      deallocate (cplmap, stat=rc)
      if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
import field list    ! 3.  Initialize 
    !
    istep_import: do istep = 1, 2
 
.eq.      if ( istep2 ) then
        allocate ( impFieldName(numImpFields), &
             impFieldStdName(numImpFields), &
             impFieldInitRqrd(numImpFields), &
             impFieldActive(numImpFields), &
             impField(numImpFields), &
             stat=rc )
        if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
        allocate ( mbgFieldName(numImpFields), &
             mbgFieldStdName(numImpFields), &
             mbgFieldActive(numImpFields), &
             mbgField(numImpFields), &
             bmskField(numImpFields), &
             stat=rc )
        if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
        if (merge_import) then
          allocate (mmskCreated(numImpFields))
          allocate (mmskField(numImpFields))
          allocate (mdtField(numImpFields))
          mmskCreated(:) = .false.
        end if
        impFieldActive(:) = .false.
      endif
      i = 0
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 1
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'seahgt'
        impFieldStdName(i)  = 'sea_surface_height_above_sea_level'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeObg
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 2
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'uucurr'
        impFieldStdName(i)  = 'surface_eastward_sea_water_velocity'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeObg
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 2
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'vvcurr'
        impFieldStdName(i)  = 'surface_northward_sea_water_velocity'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeObg
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 3
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'uutrue'
        impFieldStdName(i)  = 'eastward_wind_at_10m_height'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeAbg
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 3
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'vvtrue'
        impFieldStdName(i)  = 'northward_wind_at_10m_height'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeAbg
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        j = 4
.lt.        impFieldActive(i)   = any(inpmap(:,j)0)
        impFieldName(i)     = 'seaice'
        impFieldStdName(i)  = 'sea_ice_concentration'
        impFieldInitRqrd(i) = .true.
.and.        mbgFieldActive(i)   = impFieldActive(i)includeIbg
      endif
 
      numImpFields = i
import    enddo istep_
 
.not.    noActiveImpFields = all(impFieldActive)
 
    do i = 1,numImpFields
      mbgFieldName(i) = 'mbg_'//trim(impFieldName(i))
      mbgFieldStdName(i) = 'mbg_'//trim(impFieldStdName(i))
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! 4.  Initialize export field list
    !
    istep_export: do istep = 1, 2
 
.eq.      if ( istep2 ) then
        allocate ( expFieldName(numExpFields), &
             expFieldStdName(numExpFields), &
             expFieldDim(numExpFields), &
             expFieldActive(numExpFields), &
             expField(numExpFields), &
             stat=rc )
        if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
        expFieldActive(:) = .false.
      endif
      i = 0
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'charno'
        expFieldStdName(i) = 'wave_induced_charnock_parameter'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'z0rlen'
        expFieldStdName(i) = 'wave_z0_roughness_length'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'uscurr'
        expFieldStdName(i) = 'eastward_stokes_drift_current'
        expFieldDim(i)     = 3
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'vscurr'
        expFieldStdName(i) = 'northward_stokes_drift_current'
        expFieldDim(i)     = 3
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'x1pstk'
        expFieldStdName(i) = 'eastward_partitioned_stokes_drift_1'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'y1pstk'
        expFieldStdName(i) = 'northward_partitioned_stokes_drift_1'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'x2pstk'
        expFieldStdName(i) = 'eastward_partitioned_stokes_drift_2'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'y2pstk'
        expFieldStdName(i) = 'northward_partitioned_stokes_drift_2'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'x3pstk'
        expFieldStdName(i) = 'eastward_partitioned_stokes_drift_3'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'y3pstk'
        expFieldStdName(i) = 'northward_partitioned_stokes_drift_3'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wbcuru'
        expFieldStdName(i) = 'eastward_wave_bottom_current'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wbcurv'
        expFieldStdName(i) = 'northward_wave_bottom_current'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wbcurp'
        expFieldStdName(i) = 'wave_bottom_current_period'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wavsuu'
        expFieldStdName(i) = 'eastward_wave_radiation_stress'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wavsuv'
        expFieldStdName(i) = 'eastward_northward_wave_radiation_stress'
        expFieldDim(i)     = 2
      endif
 
      i = i + 1
.eq.      if ( istep2 ) then
        expFieldName(i)    = 'wavsvv'
        expFieldStdName(i) = 'northward_wave_radiation_stress'
        expFieldDim(i)     = 2
      endif
 
      if (med_present) then
        i = i + 1
.eq.        if ( istep2 ) then
          expFieldName(i)    = trim(flds_scalar_name)
          expFieldStdName(i) = trim(flds_scalar_name)
          expFieldDim(i)     = 1
        endif
      endif
 
      numExpFields = i
    enddo istep_export
 
.not.    noActiveExpFields = all(expFieldActive)
    !
    ! -------------------------------------------------------------------- /
import fields    ! 5.  Advertise 
    !
import fields    ! 5.a Advertise active 
    !
    n = 0
    do i = 1,numImpFields
.not.      if (impFieldActive(i)) cycle
      n = n + 1
      call NUOPC_Advertise(impState, &
           trim(impFieldStdName(i)), name=trim(impFieldName(i)), rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
.not.      if (mbgFieldActive(i)) cycle
      n = n + 1
      call NUOPC_Advertise(impState, &
           trim(mbgFieldStdName(i)), name=trim(mbgFieldName(i)), rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    enddo
    !
import fields    ! 5.b Report advertised 
    !
    write(msg,'(a,i0,a)') trim(cname)// &
         ': list of advertised import fields(',n,'):'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,a5,a,a10,a3,a)') trim(cname)// &
         ': ','index',' ','name',' ','standardname'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    n = 0
    do i = 1,numImpFields
.not.      if (impFieldActive(i)) cycle
      n = n + 1
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',n, &
           ' ',trim(impFieldName(i)),' ',trim(impFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
.not.      if (mbgFieldActive(i)) cycle
      n = n + 1
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',n, &
           ' ',trim(mbgFieldName(i)),' ',trim(mbgFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! 6.  Advertise export fields
    !
    ! 6.a Advertise all export fields
    !
    do i = 1,numExpFields
      call NUOPC_Advertise(expState, &
           trim(expFieldStdName(i)), name=trim(expFieldName(i)), rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    enddo
    !
    ! 6.b Report advertised export fields
    !
    write(msg,'(a,i0,a)') trim(cname)// &
         ': list of advertised export fields(',numExpFields,'):'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,a5,a,a10,a3,a)') trim(cname)// &
         ': ','index',' ','name',' ','standardname'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    do i = 1,numExpFields
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',i, &
           ' ',trim(expFieldName(i)),' ',trim(expFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving initializep1', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! Formats
    !
900 format (/15x,'     *** wavewatch iii multi-grid shell ***    '/ &
         15x,'================================================='/)
    !/
    !/ End of InitializeP1 ----------------------------------------------- /
    !/
  end subroutine InitializeP1
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "InitializeP3"
  !>
  !> @brief Initialize wave model (phase 3).
  !>
import and export states.  !> @details Realize fields in 
  !>
  !> @param gcomp    Gridded component.
Import state.  !> @param impState 
  !> @param expState Export state.
  !> @param extClock External clock.
  !> @param[out] rc  Return code.
  !>
  !> @author T. J. Campbell
  !> @author A. J. van der Westhuysen
  !> @date 09-Aug-2017
  !>
  subroutine InitializeP3 ( gcomp, impState, expState, extClock, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |     A. J. van der Westhuysen      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    09-Aug-2017 : Update 3D export field setup        ( version 6.03 )
    !/    28-Feb-2018 : Modifications for unstruc meshes    ( version 6.06 )
    !/
    !  1. Purpose :
    !
    !     Initialize wave model (phase 3)
import and export states.    !     * Realize fields in 
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
Import state    !       impState  Type   I/O 
    !       expState  Type   I/O Export state
    !       extClock  Type   I   External clock
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name            Type  Module   Description
    !     ----------------------------------------------------------------
    !      WMINIT          Subr. WMINITMD Wave model initialization
    !      WMINITNML       Subr. WMINITMD Wave model initialization
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    type(ESMF_State)    :: impState
    type(ESMF_State)    :: expState
    type(ESMF_Clock)    :: extClock
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=3
    real(8) :: wstime, wftime
    integer :: i1, i2, i3, i, n
    logical :: isConnected
    type(ESMF_DistGrid) :: distgrid
    type(ESMF_Grid) :: grid_scalar
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered initializep3', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
import fields    ! 1.  Realize active 
    !
import fields    ! 1.a Create ESMF grid for 
    !
.eq..or..eq.    if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
      write(msg,'(a)') trim(cname)// &
           ': creating import grid for reg/curvilinear mode'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
      call CreateImpGrid( gcomp, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
.eq.    elseif (GTYPEUNGTYPE) then
      write(msg,'(a)') trim(cname)// &
           ': creating import mesh for unstructured mode'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
      call CreateImpMesh( gcomp, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
import fields and realize    ! 1.b Create 
    !
    n = 0
    do i = 1,numImpFields
.not.      if (impFieldActive(i)) cycle
      n = n + 1
.eq..or..eq.      if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
        impField(i) = ESMF_FieldCreate( impGrid, impArraySpec2D, &
             totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
             staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
             name=trim(impFieldName(i)), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call FieldFill( impField(i), zeroValue, rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
.eq.      elseif (GTYPEUNGTYPE) then
        impField(i) = ESMF_FieldCreate( impMesh, &
             typekind=ESMF_TYPEKIND_RX, name=trim(impFieldName(i)), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call FieldFill( impField(i), zeroValue, rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      call NUOPC_Realize( impState, impField(i), rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
.eq..or..eq.      if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
        if (merge_import) then
          mmskField(i) = ESMF_FieldCreate( impGrid, impArraySpec2D, &
               totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
               staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
               name='mmsk_'//trim(impFieldName(i)), rc=rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          call FieldFill( mmskField(i), zeroValue, rc=rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          mdtField(i) = ESMF_FieldCreate( impGrid, impArraySpec2D, &
               totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
               staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
               name='mdt_'//trim(impFieldName(i)), rc=rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          call FieldFill( mdtField(i), zeroValue, rc=rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        end if
.not.        if (mbgFieldActive(i)) cycle
        n = n + 1
        mbgField(i) = ESMF_FieldCreate( impGrid, impArraySpec2D, &
             totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
             staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
             name=trim(mbgFieldName(i)), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call FieldFill( mbgField(i), zeroValue, rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call NUOPC_Realize( impState, mbgField(i), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        bmskField(i) = ESMF_FieldCreate( impGrid, impArraySpec2D, &
             totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
             staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
             name='bmsk_'//trim(impFieldName(i)), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
    enddo
    !
import fields    ! 1.c Report realized 
    !
    write(msg,'(a,i0,a)') trim(cname)// &
         ': list of realized import fields(',n,'):'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,a5,a,a10,a3,a)') trim(cname)// &
         ': ','index',' ','name',' ','standardname'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    n = 0
    do i = 1,numImpFields
.not.      if (impFieldActive(i)) cycle
      n = n + 1
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',n, &
           ' ',trim(impFieldName(i)),' ',trim(impFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
.not.      if (mbgFieldActive(i)) cycle
      n = n + 1
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',n, &
           ' ',trim(mbgFieldName(i)),' ',trim(mbgFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Realize active export fields
    !
    ! 2.a Set connected export fields as active and remove unconnected
    !     If realizeAllExport, then set all fields as active and realize.
    !
    do i = 1,numExpFields
      isConnected = NUOPC_IsConnected(expState, &
           expFieldName(i), rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
.or.      expFieldActive(i) = isConnected  realizeAllExport
      if (expFieldActive(i)) noActiveExpFields = .false.
.not.      if (expFieldActive(i)) then
        call ESMF_StateRemove(expState, (/expFieldName(i)/), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        write(msg,fmt="(a,l)") trim(cname)//': '//trim(expFieldName(i)), expFieldActive(i)
        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
      endif
    enddo
    !
    ! 2.b Create ESMF grid for export fields
    !
.eq..or..eq.    if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
      write(msg,'(a)') trim(cname)// &
           ': creating export grid for reg/curvilinear mode'
      call CreateExpGrid( gcomp, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
.eq.    elseif (GTYPEUNGTYPE) then
      write(msg,'(a)') trim(cname)// &
           ': creating export mesh for unstructured mode'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
      call CreateExpMesh( gcomp, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! 2.c Create active export fields and realize
    !
    n = 0
    do i = 1,numExpFields
.not.      if (expFieldActive(i)) cycle
      n = n + 1
.eq..or..eq.      if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
        if (trim(expFieldName(i)) == trim(flds_scalar_name)) then
          distgrid = ESMF_DistGridCreate(minIndex=(/1/), maxIndex=(/1/), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          grid_scalar = ESMF_GridCreate(distgrid, rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          expField(i) = ESMF_FieldCreate(grid_scalar, typekind=ESMF_TYPEKIND_R8, &
               name=trim(expFieldName(i)), ungriddedLBound=(/1/), &
               ungriddedUBound=(/flds_scalar_num/), gridToFieldMap=(/2/), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        else
.eq.          if ( expFieldDim(i)3 ) then
            expField(i) = ESMF_FieldCreate( expGrid, expArraySpec3D, &
                 totalLWidth=expHaloLWidth, totalUWidth=expHaloUWidth, &
                 gridToFieldMap=(/2,3/), ungriddedLBound=(/1/), ungriddedUBound=(/nz/), &
                 staggerLoc=expStaggerLoc, name=trim(expFieldName(i)), rc=rc )
            if (ESMF_LogFoundError(rc, PASSTHRU)) return
          else
            expField(i) = ESMF_FieldCreate( expGrid, expArraySpec2D, &
                 totalLWidth=expHaloLWidth, totalUWidth=expHaloUWidth, &
                 staggerLoc=expStaggerLoc, name=trim(expFieldName(i)), rc=rc )
            if (ESMF_LogFoundError(rc, PASSTHRU)) return
          endif
          call FieldFill( expField(i), zeroValue, rc=rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        endif
.eq.      elseif (GTYPEUNGTYPE) then
        expField(i) = ESMF_FieldCreate( expMesh, &
             typekind=ESMF_TYPEKIND_RX, name=trim(expFieldName(i)), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call FieldFill( expField(i), zeroValue, rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      call NUOPC_Realize( expState, expField(i), rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    enddo
    !
    ! 2.d Report realized export fields
    !
    write(msg,'(a,i0,a)') trim(cname)// &
         ': list of realized export fields(',n,'):'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,a5,a,a10,a3,a)') trim(cname)// &
         ': ','index',' ','name',' ','standardname'
    call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    n = 0
    do i = 1,numExpFields
.not.      if (expFieldActive(i)) cycle
      n = n + 1
      write(msg,'(a,i5,a,a10,a3,a)') trim(cname)//': ',n, &
           ' ',trim(expFieldName(i)),' ',trim(expFieldStdName(i))
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    enddo
    !
    ! 2.e Set W3OUTG flags needed for calculating export fields
    !
#ifdef USE_W3OUTG_FOR_EXPORT
    call w3seto ( expGridID, mdse, mdst )
 
    i1 = FieldIndex( expFieldName, 'uscurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'vscurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
         expFieldActive(i2) ) then
      flogr2(6,8) = .true.  !Spectrum of surface Stokes drift
.le.      if ( us3df(1)  0 ) then
        msg = trim(cname)//': stokes drift export using w3outg'// &
             ' requires setting us3d=1 (ww3_grid.inp: outs namelist)'
        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
        rc = ESMF_FAILURE
        return
      endif
    endif
 
    i1 = FieldIndex( expFieldName, 'wbcuru', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'wbcurv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i3 = FieldIndex( expFieldName, 'wbcurp', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
.and.         expFieldActive(i2)  &
         expFieldActive(i3) ) then
      flogr2(7,1) = .true.  !Near bottom rms amplitudes
      flogr2(7,2) = .true.  !Near bottom rms velocities
    endif
 
    i1 = FieldIndex( expFieldName, 'wavsuu', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'wavsuv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i3 = FieldIndex( expFieldName, 'wavsvv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
.and.         expFieldActive(i2)  &
         expFieldActive(i3) ) then
      flogr2(6,1) = .true.  !Radiation stresses
    endif
#endif
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving initializep3', ESMF_LOGMSG_INFO)
    !/
    !/ End of InitializeP3 ----------------------------------------------- /
    !/
  end subroutine InitializeP3
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "Finalize"
  !>
  !> @brief Finalize wave model.
  !>
  !> @param gcomp   Gridded component.
  !> @param[out] rc Return code.
  !>
  !> @author T. J. Campbell  @date 09-Aug-2017
  !>
  subroutine Finalize ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    09-Aug-2017 : Add clean up of local allocations   ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Finalize wave model
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      WMFINL    Subr. WMFINLMD Wave model finalization
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=6
    real(8) :: wstime, wftime
    integer :: i
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered finalize', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Finalize the wave model
    !
    call wmfinl
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Clean up ESMF data structures
    !
Import field and grid stuff    ! 2.a 
    !
.not.    if ( noActiveImpFields ) then
 
      do i = 1,numImpFields
.not.        if (impFieldActive(i)) cycle
        call ESMF_FieldDestroy(impField(i), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        if (merge_import) then
          call ESMF_FieldDestroy(mdtField(i), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          call ESMF_FieldDestroy(mmskField(i), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        end if
.not.        if (mbgFieldActive(i)) cycle
        call ESMF_FieldDestroy(mbgField(i), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call ESMF_FieldDestroy(bmskField(i), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      enddo
 
.eq..or..eq.      if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
        call ESMF_FieldHaloRelease(impHaloRH, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
        call ESMF_GridDestroy(impGrid, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
.eq.      elseif (GTYPEUNGTYPE) then
        !AW          call ESMF_GridDestroy(impMesh, rc=rc)
        !AW          if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
 
    endif
 
    deallocate (impFieldName, &
         impFieldStdName, &
         impFieldInitRqrd, &
         impFieldActive, &
         impField, &
         stat=rc)
    if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
 
    deallocate (mbgFieldName, &
         mbgFieldStdName, &
         mbgFieldActive, &
         mbgField, &
         bmskField, &
         stat=rc)
    if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
 
    if (merge_import) then
      deallocate(mmskCreated, &
           mmskField, &
           mdtField, &
           stat=rc)
      if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    end if
    !
    ! 2.b Export field and grid stuff
    !
.not.    if ( noActiveExpFields ) then
 
      do i = 1,numExpFields
.not.        if (expFieldActive(i)) cycle
        call ESMF_FieldDestroy(expField(i), rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      enddo
 
.eq..or..eq.      if ( (GTYPERLGTYPE)(GTYPECLGTYPE) ) then
        call ESMF_FieldHaloRelease(expHaloRH, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
        call ESMF_GridDestroy(expGrid, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
.eq.      elseif (GTYPEUNGTYPE) then
        !AW          call ESMF_GridDestroy(expMesh, rc=rc)
        !AW          if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
 
    endif
 
    deallocate (expFieldName, &
         expFieldStdName, &
         expFieldDim, &
         expFieldActive, &
         expField, &
         stat=rc)
    if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    !
    ! 2.c Native field and grid stuff
    !
.not.    if ( noActiveExpFields ) then
 
      call ESMF_FieldRedistRelease(n2eRH, rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
      call ESMF_GridDestroy(natGrid, rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
    endif
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Clean up locally allocated data structures
    !
    if (allocated(zl)) then
      deallocate (zl, stat=rc)
      if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
    call PrintTimers(trim(cname), wtnam, wtcnt, wtime)
    rc = ESMF_SUCCESS
.eq.    if ( improc  nmpscr ) write (*,999)
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving finalize', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! Formats
    !
999 format(//'  end of program '/                                   &
         ' ========================================'/           &
         '          wavewatch iii multi-grid shell '/)
    !/
    !/ End of Finalize --------------------------------------------------- /
    !/
  end subroutine Finalize
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "DataInitialize"
  !>
  !> @brief  Initialize wave model export data
  !>
  !> @param gcomp   Gridded component.
  !> @param[out] rc Return code.
  !>
  !> @author T. J. Campbell  @date 20-Jan-2017
  !>
  subroutine DataInitialize ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Initialize wave model export data
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name            Type  Module   Description
    !     ----------------------------------------------------------------
Import       subr. wmesmfmd wave model get import fields    !      Get
    !      SetExport       Subr. WMESMFMD Wave model set export fields
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=4
    real(8) :: wstime, wftime
    type(ESMF_Clock) :: clock
    type(ESMF_Time) :: currTime
    logical :: fldUpdated, allUpdated
    integer :: i, imod
    logical :: local
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered datainitialize', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
import fields show correct time stamp    ! 1.  Check that required 
    !
    if (med_present) then
      allUpdated = .true.
    else
      call ESMF_GridCompGet(gcomp, clock=clock, rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call ESMF_ClockGet(clock, currTime=currTime, rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
      allUpdated = .true.
      do i = 1,numImpFields
.not.        if (impFieldActive(i)) cycle
        if (impFieldInitRqrd(i)) then
          fldUpdated = NUOPC_IsAtTime(impField(i), currTime, rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          if (fldUpdated) then
            write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
                 trim(impFieldName(i)),': inter-model data dependency: ', &
                 'satisfied'
          else
            allUpdated = .false.
            write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
                 trim(impFieldName(i)),': inter-model data dependency: ', &
                 'not satisfied'
          endif
        else
          write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
               trim(impFieldName(i)),': inter-model data dependency: ', &
               'not  required'
        endif
.gt.        if (verbosity0) call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
.eq.        if (improcnmpscr) write(*,'(a)') trim(msg)
        ! background
.not.        if (mbgFieldActive(i)) cycle
        if (impFieldInitRqrd(i)) then
          fldUpdated = NUOPC_IsAtTime(mbgField(i), currTime, rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          if (fldUpdated) then
            write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
                 trim(mbgFieldName(i)),': inter-model data dependency: ', &
                 'satisfied'
          else
            allUpdated = .false.
            write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
                 trim(mbgFieldName(i)),': inter-model data dependency: ', &
                 'not satisfied'
          endif
        else
          write(msg,'(a,a10,a,a13)') trim(cname)//': ', &
               trim(mbgFieldName(i)),': inter-model data dependency: ', &
               'not  required'
        endif
.gt.        if (verbosity0) call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
.eq.        if (improcnmpscr) write(*,'(a)') trim(msg)
      enddo
    endif
    !
import dependencies are satisfied, then return    ! If not all 
    !
.not.    if (allUpdated) goto 1
    !
    ! -------------------------------------------------------------------- /
import dependencies are satisfied, so finish initialization    ! 2.  All 
    !
import dependencies are satisfied    ! 2.a Report all 
    !
    write(msg,'(a)') trim(cname)// &
         ': all inter-model data dependencies satisfied'
.gt.    if (verbosity0) call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
.eq.    if (improcnmpscr) write(*,'(a)') trim(msg)
    !
import field    ! 2.b Setup background blending mask for each 
    !
    do i = 1,numImpFields
.not.      if (impFieldActive(i)) cycle
.not.      if (mbgFieldActive(i)) cycle
      call SetupImpBmsk(bmskField(i), impField(i), missingValue, rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    enddo
    !
import fields    ! 2.c Get 
    !
    call GetImport(gcomp, rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 2.d Finish initialization (compute initial state), if not restart
    !
    do imod = 1,nrgrd
      call w3setg ( imod, mdse, mdst )
      call w3setw ( imod, mdse, mdst )
      call w3seta ( imod, mdse, mdst )
      call w3seti ( imod, mdse, mdst )
      call w3seto ( imod, mdse, mdst )
      call wmsetm ( imod, mdse, mdst )
.gt..and..le.      local = iaproc  0  iaproc  naproc
.and..and.      if ( local  flcold  fliwnd ) call w3uini( va )
    enddo
    !
    ! 2.e Set export fields
    !
    call SetExport(gcomp, rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 2.f Set Updated Field Attribute to "true", indicating to the
    !     generic code to set the timestamp for these fields
    !
    do i = 1,numExpFields
.not.      if (expFieldActive(i)) cycle
      call NUOPC_SetAttribute(expField(i), name="Updated", &
           value="true", rc=rc)
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    enddo
    !
    ! 2.g Set InitializeDataComplete Attribute to "true", indicating to the
    !     generic code that all inter-model data dependencies are satisfied
    !
    call NUOPC_CompAttributeSet(gcomp, name="InitializeDataComplete", &
         value="true", rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
1   rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving datainitialize', ESMF_LOGMSG_INFO)
    !/
    !/ End of DataInitialize --------------------------------------------- /
    !/
  end subroutine DataInitialize
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "ModelAdvance"
  !>
  !> @brief Advance wave model in time.
  !>
  !> @param      gcomp Gridded component.
  !> @param[out] rc    Return code.
  !>
  !> @author T. J. Campbell   @date 20-Jan-2017
  !>
  subroutine ModelAdvance ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Advance wave model in time
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
Import subr. wmesmfmd wave model get import fields    !      Get
    !      SetExport Subr. WMESMFMD Wave model set export fields
    !      WMWAVE    Subr. WMWAVEMD Wave model run
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=5
    real(8) :: wstime, wftime
    integer :: i, stat, imod, tcur(2)
    integer, allocatable :: tend(:,:)
    integer(ESMF_KIND_I4) :: yy,mm,dd,h,m,s
    type(ESMF_Clock) :: clock
    type(ESMF_Time) :: currTime, stopTime
    real :: delt
    logical :: lerr
 
    type(ESMF_Time) :: startTime
    type(ESMF_TimeInterval)    :: timeStep
    character(len=256)         :: msgString
 
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered modeladvance', ESMF_LOGMSG_INFO)
    if(profile_memory) call ESMF_VMLogMemInfo('entering ww3 '// &
         'model_advance: ')
 
    allocate (tend(2,nrgrd), stat=rc)
    if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Advance model to requested end time
    !
    ! 1.a Get component clock
    !
    call ESMF_GridCompGet(gcomp, clock=clock, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 1.b Report
    !
.eq.    if ( improc  nmpscr ) then
      write(*,'(///)')
      call ESMF_ClockPrint(clock, options="currTime", &
           preString="-->Advancing "//TRIM(cname)//" from: ")
      call ESMF_ClockPrint(clock, options="stopTime", &
           preString="-----------------> to: ")
    endif
 
    if (profile_memory) then
      call ESMF_ClockPrint(clock, options="currTime", &
           preString="------>Advancing WAV from: ", &
           unit=msgString, rc=rc)
      call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO)
 
      call ESMF_ClockGet(clock, startTime=startTime, &
           currTime=currTime, &
           timeStep=timeStep, rc=rc)
 
      call ESMF_TimePrint(currTime + timeStep, &
           preString="--------------------------------> to: ", &
           unit=msgString, rc=rc)
      call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO)
    endif
    !
    ! 1.c Check internal current time with component current time
    !
    call ESMF_ClockGet(clock, currTime=currTime, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_TimeGet(currTime, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    lerr=.false.
    do imod = 1,nrgrd
      tcur(1) = 10000*yy + 100*mm + dd
      tcur(2) = 10000*h  + 100*m  + s
      call w3setw ( imod, mdse, mdst )
      delt = dsec21 ( time, tcur )
.gt.      if ( abs(delt)0 ) then
        lerr=.true.
        write(msg,'(a,i2,a,2(a,i8,a,i8,a))') &
             'wave model grid ',imod,': ', &
             'internal time (',time(1),'.',time(2),') /= ', &
             'component time (',tcur(1),'.',tcur(2),')'
        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
      endif
    enddo
    if (lerr) then
      rc = ESMF_FAILURE
      return
    endif
    !
    ! 1.d Set end time of this advance
    !
    call ESMF_ClockGet(clock, stopTime=stopTime, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_TimeGet(stopTime, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    do imod = 1,nrgrd
      tend(1,imod) = 10000*yy + 100*mm + dd
      tend(2,imod) = 10000*h  + 100*m  + s
    enddo
    !
import fields    ! 1.e Get 
    !
    call GetImport(gcomp, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 1.f Advance model
    !
    if(profile_memory) call ESMF_VMLogMemInfo("Entering WW3 Run : ")
    call wmwave ( tend )
    if(profile_memory) call ESMF_VMLogMemInfo("Entering WW3 Run : ")
    !
    ! 1.g Set export fields
    !
    call SetExport(gcomp, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
    deallocate (tend, stat=rc)
    if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving modeladvance', ESMF_LOGMSG_INFO)
    if(profile_memory) call ESMF_VMLogMemInfo('leaving ww3 '// &
         'model_advance: ')
    !/
    !/ End of ModelAdvance ----------------------------------------------- /
    !/
  end subroutine ModelAdvance
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "GetImport"
  !>
import fields and put in internal data structures.  !> @brief Get 
  !>
  !> @param      gcomp Gridded component.
  !> @param[out] rc    Return code.
  !>
  !> @author T. J. Campbell  @date 20-Jan-2017
  !>
Import ( gcomp, rc )  subroutine Get
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
import fields and put in internal data structures    !     Get 
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
#ifdef W3_MPI
    USE WMMDATMD, ONLY: IMPROC
#endif
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    !AW ---TEST-TEST-TEST---------------------------
    character(500) :: msg
    integer :: k
    !AW ---TEST-TEST-TEST---------------------------
    integer, parameter :: iwt=7
    real(8) :: wstime, wftime
    integer :: i1, i2, i3, j, imod, jmod
    logical, save :: firstCall = .true.
    integer :: tcur(2), tend(2)
    integer(ESMF_KIND_I4) :: yy,mm,dd,h,m,s
    type(ESMF_Clock) :: clock
    type(ESMF_Time) :: currTime, stopTime
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_GETIMPORT)
    type(ESMF_State) :: dumpState
    integer, save :: timeSlice = 1
#endif
    real(ESMF_KIND_RX), pointer :: rptr(:,:)
    integer :: iy, ix
    integer :: elb(2), eub(2)
    character(len=3) :: fieldName
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    if ( noActiveImpFields ) return
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered getimport', ESMF_LOGMSG_INFO)
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_GETIMPORT)
    call NUOPC_ModelGet(gcomp, importState=dumpState, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call NUOPC_Write(dumpState, overwrite=.true., &
         fileNamePrefix="field_"//trim(cname)//"_import1_", &
         timeslice=timeSlice, relaxedFlag=.true., rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
#endif
    !
    ! -------------------------------------------------------------------- /
    ! Set time stamps using currTime and stopTime
    !
    call ESMF_GridCompGet(gcomp, clock=clock, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_ClockGet(clock, currTime=currTime, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_TimeGet(currTime, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    tcur(1) = 10000*yy + 100*mm + dd
    tcur(2) = 10000*h  + 100*m  + s
    call ESMF_ClockGet(clock, stopTime=stopTime, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_TimeGet(stopTime, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    tend(1) = 10000*yy + 100*mm + dd
    tend(2) = 10000*h  + 100*m  + s
    !
    ! -------------------------------------------------------------------- /
    ! Water levels
    !
    j = 1
    i1 = FieldIndex( impFieldName, 'seahgt', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = i1
    if ( impFieldActive(i1) ) then
      call w3setg ( impGridID, mdse, mdst )
      call w3seti ( impGridID, mdse, mdst )
      if (firstCall) then
        tln = tcur
      else
        tln = tend
      endif
      tfn(:,j) = tln
      if ( mbgFieldActive(i1) ) then
        call BlendImpField( impField(i1), mbgField(i1), bmskField(i1), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      call FieldGather( impField(i1), nx, ny, wlev, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do imod = 1,nrgrd
        call w3setg ( imod, mdse, mdst )
        call w3setw ( imod, mdse, mdst )
        call w3seti ( imod, mdse, mdst )
        call w3seto ( imod, mdse, mdst )
        call wmsetm ( imod, mdse, mdst )
#ifdef W3_MPI
.eq.        if ( mpi_comm_grd  mpi_comm_null ) cycle
#endif
        jmod = inpmap(imod,j)
.lt..and..ne.        if ( jmod0  jmod-999 ) then
          call wmupd2( imod, j, jmod, rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        endif
      enddo
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Currents
    !
    j = 2
    i1 = FieldIndex( impFieldName, 'uucurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( impFieldName, 'vvcurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impFieldActive(i1) ) then
      call w3setg ( impGridID, mdse, mdst )
      call w3seti ( impGridID, mdse, mdst )
      if (firstCall) then
        tcn = tcur
      else
        tc0 = tcn
        cx0 = cxn
        cy0 = cyn
        tcn = tend
      endif
      tfn(:,j) = tcn
      if ( mbgFieldActive(i1) ) then
        call BlendImpField( impField(i1), mbgField(i1), bmskField(i1), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call BlendImpField( impField(i2), mbgField(i2), bmskField(i2), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      call FieldGather( impField(i1), nx, ny, cxn, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call FieldGather( impField(i2), nx, ny, cyn, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      if (firstCall) then
        tc0 = tcn
        cx0 = cxn
        cy0 = cyn
      endif
      do imod = 1,nrgrd
        call w3setg ( imod, mdse, mdst )
        call w3setw ( imod, mdse, mdst )
        call w3seti ( imod, mdse, mdst )
        call wmsetm ( imod, mdse, mdst )
#ifdef W3_MPI
.eq.        if ( mpi_comm_grd  mpi_comm_null ) cycle
#endif
        jmod = inpmap(imod,j)
.lt..and..ne.        if ( jmod0  jmod-999 ) then
          call wmupd2( imod, j, jmod, rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        endif
      enddo
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Winds
    !
    j = 3
    i1 = FieldIndex( impFieldName, 'uutrue', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( impFieldName, 'vvtrue', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impFieldActive(i1) ) then
      call w3setg ( impGridID, mdse, mdst )
      call w3seti ( impGridID, mdse, mdst )
 
      if (firstCall) then
        twn = tcur
      else
        tw0 = twn
        wx0 = wxn
        wy0 = wyn
        twn = tend
      endif
      tfn(:,j) = twn
      if ( mbgFieldActive(i1) ) then
        call BlendImpField( impField(i1), mbgField(i1), bmskField(i1), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call BlendImpField( impField(i2), mbgField(i2), bmskField(i2), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      if (merge_import) then
        ! read wave input
        fieldName = 'wnd'
        call ReadFromFile(fieldName, mdtField(i1), mdtField(i2), tcur, tend, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        ! create merge mask
.not.        if ( firstCall) then
          call SetupImpMmsk(mmskField(i1), impField(i1), fillValue, mmskCreated(i1), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
          call SetupImpMmsk(mmskField(i2), impField(i2), fillValue, mmskCreated(i2), rc=rc)
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        end if
        ! blend data, mask is all zero initially (use all data)
        call BlendImpField( impField(i1), mdtField(i1), mmskField(i1), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
        call BlendImpField( impField(i2), mdtField(i2), mmskField(i2), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      end if
      call FieldGather( impField(i1), nx, ny, wxn, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call FieldGather( impField(i2), nx, ny, wyn, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      if (firstCall) then
        tw0 = twn
        wx0 = wxn
        wy0 = wyn
#ifdef W3_WRST
        ! The WRST switch saves the values of wind in the
        ! restart file and then uses the wind for the first
        ! time step here.  This is needed when coupling with
        ! an atm model that does not have 10m wind speeds at
        ! initialization.  If there is no restart, wind is zero
        wxn = WXNwrst !replace with values from restart
        wyn = WYNwrst
        wx0 = WXNwrst
        wy0 = WYNwrst
        do imod = 1,nrgrd
          call w3setg ( imod, mdse, mdst )
          call w3setw ( imod, mdse, mdst )
          call w3seti ( imod, mdse, mdst )
          call wmsetm ( imod, mdse, mdst )
#ifdef W3_MPI
.eq.          if ( mpi_comm_grd  mpi_comm_null ) cycle
#endif
          INPUTS(IMOD)%TW0(:) = INPUTS(impGridID)%TW0(:)
          INPUTS(IMOD)%TFN(:,3) = INPUTS(impGridID)%TFN(:,3)
          wxn = WXNwrst !replace with values from restart
          wyn = WYNwrst
          wx0 = WXNwrst
          wy0 = WYNwrst
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        enddo
#endif
      endif
 
#ifdef W3_WRST
.le.      if ( ((twn(1)-tw0(1))*1000000+((twn(2)-tw0(2))))  0  ) then
        !If the time of the field is still initial time, replace
        !with restart field
        wxn = WXNwrst !replace with values from restart
        wyn = WYNwrst
      else    !twn>tw0
#endif
        do imod = 1,nrgrd
          call w3setg ( imod, mdse, mdst )
          call w3setw ( imod, mdse, mdst )
          call w3seti ( imod, mdse, mdst )
          call wmsetm ( imod, mdse, mdst )
#ifdef W3_MPI
.eq.          if ( mpi_comm_grd  mpi_comm_null ) cycle
#endif
          jmod = inpmap(imod,j)
.lt..and..ne.          if ( jmod0  jmod-999 ) then
            call wmupd2( imod, j, jmod, rc )
            if (ESMF_LogFoundError(rc, PASSTHRU)) return
          endif
        enddo
#ifdef W3_WRST
.le.      endif  !if ( twn-tw0  0  )
#endif
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Sea ice concentration
    !
    j = 4
    i1 = FieldIndex( impFieldName, 'seaice', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = i1
    if ( impFieldActive(i1) ) then
      call w3setg ( impGridID, mdse, mdst )
      call w3seti ( impGridID, mdse, mdst )
      if (firstCall) then
        tin = tcur
      else
        tin = tend
      endif
      tfn(:,j) = tin
      if ( mbgFieldActive(i1) ) then
        call BlendImpField( impField(i1), mbgField(i1), bmskField(i1), rc=rc )
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
      endif
      call FieldGather( impField(i1), nx, ny, icei, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do imod = 1,nrgrd
        call w3setg ( imod, mdse, mdst )
        call w3setw ( imod, mdse, mdst )
        call w3seti ( imod, mdse, mdst )
        call wmsetm ( imod, mdse, mdst )
#ifdef W3_MPI
.eq.        if ( mpi_comm_grd  mpi_comm_null ) cycle
#endif
        jmod = inpmap(imod,j)
.lt..and..ne.        if ( jmod0  jmod-999 ) then
          call wmupd2( imod, j, jmod, rc )
          if (ESMF_LogFoundError(rc, PASSTHRU)) return
        endif
      enddo
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_GETIMPORT)
    call NUOPC_Write(dumpState, overwrite=.true., &
         fileNamePrefix="field_"//trim(cname)//"_import2_", &
         timeslice=timeSlice, relaxedFlag=.true., rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    timeSlice = timeSlice + 1
#endif
    firstCall = .false.
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving getimport', ESMF_LOGMSG_INFO)
    !/
Import -------------------------------------------------- /    !/ End of Get
    !/
Import  end subroutine Get
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "SetExport"
  !>
  !> @brief Set export fields from internal data structures.
  !>
  !> @param      gcomp Gridded component
  !> @param[out] rc    Return code
  !>
  !> @author T. J. Campbell  @date 09-Aug-2017
  !>
  subroutine SetExport ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    09-Aug-2017 : Add ocean forcing export fields     ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Set export fields from internal data structures
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer, parameter :: iwt=8
    real(8) :: wstime, wftime
    integer :: i1, i2, i3, i4, i5, i6
    logical :: flpart = .false., floutg = .false., floutg2 = .true.
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETEXPORT)
    type(ESMF_State) :: dumpState
    integer, save :: timeSlice = 1
#endif
    real(ESMF_KIND_R8), pointer :: farrayptr(:,:)
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    if ( noActiveExpFields ) return
    call ESMF_VMWtime(wstime)
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered setexport', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! Setup
    !
    call w3setg ( expGridID, mdse, mdst )
    call w3setw ( expGridID, mdse, mdst )
    call w3seta ( expGridID, mdse, mdst )
    call w3seti ( expGridID, mdse, mdst )
    call w3seto ( expGridID, mdse, mdst )
    call wmsetm ( expGridID, mdse, mdst )
#ifdef USE_W3OUTG_FOR_EXPORT
    if ( natGridIsLocal ) call w3outg( va, flpart, floutg, floutg2 )
#endif
    !
    ! -------------------------------------------------------------------- /
    ! Charnock
    !
    i1 = FieldIndex( expFieldName, 'charno', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expFieldActive(i1) ) then
      call CalcCharnk( expField(i1), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Surface Roughness
    !
    i1 = FieldIndex( expFieldName, 'z0rlen', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expFieldActive(i1) ) then
      call CalcRoughl( expField(i1), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Stokes Drift 3D
    !
    i1 = FieldIndex( expFieldName, 'uscurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'vscurr', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
         expFieldActive(i2) ) then
      call CalcStokes3D( va, expField(i1), expField(i2), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Partitioned Stokes Drift 3 2D fields
    !
    i1 = FieldIndex( expFieldName, 'x1pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'y1pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i3 = FieldIndex( expFieldName, 'x2pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i4 = FieldIndex( expFieldName, 'y2pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i5 = FieldIndex( expFieldName, 'x3pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i6 = FieldIndex( expFieldName, 'y3pstk', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
.and.         expFieldActive(i2)  &
.and.         expFieldActive(i3)  &
.and.         expFieldActive(i4)  &
.and.         expFieldActive(i5)  &
         expFieldActive(i6)  ) then
      call CalcPStokes( va, expField(i1), expField(i2), expField(i3), &
           expField(i4), expField(i5), expField(i6), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Bottom Currents
    !
    i1 = FieldIndex( expFieldName, 'wbcuru', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'wbcurv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i3 = FieldIndex( expFieldName, 'wbcurp', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
.and.         expFieldActive(i2)  &
         expFieldActive(i3) ) then
      call CalcBotcur( va, expField(i1), expField(i2), expField(i3), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Radiation stresses 2D
    !
    i1 = FieldIndex( expFieldName, 'wavsuu', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i2 = FieldIndex( expFieldName, 'wavsuv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    i3 = FieldIndex( expFieldName, 'wavsvv', rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.    if ( expFieldActive(i1)  &
.and.         expFieldActive(i2)  &
         expFieldActive(i3) ) then
      call CalcRadstr2D( va, expField(i1), expField(i2), expField(i3), rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! cpl_scalars - grid sizes
    !
    if (med_present) then
      i1 = FieldIndex( expFieldName, trim(flds_scalar_name), rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      if ( expFieldActive(i1) ) then
        call ESMF_FieldGet(expField(i1), farrayPtr=farrayptr, rc=rc)
        if (ESMF_LogFoundError(rc, PASSTHRU)) return
.and.        if (flds_scalar_index_nx > 0  flds_scalar_index_nx < flds_scalar_num) then
          farrayptr(flds_scalar_index_nx,1) = dble(nx)
        endif
.and.        if (flds_scalar_index_ny > 0  flds_scalar_index_ny < flds_scalar_num) then
          farrayptr(flds_scalar_index_ny,1) = dble(ny)
        endif
      endif
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETEXPORT)
    call NUOPC_ModelGet(gcomp, exportState=dumpState, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call NUOPC_Write(dumpState, overwrite=.true., &
         fileNamePrefix="field_"//trim(cname)//"_export_", &
         timeslice=timeSlice, relaxedFlag=.true., rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    timeSlice = timeSlice + 1
#endif
    rc = ESMF_SUCCESS
    call ESMF_VMWtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving setexport', ESMF_LOGMSG_INFO)
    !/
    !/ End of SetExport -------------------------------------------------- /
    !/
  end subroutine SetExport
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CreateImpGrid"
  !>
import fields.  !> @brief Create ESMF grid for 
  !>
  !> @param      gcomp Gridded component
  !> @param[out] rc    Return code
  !>
  !> @author T. J. Campbell  @date 20-Jan-2017
  !>
  subroutine CreateImpGrid ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
import fields    !     Create ESMF grid for 
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    type(ESMF_VM) :: vm
    character(ESMF_MAXSTR) :: cname
    integer :: nproc, nxproc, nyproc
    integer, parameter :: lde = 0
    integer :: lpet, ldecnt
    integer :: ix, iy, isea, jsea, irec, ubx, uby
    integer :: elb(2), eub(2), elbc(2), eubc(2)
    integer :: tlb(2), tub(2)
    integer(ESMF_KIND_I4), pointer :: iptr(:,:)
    real(ESMF_KIND_RX), pointer :: rptrx(:,:), rptry(:,:)
    real(ESMF_KIND_RX), pointer :: rptr(:,:)
    real(ESMF_KIND_R8), allocatable :: xgrd_center(:)
    real(ESMF_KIND_R8), allocatable :: ygrd_center(:)
    real(ESMF_KIND_R8), allocatable :: xgrd_corner(:,:)
    real(ESMF_KIND_R8), allocatable :: ygrd_corner(:,:)
    logical, allocatable :: land_sea(:)
    integer, allocatable :: grid_dims(:)
    integer :: grid_size, grid_corners, grid_rank
    type(ESMF_Field) :: tmpField
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    if ( noActiveImpFields ) return
    call ESMF_GridCompGet(gcomp, name=cname, vm=vm, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_VMGet(vm, localPet=lpet, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered createimpgrid', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Setup
    !
    ! 1.a Set grid pointers
    !
    impGridID = minval(inpmap)
.eq.    if ( impGridID-999 ) impGridID = 1
    call w3setg ( impGridID, mdse, mdst )
    call w3seti ( impGridID, mdse, mdst )
    call w3seto ( impGridID, mdse, mdst )
.gt.    if ( impGridID0 ) then
      call wmsetm ( impGridID, mdse, mdst )
      nproc = naproc
    else
      nproc = nmproc
    endif
    !
    ! 1.b Compute a 2D subdomain layout based on nproc
    !
    call CalcDecomp( nx, ny, nproc, impHaloWidth, .true., nxproc, nyproc, rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
import fields    ! 1.c Set arraySpec, staggerLoc, and indexFlag for 
    !
    call ESMF_ArraySpecSet( impArraySpec2D, rank=2, &
         typekind=ESMF_TYPEKIND_RX, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    impStaggerLoc = ESMF_STAGGERLOC_CENTER
    impIndexFlag = ESMF_INDEX_GLOBAL
    !
    ! -------------------------------------------------------------------- /
import with 2d subdomain layout    ! 2.  Create ESMF grid for 
    !     Note that the ESMF grid layout is dim1=X, dim2=Y
    !
import grid    ! 2.a Create ESMF 
    !
    select case (iclose)
    case (iclose_none)
      impGrid = ESMF_GridCreateNoPeriDim( &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/1,2/), &
           coordDep2=(/1,2/), &
           regDecomp=(/nxproc,nyproc/), &
           decompFlag=(/ESMF_DECOMP_BALANCED,ESMF_DECOMP_BALANCED/), &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           indexFlag=impIndexFlag, &
           name=trim(cname)//"_import_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case (iclose_smpl)
      impGrid = ESMF_GridCreate1PeriDim( &
           periodicDim=1, &
           poleDim=2, &
           poleKindFlag=(/ESMF_POLEKIND_NONE,ESMF_POLEKIND_NONE/), &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/1,2/), &
           coordDep2=(/1,2/), &
           regDecomp=(/nxproc,nyproc/), &
           decompFlag=(/ESMF_DECOMP_BALANCED,ESMF_DECOMP_BALANCED/), &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           indexFlag=impIndexFlag, &
           name=trim(cname)//"_import_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case default
      write(msg,'(a,i1,a)') 'index closure ',iclose, &
           ' not supported for import grid'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
      rc = ESMF_FAILURE
      return
    endselect
    !
import grid    ! 2.b Add coordinate arrays and land/sea mask to 
    !
    call ESMF_GridAddCoord( impGrid, staggerLoc=impStaggerLoc, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_GridAddItem( impGrid, ESMF_GRIDITEM_MASK, &
         staggerLoc=impStaggerLoc, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
import grid storage    ! 2.c Set flag to indicate that this processor has local 
    !
    call ESMF_GridGet( impGrid, localDECount=ldecnt, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    impGridIsLocal = ldecnt0
    !
import grid    ! 2.d Get exclusive bounds (global index) for 
    !
    if ( impGridIsLocal ) then
      call ESMF_GridGet( impGrid, impStaggerLoc, lde, &
           exclusiveLBound=elb, exclusiveUBound=eub, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
import fields    ! 2.e Set halo widths for 
    !
    if ( impGridIsLocal ) then
      impHaloLWidth = (/impHaloWidth,impHaloWidth/)
      impHaloUWidth = (/impHaloWidth,impHaloWidth/)
      select case (iclose)
      case (iclose_none)
.eq.        if ( elb(1)1  ) impHaloLWidth(1) = 0
.eq.        if ( elb(2)1  ) impHaloLWidth(2) = 0
.eq.        if ( eub(1)nx ) impHaloUWidth(1) = 0
.eq.        if ( eub(2)ny ) impHaloUWidth(2) = 0
      case (iclose_smpl)
.eq.        if ( elb(2)1  ) impHaloLWidth(2) = 0
.eq.        if ( eub(2)ny ) impHaloUWidth(2) = 0
      endselect
    else
      impHaloLWidth = (/0,0/)
      impHaloUWidth = (/0,0/)
    endif
    !
import grid coordinates    ! 2.f Set ESMF 
    !
    if ( impGridIsLocal ) then
      call ESMF_GridGetCoord( impGrid, 1, localDE=lde, &
           staggerLoc=impStaggerLoc, farrayPtr=rptrx, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call ESMF_GridGetCoord( impGrid, 2, localDE=lde, &
           staggerLoc=impStaggerLoc, farrayPtr=rptry, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptrx(ix,iy) = xgrd(iy,ix)
          rptry(ix,iy) = ygrd(iy,ix)
        enddo
      enddo
      nullify(rptrx)
      nullify(rptry)
    endif
    !
import grid land/sea mask values.    ! 2.g Set ESMF 
    !     Land/sea mask is fixed in time and based on excluded points only.
    !
    if ( impGridIsLocal ) then
      call ESMF_GridGetItem( impGrid, ESMF_GRIDITEM_MASK, localDE=lde, &
           staggerLoc=impStaggerLoc, farrayPtr=iptr, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
.ne.          if ( mapsta(iy,ix)0 ) then
            iptr(ix,iy) = maskValueWater
          else
            iptr(ix,iy) = maskValueLand
          endif
        enddo
      enddo
    endif
    !
import grid corner coordinates    ! 2.h Set ESMF 
    !
#ifdef W3_SCRIP
    call ESMF_GridAddCoord( impGrid, &
         staggerLoc=ESMF_STAGGERLOC_CORNER, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
    ! Calculate grid coordinates with help of SCRIP module
    ! It does not return coordinates of top-most row and
    ! right-most column but ESMF expects it. So, top-most row
    ! and right-most column are theated specially in below
    call get_scrip_info_structured(impGridID, &
         xgrd_center, ygrd_center, xgrd_corner, ygrd_corner, &
         land_sea, grid_dims, grid_size, grid_corners, grid_rank)
 
    ! Add corner coordinates
    if ( impGridIsLocal ) then
      ! Retrieve pointers
      call ESMF_GridGetCoord( impGrid, 1, localDE=lde, &
           staggerLoc=ESMF_STAGGERLOC_CORNER, farrayPtr=rptrx, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call ESMF_GridGetCoord( impGrid, 2, localDE=lde, &
           staggerLoc=ESMF_STAGGERLOC_CORNER, farrayPtr=rptry, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
import grid      ! Get exclusive bounds (global index) for 
      ! corner coordinates
      call ESMF_GridGet( impGrid, ESMF_STAGGERLOC_CORNER, lde, &
           exclusiveLBound=elbc, exclusiveUBound=eubc, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
      ! Adjust upper bounds for specific PEs
      ubx = 0
      uby = 0
      if (eubc(1) == grid_dims(1)+1) ubx = -1
      if (eubc(2) == grid_dims(2)+1) uby = -1
 
      ! Fill coordinates
      do iy = elbc(2),eubc(2)+uby
        do ix = elbc(1),eubc(1)+ubx
          irec = (iy-1)*grid_dims(1)+ix
          rptrx(ix,iy) = real(xgrd_corner(1,irec), kind=ESMF_KIND_RX)
          rptry(ix,iy) = real(ygrd_corner(1,irec), kind=ESMF_KIND_RX)
        enddo
      enddo
 
      ! Fill data on top-most row
      if (eubc(2) == grid_dims(2)+1) then
        do ix = elbc(1),eubc(1)+ubx
          rptrx(ix,grid_dims(2)+1) = rptrx(ix,grid_dims(2))+ &
               (rptrx(ix,grid_dims(2))-rptrx(ix,grid_dims(2)-1))
          rptry(ix,grid_dims(2)+1) = rptry(ix,grid_dims(2))+ &
               (rptry(ix,grid_dims(2))-rptry(ix,grid_dims(2)-1))
        end do
      end if
 
      ! Fill data on right-most column
      if (eubc(1) == grid_dims(1)+1) then
        do iy = elbc(2),eubc(2)+uby
          rptrx(grid_dims(1)+1,iy) = rptrx(grid_dims(1),iy)+ &
               (rptrx(grid_dims(1),iy)-rptrx(grid_dims(1)-1,iy))
          rptry(grid_dims(1)+1,iy) = rptry(grid_dims(1),iy)+ &
               (rptry(grid_dims(1),iy)-rptry(grid_dims(1)-1,iy))
        end do
      end if
 
      ! Fill data on top-right corner, single point
.and.      if (eubc(1) == grid_dims(1)+1  eubc(2) == grid_dims(2)+1) then
        rptrx(grid_dims(1)+1,grid_dims(2)+1) = &
             rptrx(grid_dims(1)+1,grid_dims(2))
        rptry(grid_dims(1)+1,grid_dims(2)+1) = &
             rptry(grid_dims(1),grid_dims(2)+1)
      end if
    endif
#endif
    !
    ! -------------------------------------------------------------------- /
import field mask and routehandle halo update    ! 3.  Create 
    !
import grid land/sea mask.    ! 3.a Create field for 
    !
    impMask = ESMF_FieldCreate( impGrid, impArraySpec2D, &
         totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
         staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
         name='mask', rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
import field halo routehandle    ! 3.b Store 
    !
    call ESMF_FieldHaloStore( impMask, routeHandle=impHaloRH, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
import field land/sea mask values and update halos    ! 3.c Set 
    !
    if ( impGridIsLocal ) then
      call ESMF_FieldGet( impMask, localDE=lde, farrayPtr=rptr, &
           totalLBound=tlb, totalUBound=tub, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = iptr(ix,iy)
        enddo
      enddo
    endif
 
    call ESMF_FieldHalo( impMask, impHaloRH, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_CREATEIMPGRID)
    write(msg,'(a,i6)')  '     impgridid: ',impGridID
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,l6)')  'impgridislocal: ',impGridIsLocal
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '        nx, ny: ',nx,ny
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') 'naproc, nmproc: ',naproc,nmproc
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,i6)')  '         nproc: ',nproc
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') 'nxproc, nyproc: ',nxproc,nyproc
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           elb: ',elb(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           eub: ',eub(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           tlb: ',tlb(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           tub: ',tub(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') ' imphalolwidth: ',impHaloLWidth(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') ' imphalouwidth: ',impHaloUWidth(:)
    call ESMF_LogWrite(trim(cname)//': createimpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    call ESMF_FieldWrite( impMask, &
         "wmesmfmd_createimpgrid_import_mask.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    tmpField = ESMF_FieldCreate( impGrid, impArraySpec2D, &
         totalLWidth=impHaloLWidth, totalUWidth=impHaloUWidth, &
         staggerLoc=impStaggerLoc, indexFlag=impIndexFlag, &
         name='temp', rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impGridIsLocal ) then
      call ESMF_FieldGet( tmpField, localDE=lde, farrayPtr=rptr, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    if ( impGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = xgrd(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createimpgrid_import_xgrd.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = ygrd(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createimpgrid_import_ygrd.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = hpfac(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createimpgrid_import_hpfac.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = hqfac(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createimpgrid_import_hqfac.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( impGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = lpet
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createimpgrid_import_dcomp.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_FieldDestroy( tmpField, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
#endif
    !
    rc = ESMF_SUCCESS
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving createimpgrid', ESMF_LOGMSG_INFO)
    !/
    !/ End of CreateImpGrid ---------------------------------------------- /
    !/
  end subroutine CreateImpGrid
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CreateExpGrid"
  !>
  !> @brief Create ESMF grid for export fields
  !>
  !> @param      gcomp Gridded component
  !> @param[out] rc    Return code
  !>
  !> @author T. J. Campbell  @date 20-Jan-2017
  !>
  subroutine CreateExpGrid ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !       !/SHRD  Switch for shared / distributed memory architecture.
    !       !/DIST  Id.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    integer :: nproc, nxproc, nyproc
    integer, parameter :: lde = 0
    integer :: ldecnt
    integer :: ix, iy, isea, jsea, irec, k, ubx, uby
    integer :: elb(2), eub(2), elbc(2), eubc(2)
    integer :: tlb(2), tub(2)
    integer(ESMF_KIND_I4), pointer :: iptr(:,:)
    real(ESMF_KIND_RX), pointer :: rptrx(:,:), rptry(:,:)
    real(ESMF_KIND_RX), pointer :: rptr(:,:)
    real(ESMF_KIND_R8), allocatable :: xgrd_center(:)
    real(ESMF_KIND_R8), allocatable :: ygrd_center(:)
    real(ESMF_KIND_R8), allocatable :: xgrd_corner(:,:)
    real(ESMF_KIND_R8), allocatable :: ygrd_corner(:,:)
    logical, allocatable :: land_sea(:)
    integer, allocatable :: grid_dims(:)
    integer :: grid_size, grid_corners, grid_rank
    integer :: arbIndexCount
    integer, allocatable :: arbIndexList(:,:)
    type(ESMF_Field) :: nField, eField
    type(ESMF_Field) :: tmpField
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    if ( noActiveExpFields ) return
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered createexpgrid', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Setup
    !
    ! 1.a Set grid pointers
    !
    !TODO: only export from one grid
    !expGridID set from config input (default = 1)
    call w3setg ( expGridID, mdse, mdst )
    call w3setw ( expGridID, mdse, mdst )
    call w3seta ( expGridID, mdse, mdst )
    call w3seti ( expGridID, mdse, mdst )
    call w3seto ( expGridID, mdse, mdst )
    call wmsetm ( expGridID, mdse, mdst )
    natGridID = expGridID
    nproc = naproc
    !
    ! 1.b Compute a 2D subdomain layout based on nproc
    !
    call CalcDecomp( nx, ny, nproc, expHaloWidth, .true., nxproc, nyproc, rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 1.c Set arraySpec, staggerLoc, and indexFlag for export fields
    !
    call ESMF_ArraySpecSet( expArraySpec2D, rank=2, &
         typekind=ESMF_TYPEKIND_RX, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_ArraySpecSet( expArraySpec3D, rank=3, &
         typekind=ESMF_TYPEKIND_RX, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    expStaggerLoc = ESMF_STAGGERLOC_CENTER
    expIndexFlag = ESMF_INDEX_GLOBAL
    !
    ! 1.d Set arraySpec, staggerLoc, and indexFlag for native fields
    !
    call ESMF_ArraySpecSet( natArraySpec2D, rank=1, &
         typekind=ESMF_TYPEKIND_RX, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_ArraySpecSet( natArraySpec3D, rank=2, &
         typekind=ESMF_TYPEKIND_RX, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    natStaggerLoc = ESMF_STAGGERLOC_CENTER
    natIndexFlag = ESMF_INDEX_DELOCAL
    !
    ! 1.e Get z-levels for 3D export fields
    !
    call GetZlevels( rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    write(msg,'(a)') 'input z-level file: '//trim(zlfile)
    call ESMF_LogWrite(trim(cname)//': '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a)') 'table of z-levels'
    call ESMF_LogWrite(trim(cname)//': '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a)') '   index         z'
    call ESMF_LogWrite(trim(cname)//': '//trim(msg), ESMF_LOGMSG_INFO)
    do k=1,nz
      write(msg,'(i8,1f10.2)') k, zl(k)
      call ESMF_LogWrite(trim(cname)//': '//trim(msg), ESMF_LOGMSG_INFO)
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Create ESMF grid for export with 2D subdomain layout
    !     Note that the ESMF grid layout is dim1=X, dim2=Y
    !
    ! 2.a Create ESMF export grid
    !
    select case (iclose)
    case (iclose_none)
      expGrid = ESMF_GridCreateNoPeriDim( &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/1,2/), &
           coordDep2=(/1,2/), &
           regDecomp=(/nxproc,nyproc/), &
           decompFlag=(/ESMF_DECOMP_BALANCED,ESMF_DECOMP_BALANCED/), &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           indexFlag=expIndexFlag, &
           name=trim(cname)//"_export_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case (iclose_smpl)
      expGrid = ESMF_GridCreate1PeriDim( &
           periodicDim=1, &
           poleDim=2, &
           poleKindFlag=(/ESMF_POLEKIND_NONE,ESMF_POLEKIND_NONE/), &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/1,2/), &
           coordDep2=(/1,2/), &
           regDecomp=(/nxproc,nyproc/), &
           decompFlag=(/ESMF_DECOMP_BALANCED,ESMF_DECOMP_BALANCED/), &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           indexFlag=expIndexFlag, &
           name=trim(cname)//"_export_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case default
      write(msg,'(a,i1,a)') 'index closure ',iclose, &
           ' not supported for export grid'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
      rc = ESMF_FAILURE
      return
    endselect
    !
    ! 2.b Add coordinate arrays and land/sea mask to export grid
    !
    call ESMF_GridAddCoord( expGrid, staggerLoc=expStaggerLoc, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_GridAddItem( expGrid, ESMF_GRIDITEM_MASK, &
         staggerLoc=expStaggerLoc, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 2.c Set flag to indicate that this processor has local export grid storage
    !
    call ESMF_GridGet( expGrid, localDECount=ldecnt, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    expGridIsLocal = ldecnt0
    !
    ! 2.d Get exclusive bounds (global index) for export grid
    !
    if ( expGridIsLocal ) then
      call ESMF_GridGet( expGrid, expStaggerLoc, lde, &
           exclusiveLBound=elb, exclusiveUBound=eub, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    !
    ! 2.e Set halo widths for export fields
    !
    if ( expGridIsLocal ) then
      expHaloLWidth = (/expHaloWidth,expHaloWidth/)
      expHaloUWidth = (/expHaloWidth,expHaloWidth/)
      select case (iclose)
      case (iclose_none)
.eq.        if ( elb(1)1  ) expHaloLWidth(1) = 0
.eq.        if ( elb(2)1  ) expHaloLWidth(2) = 0
.eq.        if ( eub(1)nx ) expHaloUWidth(1) = 0
.eq.        if ( eub(2)ny ) expHaloUWidth(2) = 0
      case (iclose_smpl)
.eq.        if ( elb(2)1  ) expHaloLWidth(2) = 0
.eq.        if ( eub(2)ny ) expHaloUWidth(2) = 0
      endselect
    else
      expHaloLWidth = (/0,0/)
      expHaloUWidth = (/0,0/)
    endif
    !
    ! 2.f Set ESMF export grid coordinate
    !
    if ( expGridIsLocal ) then
      call ESMF_GridGetCoord( expGrid, 1, localDE=lde, &
           staggerLoc=expStaggerLoc, farrayPtr=rptrx, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call ESMF_GridGetCoord( expGrid, 2, localDE=lde, &
           staggerLoc=expStaggerLoc, farrayPtr=rptry, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptrx(ix,iy) = xgrd(iy,ix)
          rptry(ix,iy) = ygrd(iy,ix)
        enddo
      enddo
    endif
    !
    ! 2.g Set ESMF export grid land/sea mask values.
    !     Land/sea mask is fixed in time and based on excluded points only.
    !
    if ( expGridIsLocal ) then
      call ESMF_GridGetItem( expGrid, ESMF_GRIDITEM_MASK, localDE=lde, &
           staggerLoc=expStaggerLoc, farrayPtr=iptr, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
.ne.          if ( mapsta(iy,ix)0 ) then
            iptr(ix,iy) = maskValueWater
          else
            iptr(ix,iy) = maskValueLand
          endif
        enddo
      enddo
    endif
    !
    ! 2.h Set ESMF export grid corner coordinates
    !
#ifdef W3_SCRIP
    call ESMF_GridAddCoord( expGrid, &
         staggerLoc=ESMF_STAGGERLOC_CORNER, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
    ! Calculate grid coordinates with help of SCRIP module
    ! It does not return coordinates of top-most row and
    ! right-most column but ESMF expects it. So, top-most row
    ! and right-most column are theated specially in below
    call get_scrip_info_structured(expGridID, &
         xgrd_center, ygrd_center, xgrd_corner, ygrd_corner, &
         land_sea, grid_dims, grid_size, grid_corners, grid_rank)
 
    ! Add corner coordinates
    if ( impGridIsLocal ) then
      ! Retrieve pointers
      call ESMF_GridGetCoord( expGrid, 1, localDE=lde, &
           staggerLoc=ESMF_STAGGERLOC_CORNER, farrayPtr=rptrx, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      call ESMF_GridGetCoord( expGrid, 2, localDE=lde, &
           staggerLoc=ESMF_STAGGERLOC_CORNER, farrayPtr=rptry, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
      ! Get exclusive bounds (global index) for export grid
      ! corner coordinates
      call ESMF_GridGet( impGrid, ESMF_STAGGERLOC_CORNER, lde, &
           exclusiveLBound=elbc, exclusiveUBound=eubc, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
 
      ! Adjust upper bounds for specific PEs
      ubx = 0
      uby = 0
      if (eubc(1) == grid_dims(1)+1) ubx = -1
      if (eubc(2) == grid_dims(2)+1) uby = -1
 
      ! Fill coordinates
      do iy = elbc(2),eubc(2)+uby
        do ix = elbc(1),eubc(1)+ubx
          irec = (iy-1)*grid_dims(1)+ix
          rptrx(ix,iy) = real(xgrd_corner(1,irec), kind=ESMF_KIND_RX)
          rptry(ix,iy) = real(ygrd_corner(1,irec), kind=ESMF_KIND_RX)
        enddo
      enddo
 
      ! Fill data on top-most row
      if (eubc(2) == grid_dims(2)+1) then
        do ix = elbc(1),eubc(1)+ubx
          rptrx(ix,grid_dims(2)+1) = rptrx(ix,grid_dims(2))+ &
               (rptrx(ix,grid_dims(2))-rptrx(ix,grid_dims(2)-1))
          rptry(ix,grid_dims(2)+1) = rptry(ix,grid_dims(2))+ &
               (rptry(ix,grid_dims(2))-rptry(ix,grid_dims(2)-1))
        end do
      end if
 
      ! Fill data on right-most column
      if (eubc(1) == grid_dims(1)+1) then
        do iy = elbc(2),eubc(2)+uby
          rptrx(grid_dims(1)+1,iy) = rptrx(grid_dims(1),iy)+ &
               (rptrx(grid_dims(1),iy)-rptrx(grid_dims(1)-1,iy))
          rptry(grid_dims(1)+1,iy) = rptry(grid_dims(1),iy)+ &
               (rptry(grid_dims(1),iy)-rptry(grid_dims(1)-1,iy))
        end do
      end if
 
      ! Fill data on top-right corner, single point
.and.      if (eubc(1) == grid_dims(1)+1  eubc(2) == grid_dims(2)+1) then
        rptrx(grid_dims(1)+1,grid_dims(2)+1) = &
             rptrx(grid_dims(1)+1,grid_dims(2))
        rptry(grid_dims(1)+1,grid_dims(2)+1) = &
             rptry(grid_dims(1),grid_dims(2)+1)
      end if
    end if
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Create export field mask and routeHandle halo update
    !
    ! 3.a Create field for export grid land/sea mask.
    !
    expMask = ESMF_FieldCreate( expGrid, expArraySpec2D, &
         totalLWidth=expHaloLWidth, totalUWidth=expHaloUWidth, &
         staggerLoc=expStaggerLoc, indexFlag=expIndexFlag, &
         name='mask', rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 3.b Store export field halo routeHandle
    !
    call ESMF_FieldHaloStore( expMask, routeHandle=expHaloRH, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 3.c Set export field land/sea mask values and update halos
    !
    if ( expGridIsLocal ) then
      call ESMF_FieldGet( expMask, localDE=lde, farrayPtr=rptr, &
           totalLBound=tlb, totalUBound=tub, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = iptr(ix,iy)
        enddo
      enddo
    endif
 
    call ESMF_FieldHalo( expMask, expHaloRH, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! 4.  Create ESMF grid with arbitrary domain decomposition to match
    !     the native domain decomposition of the non-excluded points
    !     Note that the native grid layout is dim1=X, dim2=Y
    !     Note that coordinates and mask are not needed since this
    !     grid is only used to define fields for a redist operation
    !
    ! 4.a Set flag to indicate that this processor has local native grid storage
    !
.gt..and..le.    natGridIsLocal = iaproc  0  iaproc  naproc
    !
    ! 4.b Setup arbitrary sequence index list
    !
    do ipass = 1,2
.eq.      if (ipass2) then
        allocate (arbIndexList(arbIndexCount,2), stat=rc)
        if (ESMF_LogFoundAllocError(rc, PASSTHRU)) return
      endif
      arbIndexCount = 0
      ! list local native grid non-excluded points
      if ( natGridIsLocal ) then
        do jsea = 1,nseal
#ifdef W3_DIST
          isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
          isea = jsea
#endif
          arbIndexCount = arbIndexCount+1
.eq.          if (ipass2) then
            ix = mapsf(isea,1)
            iy = mapsf(isea,2)
            ! native grid layout: dim1=X, dim2=Y
            arbIndexList(arbIndexCount,1) = ix
            arbIndexList(arbIndexCount,2) = iy
          endif
        enddo
      endif
      ! list local export grid excluded points
      if ( expGridIsLocal ) then
        do iy = elb(2),eub(2)
          do ix = elb(1),eub(1)
.ne.            if ( mapsta(iy,ix)0 ) cycle ! skip non-excluded point
            arbIndexCount = arbIndexCount+1
.eq.            if (ipass2) then
              ! native grid layout: dim1=X, dim2=Y
              arbIndexList(arbIndexCount,1) = ix
              arbIndexList(arbIndexCount,2) = iy
            endif
          enddo
        enddo
      endif
    enddo !ipass
    !
    ! 4.c Create ESMF native grid
    !
    select case (iclose)
    case (iclose_none)
      natGrid = ESMF_GridCreateNoPeriDim( &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/ESMF_DIM_ARB,ESMF_DIM_ARB/), &
           coordDep2=(/ESMF_DIM_ARB,ESMF_DIM_ARB/), &
           arbIndexCount=arbIndexCount, &
           arbIndexList=arbIndexList, &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           name=trim(cname)//"_native_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case (iclose_smpl)
      natGrid = ESMF_GridCreate1PeriDim( &
           periodicDim=1, &
           poleDim=2, &
           poleKindFlag=(/ESMF_POLEKIND_NONE,ESMF_POLEKIND_NONE/), &
           minIndex=(/ 1, 1/), &
           maxIndex=(/nx,ny/), &
           coordDep1=(/ESMF_DIM_ARB,ESMF_DIM_ARB/), &
           coordDep2=(/ESMF_DIM_ARB,ESMF_DIM_ARB/), &
           arbIndexCount=arbIndexCount, &
           arbIndexList=arbIndexList, &
           coordTypeKind=ESMF_TYPEKIND_RX, &
           coordSys=ESMF_COORDSYS_SPH_DEG, &
           name=trim(cname)//"_native_grid", rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    case default
      write(msg,'(a,i1,a)') 'index closure ',iclose, &
           ' not supported for native grid'
      call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR)
      rc = ESMF_FAILURE
      return
    endselect
    !
    ! 4.d Deallocate arbitrary sequence index list
    !
    deallocate (arbIndexList, stat=rc)
    if (ESMF_LogFoundDeallocError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! 5.  Create route handle for redist between native grid domain
    !     decomposition and the export grid domain decomposition
    !
    ! 5.a Create temporary fields
    !
    nField = ESMF_FieldCreate( natGrid, natArraySpec2D, &
         staggerLoc=natStaggerLoc, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    eField = ESMF_FieldCreate( expGrid, expArraySpec2D, &
         totalLWidth=expHaloLWidth, totalUWidth=expHaloUWidth, &
         staggerLoc=expStaggerLoc, indexFlag=expIndexFlag, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 5.b Store route handle
    !
    call ESMF_FieldRedistStore( nField, eField, n2eRH, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! 5.c Clean up
    !
    call ESMF_FieldDestroy( nField, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_FieldDestroy( eField, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !
    ! -------------------------------------------------------------------- /
    ! Post
    !
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_CREATEEXPGRID)
    write(msg,'(a,i6)')  '     expgridid: ',expGridID
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,l6)')  'expgridislocal: ',expGridIsLocal
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '        nx, ny: ',nx,ny
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') 'naproc, nmproc: ',naproc,nmproc
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,i6)')  '         nproc: ',nproc
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') 'nxproc, nyproc: ',nxproc,nyproc
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           elb: ',elb(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           eub: ',eub(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           tlb: ',tlb(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') '           tub: ',tub(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') ' exphalolwidth: ',expHaloLWidth(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    write(msg,'(a,2i6)') ' exphalouwidth: ',expHaloUWidth(:)
    call ESMF_LogWrite(trim(cname)//': createexpgrid: '//trim(msg), ESMF_LOGMSG_INFO)
    call ESMF_FieldWrite( expMask, &
         "wmesmfmd_createexpgrid_export_mask.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    tmpField = ESMF_FieldCreate( expGrid, expArraySpec2D, &
         totalLWidth=expHaloLWidth, totalUWidth=expHaloUWidth, &
         staggerLoc=expStaggerLoc, indexFlag=expIndexFlag, &
         name='temp', rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expGridIsLocal ) then
      call ESMF_FieldGet( tmpField, localDE=lde, farrayPtr=rptr, rc=rc )
      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    endif
    if ( expGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = xgrd(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createexpgrid_export_xgrd.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = ygrd(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createexpgrid_export_ygrd.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = hpfac(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createexpgrid_export_hpfac.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = hqfac(iy,ix)
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createexpgrid_export_hqfac.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    if ( expGridIsLocal ) then
      do iy = elb(2),eub(2)
        do ix = elb(1),eub(1)
          rptr(ix,iy) = lpet
        enddo
      enddo
    endif
    call ESMF_FieldWrite( tmpField, &
         "wmesmfmd_createexpgrid_export_dcomp.nc", overwrite=.true., rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
    call ESMF_FieldDestroy( tmpField, rc=rc )
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
#endif
    !
    rc = ESMF_SUCCESS
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': leaving createexpgrid', ESMF_LOGMSG_INFO)
    !/
    !/ End of CreateExpGrid ---------------------------------------------- /
    !/
  end subroutine CreateExpGrid
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CreateImpMesh"
  !>
import fields.  !> @brief Create ESMF mesh (unstructured) for 
  !>
import using the unstructured mesh  !> @details Create an ESMF Mesh for 
import mesh is not  !>  description in W3GDATMD. At present, this 
  !>  domain decomposed, but instead is defined on PET 0 only. (In
  !>  future, when the unstructured mesh will run on domain decomposition,
  !>  we will use that decomposition.)
  !>
  !> @param      gcomp Gridded component
  !> @param[out] rc    Return code
  !>
  !> @author A. J. van der Westhuysen  @date 28-Feb-2018
  !>
  subroutine CreateImpMesh ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |     A. J. van der Westhuysen      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         28-FEB_2018 |
    !/                  +-----------------------------------+
    !/
    !/    28-Feb-2018 : Origination.                        ( version 6.06 )
    !/
    !  1. Purpose :
    !
import fields    !     Create ESMF mesh (unstructured) for 
    !
    !  2. Method :
    !
import using the unstructured mesh description    !     Create an ESMF Mesh for 
import mesh is not domain decomposed,    !     in W3GDATMD. At present, this 
    !     but instead is defined on PET 0 only. (In future, when the unstructured
    !     mesh will run on domain decomposition, we will use that decomposition.)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_PDLIB
    use yowNodepool, only: npa, iplg, nodes_global
    use yowElementpool, only: ne, ielg, INE
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(ESMF_GridComp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    character(128) :: msg
    integer :: nproc, nxproc, nyproc, n, nfac, irp
    real    :: gr, rp, pr, diff
    integer, parameter :: lde = 0
    integer :: ldecnt
    integer :: i, j, pos, ix, iy
    integer(ESMF_KIND_I4), pointer :: iptr(:,:)
    real(ESMF_KIND_RX), pointer :: rptrx(:,:), rptry(:,:)
    real(ESMF_KIND_RX), pointer :: rptr(:,:)
    type(ESMF_Field) :: tmpField
    integer(ESMF_KIND_I4), allocatable :: nodeIds(:)
    real(ESMF_KIND_R8), allocatable    :: nodeCoords(:)
    integer(ESMF_KIND_I4), allocatable :: nodeOwners(:)
    integer(ESMF_KIND_I4), allocatable :: elemIds(:)
    integer(ESMF_KIND_I4), allocatable :: elemTypes(:)
    integer(ESMF_KIND_I4), allocatable :: elemConn(:)
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = ESMF_SUCCESS
    if ( noActiveImpFields ) return
    call ESMF_GridCompGet(gcomp, name=cname, rc=rc)
    if (ESMF_LogFoundError(rc, PASSTHRU)) return
.gt.    if (verbosity0) call ESMF_LogWrite(trim(cname)// &
         ': entered createimpmesh', ESMF_LOGMSG_INFO)
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Setup
    !
    ! 1.a Set grid pointers
    !
    impGridID = minval(inpmap)
.eq.    if ( impGridID-999 ) impGridID = 1
    call w3setg ( impGridID, mdse, mdst )
    call w3seti ( impGridID, mdse, mdst )
    call w3seto ( impGridID, mdse, mdst )
.gt.    if ( impGridID0 ) then
      call wmsetm ( impGridID, mdse, mdst )
      nproc = naproc
    else
      nproc = nmproc
    endif
    !
import fields    ! 1.b Set arraySpec, staggerLoc, and indexFlag for 
    !
    !      call ESMF_ArraySpecSet( impArraySpec2D, rank=2, &
    !        typekind=ESMF_TYPEKIND_RX, rc=rc )
    !      if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !      impStaggerLoc = ESMF_STAGGERLOC_CENTER
    !      impIndexFlag = ESMF_INDEX_GLOBAL
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Create ESMF mesh for import, currently without domain decomposition
    !     Note that the ESMF grid layout is dim1=X, dim2=Y
    !
import mesh    ! 2.a Create ESMF 
    !
    ! Allocate and fill the node id array.
#ifdef W3_PDLIB
.EQV.    if ( LPDLIB  .FALSE. ) then
#endif
      allocate(nodeIds(NX))
      do i = 1,NX
        nodeIds(i)=i
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: The global id's of the nodes resident on this processor      !        ESMF Definition: The global id'
      !        -------------------------------------------------------------------
      !        Allocate global node ids, including ghost nodes (npa=np+ng)
      allocate(nodeIds(npa))
      do i = 1,npa
        nodeIds(i)=iplg(i)
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodeids=', &
    !          ESMF_LOGMSG_INFO)
    !         do i = 1,npa
    !         write(msg,*) trim(cname)//': nodeids(i)',i, &
    !          ' ',nodeIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodeids=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NX
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill node coordinate array.
    ! Since this is a 2D Mesh the size is 2x the
    ! number of nodes.
#ifdef W3_PDLIB
.EQV.    if ( LPDLIB  .FALSE. ) then
#endif
      allocate(nodeCoords(2*NX))
      do i = 1,NX
        do j = 1,2
          pos=2*(i-1)+j
          if (j == 1) then
            nodeCoords(pos) = xgrd(1,i)
          else
            nodeCoords(pos) = ygrd(1,i)
          endif
        enddo
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Physical coordinates of the nodes
      !        -------------------------------------------------------------------
      allocate(nodeCoords(2*npa))
      do i = 1,npa
        do j = 1,2
          pos=2*(i-1)+j
          if ( j == 1) then
            nodeCoords(pos) = xgrd(1,iplg(i))
          else
            nodeCoords(pos) = ygrd(1,iplg(i))
          endif
        enddo
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodecoords=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(2*npa)
    !        write(msg,*) trim(cname)//': nodecoords(i)',i, &
    !          ' ',nodeCoords(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodecoords=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(2*NX)
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeCoords(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the node owner array.
    ! Since this mesh is all on PET 0, it’s just set to all 0.
#ifdef W3_PDLIB
.EQV.    if ( LPDLIB  .FALSE. ) then
#endif
      allocate(nodeOwners(NX))
      nodeOwners=0 ! everything on PET 0
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Processor that owns the node
      !        -------------------------------------------------------------------
      allocate(nodeOwners(npa))
      nodeOwners=nodes_global(iplg(1:npa))%domainID-1
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodeowners=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,npa
    !        write(msg,*) trim(cname)//': nodeowners(i)',i, &
    !          ' ',nodeOwners(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': in createimpmesh, nodeowners=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NX
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeOwners(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element id array.
#ifdef W3_PDLIB
.EQV.    if ( LPDLIB  .FALSE. ) then
#endif
      allocate(elemIds(NTRI))
      do i = 1,NTRI
        elemIds(i)=i
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: The global id's of the elements resident on this processor
      !        -------------------------------------------------------------------
      allocate(elemids(ne))
      do i = 1,ne
        elemids(i)=ielg(i)
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpMesh, elemIds=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,ne
    !        write(msg,*) trim(cname)//': elemIds(i)',i, &
    !          ' ',elemIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpMesh, elemIds=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NTRI
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element topology type array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(elemtypes(ntri))
      do i = 1,ntri
        elemtypes(i)=esmf_meshelemtype_tri
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Topology of the given element (one of ESMF_MeshElement)
      !        -------------------------------------------------------------------
      allocate(elemtypes(ne))
      do i = 1,ne
        elemtypes(i)=esmf_meshelemtype_tri
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpM, elemTypes=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,ne
    !        write(msg,*) trim(cname)//': elemTypes(i)',i, &
    !          ' ',elemTypes(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpM, elemTypes=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NTRI
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemTypes(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element connection type array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(elemconn(3*ntri))
      do i = 1,ntri
        do j = 1,3
          pos=3*(i-1)+j
          elemconn(pos)=trigp(j,i)
        enddo
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Connectivity table. The number of entries should
      !        be equal to the number of nodes in the given topology. The indices
      !        should be the local index (1 based) into the array of nodes that
      !        was declared with MeshAddNodes.
      !        -------------------------------------------------------------------
      !        > INE is local element array. it stores the local node IDs
      !        > first index from 1 to 3.
      !        > second index from 1 to ne.
      allocate(elemconn(3*ne))
      do i = 1,ne
        do j = 1,3
          pos=3*(i-1)+j
          elemconn(pos)=ine(j,i)
        enddo
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpMesh, elemConn=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(3*ne)
    !        write(msg,*) trim(cname)//': elemConn(i)',i, &
    !          ' ',elemConn(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateImpMesh, elemConn=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(3*NTRI)
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemConn(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    impmesh = esmf_meshcreate( parametricdim=2,spatialdim=2, &
         nodeids=nodeids, nodecoords=nodecoords, &
         nodeowners=nodeowners, elementids=elemids,&
         elementtypes=elemtypes, elementconn=elemconn, &
         rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    deallocate(nodeids)
    deallocate(nodecoords)
    deallocate(nodeowners)
    deallocate(elemids)
    deallocate(elemtypes)
    deallocate(elemconn)
 
    call esmf_logwrite(trim(cname)//': In CreateImpMesh, created impMesh', &
         esmf_logmsg_info)
    !
    rc = esmf_success
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': leaving CreateImpMesh', esmf_logmsg_info)
    !/
    !/ End of CreateImpMesh ---------------------------------------------- /
    !/
  end subroutine createimpmesh
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CreateExpMesh"
 
  subroutine createexpmesh ( gcomp, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |     A. J. van der Westhuysen      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         28-FEB-2018 |
    !/                  +-----------------------------------+
    !/
    !/    28-Feb-2018 : Origination.                        ( version 6.06 )
    !/
    !  1. Purpose :
    !
    !     Create ESMF mesh (unstructured) for export fields
    !
    !  2. Method :
    !
    !     Create an ESMF Mesh for export using the unstructured mesh description
    !     in W3GDATMD. At present, this export mesh is not domain decomposed,
    !     but instead is defined on PET 0 only. (In future, when the unstructured
    !     mesh will run on domain decomposition, we will use that decomposition.)
    !
    !     Since the internal parallel data is currently stored accross grid points
    !     in a "card deck" fashion, we will define an intermediate native grid, as
    !     is done for regular/curvilinear grids, and perform an ESMF regrid to the
    !     export mesh. This code segment is taken from T. J. Campbell, and
    !     modified to 1D, because the internal data structure for unstructred
    !     meshes is an array with dimensions [NX,NY=1].
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       gcomp     Type   I/O Gridded component
    !       rc        Int.   O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !       !/SHRD  Switch for shared / distributed memory architecture.
    !       !/DIST  Id.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_PDLIB
    use yownodepool, only: npa, iplg, nodes_global
    use yowelementpool, only: ne, ielg, ine
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_gridcomp) :: gcomp
    integer,intent(out) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    character(128) :: msg
    integer :: nproc, nxproc, nyproc, n, nfac, irp
    real    :: gr, rp, pr, diff
    integer, parameter :: lde = 0
    integer :: ldecnt
    integer :: i, j, pos, ix, iy, isea, jsea, iproc
    integer :: elb(2), eub(2)
    integer(ESMF_KIND_I4), pointer :: iptr(:,:)
    real(ESMF_KIND_RX), pointer :: rptrx(:,:), rptry(:,:)
    real(ESMF_KIND_RX), pointer :: rptr(:,:)
    integer :: arbIndexCount
    integer, allocatable :: arbIndexList(:,:)
    type(esmf_field) :: nField, eField
    type(esmf_field) :: tmpField
    integer(ESMF_KIND_I4), allocatable :: nodeIds(:)
    real(ESMF_KIND_R8), allocatable    :: nodeCoords(:)
    integer(ESMF_KIND_I4), allocatable :: nodeOwners(:)
    integer(ESMF_KIND_I4), allocatable :: elemIds(:)
    integer(ESMF_KIND_I4), allocatable :: elemTypes(:)
    integer(ESMF_KIND_I4), allocatable :: elemConn(:)
    !
    ! -------------------------------------------------------------------- /
    ! Prep
    !
    rc = esmf_success
    if ( noactiveexpfields ) return
    call esmf_gridcompget(gcomp, name=cname, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': entered CreateExpMesh', esmf_logmsg_info)
    !
    ! -------------------------------------------------------------------- /
    !     Set flag to indicate that this processor has local native grid storage
    !
    natgridislocal = iaproc .gt. 0 .and. iaproc .le. naproc
 
    ! 1.  Setup
    !
    ! 1.a Set grid pointers
    !
    expgridid = 1 !TODO: only export from grid 1
    call w3setg ( expgridid, mdse, mdst )
    call w3setw ( expgridid, mdse, mdst )
    call w3seta ( expgridid, mdse, mdst )
    call w3seti ( expgridid, mdse, mdst )
    call w3seto ( expgridid, mdse, mdst )
    call wmsetm ( expgridid, mdse, mdst )
    natgridid = expgridid
    nproc = naproc
    !
    ! 1.b Set arraySpec, staggerLoc, and indexFlag for native fields.
    !     NOTE: For unstructured meshes the native grid is a 1D array (NY=1)
    !
    call esmf_arrayspecset( natarrayspec1d, rank=1, &
         typekind=esmf_typekind_rx, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    natstaggerloc = esmf_staggerloc_center
    natindexflag = esmf_index_delocal
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Create ESMF mesh for export
    !
    ! 2.a Create ESMF export mesh
    !
    ! Allocate and fill the node id array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(nodeids(nx))
      do i = 1,nx
        nodeids(i)=i
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: The global id's of the nodes resident on this processor
      !        -------------------------------------------------------------------
      !        Allocate global node ids, including ghost nodes (npa=np+ng)
      allocate(nodeids(npa))
      do i = 1,npa
        nodeids(i)=iplg(i)
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeIds=', &
    !          ESMF_LOGMSG_INFO)
    !         do i = 1,npa
    !         write(msg,*) trim(cname)//': nodeIds(i)',i, &
    !          ' ',nodeIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeIds=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NX
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill node coordinate array.
    ! Since this is a 2D Mesh the size is 2x the
    ! number of nodes.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(nodecoords(2*nx))
      do i = 1,nx
        do j = 1,2
          pos=2*(i-1)+j
          if (j == 1) then
            nodecoords(pos) = xgrd(1,i)
          else
            nodecoords(pos) = ygrd(1,i)
          endif
        enddo
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Physical coordinates of the nodes
      !        -------------------------------------------------------------------
      allocate(nodecoords(2*npa))
      do i = 1,npa
        do j = 1,2
          pos=2*(i-1)+j
          if ( j == 1) then
            nodecoords(pos) = xgrd(1,iplg(i))
          else
            nodecoords(pos) = ygrd(1,iplg(i))
          endif
        enddo
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeCoords=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(2*npa)
    !        write(msg,*) trim(cname)//': nodeCoords(i)',i, &
    !          ' ',nodeCoords(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeCoords=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(2*NX)
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeCoords(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the node owner array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(nodeowners(nx))
      nodeowners=0 ! TODO: For now, export everything via PET 0
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Processor that owns the node
      !        -------------------------------------------------------------------
      allocate(nodeowners(npa))
      nodeowners=nodes_global(iplg(1:npa))%domainID-1
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeOwners=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,npa
    !        write(msg,*) trim(cname)//': nodeOwners(i)',i, &
    !          ' ',nodeOwners(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, nodeOwners=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NX
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',nodeOwners(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element id array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(elemids(ntri))
      do i = 1,ntri
        elemids(i)=i
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: The global id's of the elements resident on this processor
      !        -------------------------------------------------------------------
      allocate(elemids(ne))
      do i = 1,ne
        elemids(i)=ielg(i)
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemIds=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,ne
    !        write(msg,*) trim(cname)//': elemIds(i)',i, &
    !          ' ',elemIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemIds=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NTRI
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemIds(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element topology type array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(elemtypes(ntri))
      do i = 1,ntri
        elemtypes(i)=esmf_meshelemtype_tri
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Topology of the given element (one of ESMF_MeshElement)
      !        -------------------------------------------------------------------
      allocate(elemtypes(ne))
      do i = 1,ne
        elemtypes(i)=esmf_meshelemtype_tri
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemTypes=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,ne
    !        write(msg,*) trim(cname)//': elemTypes(i)',i, &
    !          ' ',elemTypes(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemTypes=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,NTRI
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemTypes(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    ! Allocate and fill the element connection type array.
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      allocate(elemconn(3*ntri))
      do i = 1,ntri
        do j = 1,3
          pos=3*(i-1)+j
          elemconn(pos)=trigp(j,i)
        enddo
      enddo
#ifdef W3_PDLIB
    else
      !        -------------------------------------------------------------------
      !        ESMF Definition: Connectivity table. The number of entries should
      !        be equal to the number of nodes in the given topology. The indices
      !        should be the local index (1 based) into the array of nodes that
      !        was declared with MeshAddNodes.
      !        -------------------------------------------------------------------
      !        > INE is local element array. it stores the local node IDs
      !        > first index from 1 to 3.
      !        > second index from 1 to ne.
      allocate(elemconn(3*ne))
      do i = 1,ne
        do j = 1,3
          pos=3*(i-1)+j
          elemconn(pos)=ine(j,i)
        enddo
      enddo
    endif
    !
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemConn=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(3*ne)
    !        write(msg,*) trim(cname)//': elemConn(i)',i, &
    !          ' ',elemConn(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
#endif
 
    !      call ESMF_LogWrite(trim(cname)//': In CreateExpMesh, elemConn=', &
    !          ESMF_LOGMSG_INFO)
    !      do i = 1,(3*NTRI)
    !        write(msg,*) trim(cname)//': ',i, &
    !          ' ',elemConn(i)
    !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
    !      enddo
 
    expmesh = esmf_meshcreate( parametricdim=2,spatialdim=2, &
         nodeids=nodeids, nodecoords=nodecoords, &
         nodeowners=nodeowners, elementids=elemids,&
         elementtypes=elemtypes, elementconn=elemconn, &
         rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    deallocate(nodeids)
    deallocate(nodecoords)
    deallocate(nodeowners)
    deallocate(elemids)
    deallocate(elemtypes)
    deallocate(elemconn)
    !
    !     Set flag to indicate that this processor has local export grid storage
    !
    !AW      if (lpet == 0) then
    !AW         call ESMF_GridGet( expMesh, localDECount=ldecnt, rc=rc )
    !AW         if (ESMF_LogFoundError(rc, PASSTHRU)) return
    !AW         expGridIsLocal = ldecnt.gt.0
    !AW      endif
 
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Create ESMF grid with arbitrary domain decomposition to match
    !     the native domain decomposition of the non-excluded points
    !     Note that the native grid layout is dim1=Y, dim2=X
    !     Note that coordinates and mask are not needed since this
    !     grid is only used to define fields for a redist operation
    !
    ! 3.a Set flag to indicate that this processor has local native grid storage
    !
    natgridislocal = iaproc .gt. 0 .and. iaproc .le. naproc
 
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      !
      ! 3.b Setup arbitrary sequence index list
      !
      do ipass = 1,2
        if (ipass.eq.2) then
          allocate (arbindexlist(arbindexcount,2), stat=rc)
          if (esmf_logfoundallocerror(rc, passthru)) return
        endif
        arbindexcount = 0
        ! list local native grid non-excluded points
        if ( natgridislocal ) then
          do jsea = 1,nseal
#ifdef W3_DIST
            isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
            isea = jsea
#endif
            arbindexcount = arbindexcount+1
            if (ipass.eq.2) then
              ix = mapsf(isea,1)
              iy = mapsf(isea,2)
              ! native grid layout: dim1=Y, dim2=X
              arbindexlist(arbindexcount,1) = iy
              arbindexlist(arbindexcount,2) = ix
            endif
          enddo
        endif
      enddo !ipass
      !
      ! 3.c Create ESMF native grid
      !
      natgrid = esmf_gridcreatenoperidim( &
           minindex=(/ 1, 1/), &
           maxindex=(/ny,nx/), &
           coorddep1=(/esmf_dim_arb,esmf_dim_arb/), &
           coorddep2=(/esmf_dim_arb,esmf_dim_arb/), &
           arbindexcount=arbindexcount, &
           arbindexlist=arbindexlist, &
           coordtypekind=esmf_typekind_rx, &
           coordsys=esmf_coordsys_sph_deg, &
           name=trim(cname)//"_native_grid", rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      !
      ! 3.d Deallocate arbitrary sequence index list
      !
      deallocate (arbindexlist, stat=rc)
      if (esmf_logfounddeallocerror(rc, passthru)) return
      !
      ! -------------------------------------------------------------------- /
      ! 4.  Create route handle for redist between native grid domain
      !     decomposition and the export grid domain decomposition
      !
      ! 4.a Create temporary fields
      !
      nfield = esmf_fieldcreate( natgrid, natarrayspec1d, &
           staggerloc=natstaggerloc, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      efield = esmf_fieldcreate(expmesh, typekind=esmf_typekind_rx, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      !
      ! 4.b Store route handle
      !
      call esmf_fieldrediststore( nfield, efield, n2erh, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      !
      ! 4.c Clean up
      !
      call esmf_fielddestroy( nfield, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fielddestroy( efield, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
#ifdef W3_PDLIB
    endif
#endif
 
    call esmf_logwrite(trim(cname)//': In CreateExpMesh, created expMesh', &
         esmf_logmsg_info)
 
    rc = esmf_success
    if (verbosity.gt.0) call esmf_logwrite(trim(cname)// &
         ': leaving CreateExpMesh', esmf_logmsg_info)
    !/
    !/ End of CreateExpMesh ---------------------------------------------- /
    !/
  end subroutine createexpmesh
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "SetupImpBmsk"
 
  subroutine setupimpbmsk( bmskField, impField, missingVal, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Setup background blending mask field for an import field
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       bmskField  Type   I/O blending mask field
    !       impField   Type   I   import field
    !       missingVal Real   I   missing value
    !       rc         Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field)      :: bmskField
    type(esmf_field)      :: impField
    real(ESMF_KIND_RX)    :: missingVal
    integer, optional     :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real(ESMF_KIND_RX), parameter :: zero = 0.0
    real(ESMF_KIND_RX), parameter :: half = 0.5
    real(ESMF_KIND_RX), parameter ::  one = 1.0
    integer, parameter :: nsig = imphalowidth-1
    integer, parameter :: niter = 10
    integer, parameter :: iter0 = 1-niter
    integer, parameter :: lde = 0
    integer :: iter, i, j, ii, jj, k, l
    integer :: elb(2), eub(2)
    integer :: tlb(2), tub(2)
    real(ESMF_KIND_RX), pointer :: mptr(:,:), dptr(:,:)
    type(esmf_field) :: cmskField
    real(ESMF_KIND_RX), pointer :: bmsk(:,:), cmsk(:,:)
    real(ESMF_KIND_RX) :: bsum, wsum
    real(ESMF_KIND_RX) :: wflt(-nsig:nsig,-nsig:nsig)
    character(ESMF_MAXSTR) :: fnm
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETUPIMPBMSK)
    integer :: timeSlice
#endif
    !
    ! -------------------------------------------------------------------- /
    ! Initialize filter
    !
    if (present(rc)) rc = esmf_success
 
    do l = -nsig,nsig
      do k = -nsig,nsig
        wflt(k,l) = exp( -half*( real(k,esmf_kind_rx)**2 &
             + real(l,esmf_kind_rx)**2 ) )
      enddo
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! Set up fields and pointers
    !
    call esmf_fieldget( impfield, name=fnm, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    cmskfield = esmf_fieldcreate( impgrid, imparrayspec2d, &
         totallwidth=imphalolwidth, totaluwidth=imphalouwidth, &
         staggerloc=impstaggerloc, indexflag=impindexflag, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( impgridislocal ) then
      call esmf_fieldget( impmask, localde=lde, farrayptr=mptr, &
           exclusivelbound=elb, exclusiveubound=eub, &
           totallbound=tlb, totalubound=tub, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( impfield,  localde=lde, farrayptr=dptr, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( bmskfield, localde=lde, farrayptr=bmsk, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( cmskfield, localde=lde, farrayptr=cmsk, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Create blending mask
    !
    if ( impgridislocal ) then
      do j = elb(2),eub(2)
        do i = elb(1),eub(1)
          if ( nint(dptr(i,j)).eq.nint(missingval) ) then
            bmsk(i,j) = zero
          else
            bmsk(i,j) = one
          endif
          cmsk(i,j) = bmsk(i,j)
        enddo
      enddo
    endif
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETUPIMPBMSK)
    timeslice = 1
    call esmf_fieldwrite( bmskfield, &
         "wmesmfmd_setupimpbmsk_"//trim(fnm)//".nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
#endif
 
    iter_loop: do iter = iter0,niter
 
      call esmf_fieldhalo( bmskfield, imphalorh, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      if ( impgridislocal ) then
 
        j_loop: do j = elb(2),eub(2)
          i_loop: do i = elb(1),eub(1)
            if ( nint(mptr(i,j)).eq.maskvalueland ) cycle i_loop
            if ( nint(dptr(i,j)).eq.nint(missingval) ) cycle i_loop
 
            if (iter.le.0) then
 
              ! initialize blending zone to zero
              l_loop0: do l = -1,1
                jj = j + l
                if ( jj.lt.tlb(2).or.jj.gt.tub(2) ) cycle l_loop0
                k_loop0: do k = -1,1
                  ii = i + k
                  if ( ii.lt.tlb(1).or.ii.gt.tub(1) ) cycle k_loop0
                  if ( nint(mptr(ii,jj)).eq.maskvalueland ) cycle k_loop0
                  if ( bmsk(ii,jj).eq.zero ) cmsk(i,j) = zero
                enddo k_loop0
              enddo l_loop0
 
            else
 
              ! iterate filter over blending zone
              bsum = zero
              wsum = zero
              l_loop: do l = -nsig,nsig
                jj = j + l
                if ( jj.lt.tlb(2).or.jj.gt.tub(2) ) cycle l_loop
                k_loop: do k = -nsig,nsig
                  ii = i + k
                  if ( ii.lt.tlb(1).or.ii.gt.tub(1) ) cycle k_loop
                  if ( nint(mptr(ii,jj)).eq.maskvalueland ) cycle k_loop
                  bsum = bsum + wflt(k,l)*bmsk(ii,jj)
                  wsum = wsum + wflt(k,l)
                enddo k_loop
              enddo l_loop
              if ( wsum.gt.zero ) cmsk(i,j) = bsum/wsum
 
            endif
 
          enddo i_loop
        enddo j_loop
 
        do j = elb(2),eub(2)
          do i = elb(1),eub(1)
            if ( nint(mptr(i,j)).eq.maskvalueland ) cycle
            bmsk(i,j) = cmsk(i,j)
          enddo
        enddo
 
      endif
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETUPIMPBMSK)
      timeslice = timeslice + 1
      call esmf_fieldwrite( bmskfield, &
           "wmesmfmd_setupimpbmsk_"//trim(fnm)//".nc", &
           overwrite=.true., timeslice=timeslice, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
#endif
 
    enddo iter_loop
    !
    ! -------------------------------------------------------------------- /
    ! Clean up
    !
    call esmf_fielddestroy( cmskfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    !/
    !/ End of SetupImpBmsk ----------------------------------------------- /
    !/
  end subroutine setupimpbmsk
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "BlendImpField"
 
  subroutine blendimpfield( impField, mbgField, bmskField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Blend import field with background field
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       impField   Type   I/O import field
    !       mbgField   Type   I   import background field
    !       bmskField  Type   I   blending mask field
    !       rc         Int    I/O Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field), intent(inout)  :: impField
    type(esmf_field), intent(in)     :: mbgField
    type(esmf_field), intent(in)     :: bmskField
    integer, optional, intent(inout) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real(ESMF_KIND_RX), parameter ::  one = 1.0
    integer, parameter :: lde = 0
    integer :: i, j
    integer :: elb(2), eub(2)
    real(ESMF_KIND_RX), pointer :: mptr(:,:), dptr(:,:), sptr(:,:)
    real(ESMF_KIND_RX), pointer :: bmsk(:,:)
    !
    ! -------------------------------------------------------------------- /
    ! Get array pointers and bounds
    !
    if (present(rc)) rc = esmf_success
 
    if ( impgridislocal ) then
      call esmf_fieldget( impmask, localde=lde, farrayptr=mptr, &
           exclusivelbound=elb, exclusiveubound=eub, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( impfield,  localde=lde, farrayptr=dptr, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( mbgfield,  localde=lde, farrayptr=sptr, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( bmskfield, localde=lde, farrayptr=bmsk, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Blend Fields
    !
    if ( impgridislocal ) then
      do j = elb(2),eub(2)
        do i = elb(1),eub(1)
          if ( nint(mptr(i,j)).eq.maskvalueland ) cycle
          dptr(i,j) = bmsk(i,j)*dptr(i,j) + (one-bmsk(i,j))*sptr(i,j)
        enddo
      enddo
    endif
    !/
    !/ End of BlendImpField ---------------------------------------------- /
    !/
  end subroutine blendimpfield
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "SetupImpMmsk"
 
  subroutine setupimpmmsk( mmskField, impField, fillVal, mskCreated, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           U. Turuncoglu           |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         18-May-2021 |
    !/                  +-----------------------------------+
    !/
    !/    18-May-2021 : Origination.                        ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Setup merging mask field for an import field for the cases that
    !     model domains does not overlap completely
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       mmskField   Type   I/O merging mask field
    !       impField    Type   I   import field
    !       fillVal     Real   I   fill value
    !       mskCreated  Log.   I/O mask is created or not
    !       rc          Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field)   :: mmskField
    type(esmf_field)   :: impField
    real(ESMF_KIND_RX) :: fillVal
    logical            :: mskCreated
    integer, optional  :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    integer, parameter :: lde = 0
    integer :: i, j
    integer :: elb(2), eub(2)
    integer :: tlb(2), tub(2)
    real(ESMF_KIND_RX), pointer :: mptr(:,:), dptr(:,:)
    real(ESMF_KIND_RX), pointer :: mmsk(:,:)
    character(ESMF_MAXSTR) :: fnm
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETUPIMPMMSK)
    integer, save :: timeSlice = 1
#endif
    !
    ! -------------------------------------------------------------------- /
    ! Check mask is created or not
    !
    if ( mskcreated ) return
    !
    ! -------------------------------------------------------------------- /
    ! Set up fields and pointers
    !
 
    call esmf_fieldget( impfield, name=fnm, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( impgridislocal ) then
      call esmf_fieldget( impmask, localde=lde, farrayptr=mptr, &
           exclusivelbound=elb, exclusiveubound=eub, &
           totallbound=tlb, totalubound=tub, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( impfield,  localde=lde, farrayptr=dptr, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( mmskfield, localde=lde, farrayptr=mmsk, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Create merging mask
    !
    if ( impgridislocal ) then
      do j = elb(2),eub(2)
        do i = elb(1),eub(1)
          mmsk(i,j) = 0.0
          if (dptr(i,j).lt.fillval) then
            mmsk(i,j) = 1.0
          end if
        enddo
      enddo
    endif
    !
    ! -------------------------------------------------------------------- /
    ! Set mask created flag
    !
    mskcreated = .true.
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_SETUPIMPMMSK)
    call esmf_fieldwrite( mmskfield, &
         "wmesmfmd_setupimpmmsk_"//trim(fnm)//".nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice+1
#endif
    !/
    !/ End of SetupImpMmsk ----------------------------------------------- /
    !/
  end subroutine setupimpmmsk
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "FieldFill"
 
  subroutine fieldfill(field, fillVal, rc)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    09-Aug-2017 : Remove mask parameter.              ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Fill ESMF Field
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       field      Type   I/O ESMF field
    !       fillVal    Real   I   fill value
    !       rc         Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field)                 :: field
    real(ESMF_KIND_RX)               :: fillVal
    integer,               optional  :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    integer :: ldecnt, lde, i, j, k
    integer :: rank
    integer :: lb1(1), ub1(1)
    integer :: lb2(2), ub2(2)
    integer :: lb3(3), ub3(3)
    real(ESMF_KIND_RX), pointer :: dptr1(:)
    real(ESMF_KIND_RX), pointer :: dptr2(:,:)
    real(ESMF_KIND_RX), pointer :: dptr3(:,:,:)
    integer, parameter :: iwt=10
    real(8) :: wstime, wftime
    !
    ! -------------------------------------------------------------------- /
    ! Fill Field
    !
    if (present(rc)) rc = esmf_success
 
    call esmf_vmwtime(wstime)
 
    call esmf_fieldget(field, localdecount=ldecnt, rank=rank, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return ! bail out
    if (rank.ne.1.and.rank.ne.2.and.rank.ne.3) then
      call esmf_logseterror(esmf_failure, rctoreturn=rc, &
           msg='FieldFill: rank must be 1, 2 or 3')
      return ! bail out
    endif
 
    do lde=0,ldecnt-1
 
      if (rank.eq.1) then
        call esmf_fieldget(field, localde=lde, farrayptr=dptr1, &
             exclusivelbound=lb1, exclusiveubound=ub1, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return ! bail out
      elseif (rank.eq.2) then
        call esmf_fieldget(field, localde=lde, farrayptr=dptr2, &
             exclusivelbound=lb2, exclusiveubound=ub2, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return ! bail out
      else
        call esmf_fieldget(field, localde=lde, farrayptr=dptr3, &
             exclusivelbound=lb3, exclusiveubound=ub3, rc=rc)
        if (esmf_logfounderror(rc, passthru)) return ! bail out
      endif
 
      if (rank.eq.1) then
        dptr1(lb1(1):ub1(1)) = fillval
      elseif (rank.eq.2) then
        dptr2(lb2(1):ub2(1),lb2(2):ub2(2)) = fillval
      else
        dptr3(lb3(1):ub3(1),lb3(2):ub3(2),lb3(3):ub3(3)) = fillval
      endif
 
    enddo
 
    call esmf_vmwtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
    !/
    !/ End of FieldFill ------------------------------------------------- /
    !/
  end subroutine fieldfill
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "FieldGather"
 
  subroutine fieldgather(field, n1, n2, fout, rc)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |     A. J. van der Westhuysen      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/    27-Feb-2018 : Modification for use with UNGTYPE   ( version 6.06 )
    !/
    !  1. Purpose :
    !
    !     All gather of ESMF field
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       field      Type   I/O ESMF field
    !       n1,n2      Int    I   Dimensions of output array
    !       fout       R.A.   O   global output array
    !       rc         Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_PDLIB
    use yownodepool, only: np, iplg
    use yowrankmodule, only: rank
#endif
    !/
    implicit none
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    type(esmf_field)                 :: field
    integer                          :: n1, n2
    real                             :: fout(n1,n2)
    integer,               optional  :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    character(500) :: msg
    integer :: i, j, k, ir, ip, count
    real(ESMF_KIND_RX) :: floc(n1,n2)
    real(ESMF_KIND_RX) :: floc1d(n1), floc1dary(n1*n2)
#ifdef W3_PDLIB
    real(ESMF_KIND_R8), pointer :: fldptr(:)
#endif
    integer, parameter :: iwt=9
    real(8) :: wstime, wftime
    !
    ! -------------------------------------------------------------------- /
    ! Gather Field
    !
    if (present(rc)) rc = esmf_success
 
    call esmf_vmwtime(wstime)
 
    if ( (gtype.eq.rlgtype).or.(gtype.eq.clgtype) ) then
      count = n1 * n2
      floc = 0.
      floc1dary = 0.
      call esmf_fieldgather( field,  floc, rootpet=0, vm=vm, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      do j=1,n2
        do i=1,n1
          floc1dary(i+(j-1)*n1) = floc(i,j)
        enddo
      enddo
      call esmf_vmbroadcast( vm, bcstdata=floc1dary, count=count, rootpet=0, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      do j=1,n2
        do i=1,n1
          fout(i,j) = floc1dary(i+(j-1)*n1)
        enddo
      enddo
    elseif (gtype.eq.ungtype) then
      count = n1
      floc1d = 0.
      call esmf_fieldgather( field,  floc1d, rootpet=0, vm=vm, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_vmbroadcast( vm, bcstdata=floc1d, count=count, rootpet=0, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
#ifdef W3_PDLIB
      if ( lpdlib .EQV. .false. ) then
#endif
        do k = 1, n1
          fout(k,1) = floc1d(k)
        enddo
#ifdef W3_PDLIB
      else
        count = 0
        do ir = 1, npet
          do ip = 1, rank(ir)%np
            count = count+1
            fout(rank(ir)%iplg(ip),1) = floc1d(count)
            !             write(msg,*) trim(cname)//': count,ir,ip =',count, &
            !                ir,ip
            !             call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
          enddo
        enddo
      endif
#endif
 
#ifdef W3_PDLIB
      !      call ESMF_LogWrite(trim(cname)//': In FieldGather, fout(k,1)=', &
      !          ESMF_LOGMSG_INFO)
      !      do k = 1, n1
      !         write(msg,*) trim(cname)//': fout(k,1) =',k, &
      !          ' ',fout(k,1)
      !         call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
      !      enddo
#endif
 
    endif
 
    call esmf_vmwtime(wftime)
    wtime(iwt) = wtime(iwt) + wftime - wstime
    wtcnt(iwt) = wtcnt(iwt) + 1
 
    !/
    !/ End of FieldGather ------------------------------------------------ /
    !/
  end subroutine fieldgather
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "FieldIndex"
 
  function fieldindex ( fnameList, fname, rc ) result (indx)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Return index associated with field name
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       fnameList StrA   I   Array of field names
    !       fname     Str    I   Field name
    !       rc        Int.   O   Return code
    !       indx      Int    I   Returned index of fname
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    character (len=*) :: fnamelist(:)
    character (len=*) :: fname
    integer           :: rc
    integer           :: indx
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    integer       :: i, check
    !
    ! -------------------------------------------------------------------- /
    ! Find name in fnameList that matches fname
    !
    check = lbound(fnamelist,1) - 1
    indx = check
    do i = lbound(fnamelist,1),ubound(fnamelist,1)
      if ( trim(fnamelist(i)).eq.trim(fname) ) then
        indx = i
        exit
      endif
    enddo
    if ( indx.eq.check ) then
      call esmf_logseterror(esmf_failure, rctoreturn=rc, &
           msg='FieldIndex: input name ('//fname//') not in list')
    endif
    !/
    !/ End of FieldIndex ------------------------------------------------- /
    !/
  end function fieldindex
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "PrintTimers"
 
  subroutine printtimers ( cname, wtnam, wtcnt, wtime )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Jan-2017 |
    !/                  +-----------------------------------+
    !/
    !/    20-Jan-2017 : Origination.                        ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Print wallclock timers to ESMF log file
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       cname     Str    I   Name of component
    !       wtnam     Str    I   Timer names
    !       wtcnt     Int    I   Timer counts
    !       wtime     R8     I   Timers
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    character(*) :: cname
    character(*) :: wtnam(:)
    integer      :: wtcnt(:)
    real(8)      :: wtime(:)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(128) :: msg
    integer :: k
    !
    ! -------------------------------------------------------------------- /
    ! Print timers to ESMF log file
    !
    write(msg,1) trim(cname),"timer","count","time"
    call esmf_logwrite(trim(msg), esmf_logmsg_info)
    do k=lbound(wtcnt,1),ubound(wtcnt,1)
      write(msg,2) trim(cname),trim(wtnam(k)),wtcnt(k),wtime(k)
      call esmf_logwrite(trim(msg), esmf_logmsg_info)
    enddo
    !
    ! -------------------------------------------------------------------- /
    ! Formats
    !
1   format(a,': wtime: ',a20,a10,a14)
2   format(a,': wtime: ',a20,i10,e14.6)
    !/
    !/ End of PrintTimers ------------------------------------------------ /
    !/
  end subroutine printtimers
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcDecomp"
 
  subroutine calcdecomp ( nx, ny, nproc, npmin, adjust, nxproc, nyproc, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Calculate a 2D processor layout
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       nx,ny     Int    I   Grid dimensions
    !       nproc     Int    I   Total processor count
    !       npmin     Int    I   Min number of grid points per tile per direction
    !       adjust    Log    I   Enable/disable adjusting proc count downward
    !       nxproc    Int    O   Processor count in x-direction
    !       nyproc    Int    O   Processor count in y-direction
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    integer, intent(in)    :: nx, ny
    integer, intent(in)    :: nproc
    integer, intent(in)    :: npmin
    logical, intent(in)    :: adjust
    integer, intent(out)   :: nxproc
    integer, intent(out)   :: nyproc
    integer, intent(inout) :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    integer, parameter :: k = 4
    integer :: mproc, n, nfac, irp
    real(k) :: gr, rp, pr, diff, npx, npy
    character(256) :: msg
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    if ( nx.gt.ny ) then
      gr = real(nx,k)/real(ny,k)
    else
      gr = real(ny,k)/real(nx,k)
    endif
 
    mproc = nproc
    mproc_loop: do
 
      irp = int(sqrt(real(mproc,k)))
      diff = huge(gr)
      nfac = mproc
      do n = irp,mproc
        if ( mod(mproc,n).ne.0 ) cycle
        pr = real(n**2,k)/real(mproc,k)
        if ( abs(gr-pr) < diff ) then
          diff = abs(gr-pr)
          nfac = n
        endif
      enddo
      if ( nx.gt.ny ) then
        nxproc = nfac
        nyproc = mproc/nfac
      else
        nxproc = mproc/nfac
        nyproc = nfac
      endif
 
      npx = nx/real(nxproc,k)
      npy = ny/real(nyproc,k)
      if (.not.adjust) exit mproc_loop
 
      if ( npx.ge.npmin .and. npy.ge.npmin ) then
        exit mproc_loop
      else
        if ( mproc.gt.1 ) then
          mproc = mproc - 1
        else
          exit mproc_loop
        endif
      endif
 
    enddo mproc_loop
 
    if ( npx.lt.npmin .or. npy.lt.npmin ) then
      write(msg,'(a,7i6)') 'proc count is too large for grid size:', &
           nx,ny,npmin,nproc,mproc,nxproc,nyproc
      call esmf_logwrite(trim(msg), esmf_logmsg_error)
      rc = esmf_failure
    endif
    !/
    !/ End of CalcDecomp ------------------------------------------------- /
    !/
  end subroutine calcdecomp
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "GetEnvValue"
 
  subroutine getenvvalue ( cenv, cval, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Get value of environment variable
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       cenv      Str    I   Name of environment variable
    !       cval      Str    O   Value of environment variable
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    character(*) :: cenv
    character(*) :: cval
    integer      :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(256) :: msg
    integer :: length, istat
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
    call get_environment_variable( name=trim(cenv), value=cval, &
         length=length, trim_name=.false., status=istat )
    if (istat.lt.0) then
      ! The VALUE argument is present and has a length less than
      ! the significant length of the environment variable value.
      write(msg,'(a,i3,a)') "Length of input variable", &
           " is less than length of environment variable " &
           //trim(cenv)//" value (",length,")."
      call esmf_logwrite(trim(msg), esmf_logmsg_error)
      rc = istat
    elseif (istat.gt.0) then
      ! 1: The specified environment variable NAME does not exist.
      ! 2: The processor does not support environment variables.
      cval=" "
    endif
    if (length.eq.0) cval=" "
    !/
    !/ End of GetEnvValue ------------------------------------------------ /
    !/
  end subroutine getenvvalue
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "GetZlevels"
 
  subroutine getzlevels ( rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Get array of z-levels from zlfile for SDC
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    integer      :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(256) :: msg
    integer :: k, iunit, ierr
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    if (len_trim(zlfile).eq.0 .or. trim(zlfile) .eq. 'none') then
 
      nz = 1
      allocate(zl(nz), stat=rc)
      if (esmf_logfoundallocerror(rc, passthru)) return
      zl(1) = 0
 
    else
 
      call esmf_utiliounitget(iunit, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      open(unit=iunit, file=trim(zlfile), form='formatted', &
           status='old', access='sequential', iostat=ierr)
      if (ierr.ne.0) then
        msg = "failed opening "//trim(zlfile)
        call esmf_logwrite(trim(msg), esmf_logmsg_error)
        rc = esmf_failure
        return
      endif
      read(iunit, fmt=*, iostat=ierr) nz
      if (ierr.ne.0) then
        msg = "read nz failed: "//trim(zlfile)
        call esmf_logwrite(trim(msg), esmf_logmsg_error)
        rc = esmf_failure
        return
      endif
      allocate(zl(nz), stat=rc)
      if (esmf_logfoundallocerror(rc, passthru)) return
      do k=1,nz
        read(iunit, fmt=*, iostat=ierr) zl(k)
        if (ierr.ne.0) then
          msg = "read zl failed: "//trim(zlfile)
          call esmf_logwrite(trim(msg), esmf_logmsg_error)
          rc = esmf_failure
          return
        endif
      enddo
      close(iunit)
 
    endif
    !/
    !/ End of GetZlevels ------------------------------------------------- /
    !/
  end subroutine getzlevels
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcCharnk"
 
  subroutine calccharnk ( chkField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Calculate Charnock for export
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       chkField  Type   I/O 2D Charnock export field
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field) :: chkfield
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real   , parameter :: zero  = 0.0
    logical, save :: firstcall = .true.
    integer :: isea, jsea
    real    :: emean, fmean, fmean1, wnmean, amax, ustar, ustdr, &
         tauwx, tauwy, cd, z0, fmeanws, dlwmean
    logical :: llws(nspec)
    type(esmf_field) :: chknfield
    real(esmf_kind_rx), pointer :: chkn(:)
    integer, save :: timeslice = 1
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    chknfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call fieldfill( chknfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( natgridislocal ) then
 
      call esmf_fieldget( chknfield, farrayptr=chkn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
        isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
        isea = jsea
#endif
        if ( firstcall ) then
          charn(jsea) = zero
#ifdef W3_ST3
          llws(:) = .true.
          ustar = zero
          ustdr = zero
          tauwx = zero
          tauwy = zero
          call w3spr3( va(:,jsea), cg(1:nk,isea), wn(1:nk,isea),   &
               emean, fmean, fmean1, wnmean, amax,         &
               u10(isea), u10d(isea), ustar, ustdr, tauwx, &
               tauwy, cd, z0, charn(jsea), llws, fmeanws )
#endif
#ifdef W3_ST4
          llws(:) = .true.
          ustar = zero
          ustdr = zero
          tauwx = zero
          tauwy = zero
          call w3spr4( va(:,jsea), cg(1:nk,isea), wn(1:nk,isea),   &
               emean, fmean, fmean1, wnmean, amax,         &
               u10(isea), u10d(isea), ustar, ustdr, tauwx, &
               tauwy, cd, z0, charn(jsea), llws, fmeanws,  &
               dlwmean )
#endif
        endif !firstCall
        chkn(jsea) = charn(jsea)
      enddo jsea_loop
 
    endif !natGridIsLocal
 
    call esmf_fieldredist( chknfield, chkfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call esmf_fielddestroy( chknfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    firstcall = .false.
 
#ifdef TEST_WMESMFMD_CHARNK
    call esmf_fieldwrite( chkfield, "wmesmfmd_charnk_chk.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcCharnk ------------------------------------------------- /
    !/
  end subroutine calccharnk
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcRoughl"
 
  subroutine calcroughl ( wrlField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Calculate 2D wave roughness length for export
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       wrlField  Type   I/O 2D roughness length export field
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    type(esmf_field) :: wrlfield
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real   , parameter :: zero  = 0.0
    logical, save :: firstcall = .true.
    integer :: isea, jsea, ix, iy
    real    :: emean, fmean, fmean1, wnmean, amax, ustar, ustdr, &
         tauwx, tauwy, cd, z0, fmeanws, dlwmean
    logical :: llws(nspec)
    type(esmf_field) :: wrlnfield
    real(esmf_kind_rx), pointer :: wrln(:)
    integer, save :: timeslice = 1
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    wrlnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call fieldfill( wrlnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( natgridislocal ) then
 
      call esmf_fieldget( wrlnfield, farrayptr=wrln, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
        isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
        isea = jsea
#endif
        ix = mapsf(isea,1)
        iy = mapsf(isea,2)
        IF ( mapsta(iy,ix) .EQ. 1 ) THEN
          if ( firstcall ) then
            charn(jsea) = zero
#ifdef W3_ST3
            llws(:) = .true.
            ustar = zero
            ustdr = zero
            tauwx = zero
            tauwy = zero
            call w3spr3( va(:,jsea), cg(1:nk,isea), wn(1:nk,isea),   &
                 emean, fmean, fmean1, wnmean, amax,         &
                 u10(isea), u10d(isea), ustar, ustdr, tauwx, &
                 tauwy, cd, z0, charn(jsea), llws, fmeanws )
#endif
#ifdef W3_ST4
            llws(:) = .true.
            ustar = zero
            ustdr = zero
            tauwx = zero
            tauwy = zero
            call w3spr4( va(:,jsea), cg(1:nk,isea), wn(1:nk,isea),   &
                 emean, fmean, fmean1, wnmean, amax,         &
                 u10(isea), u10d(isea), ustar, ustdr, tauwx, &
                 tauwy, cd, z0, charn(jsea), llws, fmeanws,  &
                 dlwmean )
#endif
          endif !firstCall
          wrln(jsea) = charn(jsea)*ust(isea)**2/grav
        endif
      enddo jsea_loop
 
    endif !natGridIsLocal
 
    call esmf_fieldredist( wrlnfield, wrlfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call esmf_fielddestroy( wrlnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    firstcall = .false.
 
#ifdef TEST_WMESMFMD_ROUGHL
    call esmf_fieldwrite( wrlfield, "wmesmfmd_roughl_wrl.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcRoughl ------------------------------------------------- /
    !/
  end subroutine calcroughl
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcBotcur"
 
  subroutine calcbotcur ( a, wbxField, wbyField, wbpField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Calculate wave-bottom currents for export
    !
    !  2. Method :
    !
    !   > Madsen, O. S. (1994), ICCE
    !
    !                          //
    !   U_bot       = sqrt( 2*|| S_ub dsig dtheta )
    !   (magnitude)          //
    !
    !                            //
    !                           || S_ub*sin(theta) dsig dtheta
    !                          //
    !   phi_b       = arctan( ---------------------------------- )
    !   (direction)              //
    !                           || S_ub*cos(theta) dsig dtheta
    !                          //
    !
    !                              //
    !                             || S_ub dsig dtheta
    !                            //
    !   U_bp        =  2*pi * ----------------------------
    !   (period)                 //
    !                           || sig*S_ub dsig dtheta
    !                          //
    !
    ! Where:
    !   S_ub(theta,k) = near-bottom orbital velocity spectrum
    !                 = ( sig^2 / sinh^2(kD) ) *  E(theta,k) / Cg
    !                 = ( sig^3 / sinh^2(kD) ) * Ac(theta,k) / Cg
    !   Ac(theta,k) = wave action density
    !   D = depth
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       a         Real   I   Input spectra (in par list to change shape)
    !       wbxField  Type   I/O WBC 2D eastward-component export field
    !       wbyField  Type   I/O WBC 2D northward-component export field
    !       wbpField  Type   I/O WBC 2D period export field
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    real             :: a(nth,nk,0:nseal)
    type(esmf_field) :: wbxfield
    type(esmf_field) :: wbyfield
    type(esmf_field) :: wbpfield
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real(8), parameter :: zero  = 0.0
    real(8), parameter :: half  = 0.5
    real(8), parameter ::  one  = 1.0
    real(8), parameter ::  two  = 2.0
    ! kdmin = 1e-7: sinh(kdmin)**2 ~ 1e-14
    real(8), parameter :: kdmin = 1e-7
    ! kdmax = 18.0:  1/sinh(kdmax)**2 ~ 1e-14
    real(8), parameter :: kdmax = 18.0
    integer :: isea, jsea, ik, ith
    real(8) :: depth
    real(8) :: kd, fack, fkd, aka, akx, aky, abr, ubr, ubx, uby, dir
    real(8), allocatable :: sig2(:)
    type(esmf_field) :: wbxnfield, wbynfield, wbpnfield
    real(esmf_kind_rx), pointer :: wbxn(:), wbyn(:), wbpn(:)
    integer, save :: timeslice = 1
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    wbxnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    wbynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    wbpnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call fieldfill( wbxnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( wbynfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( wbpnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( natgridislocal ) then
 
      call esmf_fieldget( wbxnfield, farrayptr=wbxn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( wbynfield, farrayptr=wbyn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( wbpnfield, farrayptr=wbpn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      allocate( sig2(1:nk) )
      sig2(1:nk) = sig(1:nk)**2
 
      jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
        isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
        isea = jsea
#endif
        if ( dw(isea).le.zero ) cycle jsea_loop
        depth = max(dmin,dw(isea))
#ifdef USE_W3OUTG_FOR_EXPORT
        if ( aba(jsea).le.zero ) cycle jsea_loop
        if ( uba(jsea).le.zero ) cycle jsea_loop
        wbxn(jsea) = uba(jsea)*cos(ubd(jsea))
        wbyn(jsea) = uba(jsea)*sin(ubd(jsea))
        wbpn(jsea) = tpi*aba(jsea)/uba(jsea)
#else
        abr = zero
        ubr = zero
        ubx = zero
        uby = zero
        ik_loop: do ik = 1,nk
          aka = zero
          akx = zero
          aky = zero
          ith_loop: do ith = 1,nth
            aka = aka + a(ith,ik,jsea)
            akx = akx + a(ith,ik,jsea)*ecos(ith)
            aky = aky + a(ith,ik,jsea)*esin(ith)
          enddo ith_loop
          fack = dden(ik)/cg(ik,isea)
          kd = max(kdmin,min(kdmax,wn(ik,isea)*depth))
          fkd = fack/sinh(kd)**2
          abr = abr + aka*fkd
          ubr = ubr + aka*sig2(ik)*fkd
          ubx = ubx + akx*sig2(ik)*fkd
          uby = uby + aky*sig2(ik)*fkd
        enddo ik_loop
        if ( abr.le.zero .or. ubr.le.zero ) cycle jsea_loop
        abr = sqrt(two*abr)
        ubr = sqrt(two*ubr)
        dir = atan2(uby,ubx)
        wbxn(jsea) = ubr*cos(dir)
        wbyn(jsea) = ubr*sin(dir)
        wbpn(jsea) = tpi*abr/ubr
#endif
      enddo jsea_loop
 
      deallocate( sig2 )
 
    endif !natGridIsLocal
 
    call esmf_fieldredist( wbxnfield, wbxfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( wbynfield, wbyfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( wbpnfield, wbpfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call esmf_fielddestroy( wbxnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( wbynfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( wbpnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
#ifdef TEST_WMESMFMD_BOTCUR
    call esmf_fieldwrite( wbxfield, "wmesmfmd_botcur_wbx.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( wbyfield, "wmesmfmd_botcur_wby.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( wbpfield, "wmesmfmd_botcur_wbp.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcBotcur ------------------------------------------------- /
    !/
  end subroutine calcbotcur
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcRadstr2D"
 
  subroutine calcradstr2d ( a, sxxField, sxyField, syyField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |     A. J. van der Westhuysen      |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/    27-Feb-2018 : Modification for use with UNGTYPE   ( version 6.06 )
    !/
    !  1. Purpose :
    !
    !     Calculate 2D radiation stresses for export
    !
    !  2. Method :
    !
    ! Radiation stresses are defined as:
    !
    !                   //
    !   Sxx = rho grav || (N*cos^2(theta) + N - 1/2) * sig*Ac(theta,k)/Cg dsig dtheta
    !                 //
    !                   //
    !   Sxy = rho grav ||  N*sin(theta)*cos(theta)   * sig*Ac(theta,k)/Cg dsig dtheta
    !                 //
    !                   //
    !   Syy = rho grav || (N*sin^2(theta) + N - 1/2) * sig*Ac(theta,k)/Cg dsig dtheta
    !                 //
    !
    ! Where:
    !   rho = density of sea water
    !   grav = acceleration due to gravity
    !   Ac(theta,k) = wave action density
    !   N = Cg/C = ratio of group velocity and phase velocity
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       a         Real   I   Input spectra (in par list to change shape)
    !       sxxField  Type   I/O RS 2D eastward-component export field
    !       sxyField  Type   I/O RS 2D eastward-northward-component export field
    !       syyField  Type   I/O RS 2D northward-component export field
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_PDLIB
    use yownodepool, only: np, iplg
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    real             :: a(nth,nk,0:nseal)
    type(esmf_field) :: sxxfield
    type(esmf_field) :: sxyfield
    type(esmf_field) :: syyfield
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    character(ESMF_MAXSTR) :: cname
    character(128) :: msg
    real(8), parameter :: zero  = 0.0
    real(8), parameter :: half  = 0.5
    real(8), parameter ::  one  = 1.0
    real(8), parameter ::  two  = 2.0
    integer :: isea, jsea, ik, ith
    real(8) :: sxxs, sxys, syys
    real(8) :: akxx, akxy, akyy, cgoc, facd, fack, facs
    type(esmf_field) :: sxxnfield, sxynfield, syynfield
    real(esmf_kind_rx), pointer :: sxxn(:), sxyn(:), syyn(:)
    integer, save :: timeslice = 1
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    !     For regular and curvilinear grids the native grid has a 2D
    !     layout, whereas for unstructured meshes it is a 1D array
    if ( (gtype.eq.rlgtype).or.(gtype.eq.clgtype) ) then
      sxxnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
           staggerloc=natstaggerloc, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      sxynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
           staggerloc=natstaggerloc, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      syynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
           staggerloc=natstaggerloc, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
    elseif (gtype.eq.ungtype) then
#ifdef W3_PDLIB
      if ( lpdlib .EQV. .false. ) then
#endif
        sxxnfield = esmf_fieldcreate( natgrid, natarrayspec1d, &
             staggerloc=natstaggerloc, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        sxynfield = esmf_fieldcreate( natgrid, natarrayspec1d, &
             staggerloc=natstaggerloc, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        syynfield = esmf_fieldcreate( natgrid, natarrayspec1d, &
             staggerloc=natstaggerloc, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
#ifdef W3_PDLIB
      endif
#endif
    endif
 
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      call fieldfill( sxxnfield, zerovalue, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call fieldfill( sxynfield, zerovalue, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call fieldfill( syynfield, zerovalue, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
#ifdef W3_PDLIB
    endif
#endif
 
    if ( natgridislocal ) then
 
#ifdef W3_PDLIB
      if ( lpdlib .EQV. .false. ) then
        !        Use auxiliary native grid/mesh to populate and redistribute data
#endif
        call esmf_fieldget( sxxnfield, farrayptr=sxxn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        call esmf_fieldget( sxynfield, farrayptr=sxyn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        call esmf_fieldget( syynfield, farrayptr=syyn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
#ifdef W3_PDLIB
      else
        !        Use single domain-decomposed native mesh to populate and communicate data
        call esmf_fieldget( sxxfield, farrayptr=sxxn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        call esmf_fieldget( sxyfield, farrayptr=sxyn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
        call esmf_fieldget( syyfield, farrayptr=syyn, rc=rc )
        if (esmf_logfounderror(rc, passthru)) return
      endif
#endif
 
      facd = dwat*grav
#ifdef W3_PDLIB
      if ( lpdlib .EQV. .false. ) then
#endif
        jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
          isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
          isea = jsea
#endif
          if ( dw(isea).le.zero ) cycle jsea_loop
#ifdef USE_W3OUTG_FOR_EXPORT
          sxxn(jsea) = sxx(jsea)
          sxyn(jsea) = sxy(jsea)
          syyn(jsea) = syy(jsea)
#else
          sxxs = zero
          sxys = zero
          syys = zero
          ik_loop: do ik = 1,nk
            akxx = zero
            akxy = zero
            akyy = zero
            cgoc = cg(ik,isea)*wn(ik,isea)/sig(ik)
            cgoc = min(one,max(half,cgoc))
            ith_loop: do ith = 1,nth
              akxx = akxx + (cgoc*(ec2(ith)+one)-half)*a(ith,ik,jsea)
              akxy = akxy + cgoc*esc(ith)*a(ith,ik,jsea)
              akyy = akyy + (cgoc*(es2(ith)+one)-half)*a(ith,ik,jsea)
            enddo ith_loop
            fack = dden(ik)/cg(ik,isea)
            sxxs = sxxs + akxx*fack
            sxys = sxys + akxy*fack
            syys = syys + akyy*fack
          enddo ik_loop
          facs = (one+fte/cg(nk,isea))*facd
          sxxn(jsea) = sxxs*facs
          sxyn(jsea) = sxys*facs
          syyn(jsea) = syys*facs
#endif
        enddo jsea_loop
#ifdef W3_PDLIB
      else
        jsea_loop2: do jsea = 1,np
          isea = iplg(jsea)
          !          if ( dw(isea).le.zero ) cycle jsea_loop
          sxxn(jsea) = sxx(jsea)
          sxyn(jsea) = sxy(jsea)
          syyn(jsea) = syy(jsea)
          !        write(msg,*) trim(cname)//' sxxn', sxxn(jsea)
          !        call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO)
        enddo jsea_loop2
      endif
#endif
 
    endif !natGridIsLocal
 
#ifdef W3_PDLIB
    if ( lpdlib .EQV. .false. ) then
#endif
      call esmf_fieldredist( sxxnfield, sxxfield, n2erh, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldredist( sxynfield, sxyfield, n2erh, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldredist( syynfield, syyfield, n2erh, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      call esmf_fielddestroy( sxxnfield, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fielddestroy( sxynfield, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fielddestroy( syynfield, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
#ifdef W3_PDLIB
    endif
#endif
 
#ifdef TEST_WMESMFMD_RADSTR2D
    call esmf_fieldwrite( sxxfield, "wmesmfmd_radstr2d_sxx.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( sxyfield, "wmesmfmd_radstr2d_sxy.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( syyfield, "wmesmfmd_radstr2d_syy.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcRadstr2D ----------------------------------------------- /
    !/
  end subroutine calcradstr2d
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcStokes3D"
 
  subroutine calcstokes3d ( a, usxField, usyField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |        T. J. Campbell, NRL        |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Aug-2017 |
    !/                  +-----------------------------------+
    !/
    !/    09-Aug-2017 : Origination.                        ( version 6.03 )
    !/
    !  1. Purpose :
    !
    !     Calculate 3D Stokes drift current for export
    !
    !  2. Method :
    !
    !    Kenyon, K. E. (1969), J. Geophys. R., Vol 74, No 28, p 6991
    !
    !    U_vec(z)
    !             //
    !      = 2 g || ( F(f,theta) k_vec/C cosh(2k(D+z))/sinh(2kD) ) dsig dtheta
    !           //
    !
    !          //
    !      =  || (Ac(k,theta) sig^2 k_vec/Cg cosh(2k(D+z))/sinh^2(kD) ) dsig dtheta
    !        //
    !
    ! Where:
    !   Ac(k,theta) = wave action density
    !   k_vec = k*[cos(theta),sin(theta)]
    !   D = depth
    !   z = height (0 = mean sea level)
    !
    ! In deep water (kD large):  cosh(2k(D+z))/sinh^2(kD) --> 2*exp(2kz)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       a         Real   I   Input spectra (in par list to change shape)
    !       usxField  Type   I/O 3D SDC eastward-component export field
    !       usyField  Type   I/O 3D SDC northward-component export field
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    implicit none
    real             :: a(nth,nk,0:nseal)
    type(esmf_field) :: usxField
    type(esmf_field) :: usyField
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    real(8), parameter :: zero  = 0.0
    real(8), parameter :: half  = 0.5
    real(8), parameter ::  one  = 1.0
    real(8), parameter ::  two  = 2.0
    ! kdmin = 1e-7: sinh(kdmin)**2 ~ 1e-14
    real(8), parameter :: kdmin = 1e-7
    ! kdmax = 18.0:  cosh(2*(kdmax+kz))/sinh(kdmax)**2 - 2*exp(2*kz) < 1e-14
    real(8), parameter :: kdmax = 18.0
    ! kdmin & kdmax settings used in w3iogomd
    real(8), parameter :: kdmin_us3d = 1e-3
    real(8), parameter :: kdmax_us3d = 6.0
    integer :: isea, jsea, ik, ith, iz
    real(8) :: depth
    real(8) :: akx, aky, kd, kz, fac1, fac2, fac3
    real(8) :: uzx(nz), uzy(nz)
    real(8), allocatable :: fack(:)
    type(esmf_field) :: usxnField, usynField
    real(ESMF_KIND_RX), pointer :: usxn(:,:), usyn(:,:)
    integer, save :: timeSlice = 1
    ! Need this workaround to deal with ESMF_FieldCreate not setting up the
    ! Fortran arrays with the ungridded dimension as the first array dimension
#define ESMF_ARBSEQ_WORKAROUND
#ifdef ESMF_ARBSEQ_WORKAROUND
    type(esmf_distgrid) :: natDistGrid
    type(esmf_array) :: usxnArray, usynArray
#endif
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
#ifdef ESMF_ARBSEQ_WORKAROUND
    call esmf_gridget( natgrid, distgrid=natdistgrid, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    usxnarray = esmf_arraycreate( natdistgrid, esmf_typekind_rx, &
         distgridtoarraymap=(/2/), undistlbound=(/1/), undistubound=(/nz/), rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    usxnfield = esmf_fieldcreate( natgrid, usxnarray, &
         gridtofieldmap=(/2,3/), ungriddedlbound=(/1/), ungriddedubound=(/nz/), &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    usynarray = esmf_arraycreate( natdistgrid, esmf_typekind_rx, &
         distgridtoarraymap=(/2/), undistlbound=(/1/), undistubound=(/nz/), rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    usynfield = esmf_fieldcreate( natgrid, usynarray, &
         gridtofieldmap=(/2,3/), ungriddedlbound=(/1/), ungriddedubound=(/nz/), &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
#else
    usxnfield = esmf_fieldcreate( natgrid, natarrayspec3d, &
         gridtofieldmap=(/2,3/), ungriddedlbound=(/1/), ungriddedubound=(/nz/), &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    usynfield = esmf_fieldcreate( natgrid, natarrayspec3d, &
         gridtofieldmap=(/2,3/), ungriddedlbound=(/1/), ungriddedubound=(/nz/), &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
#endif
 
    call fieldfill( usxnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( usynfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( natgridislocal ) then
 
      call esmf_fieldget( usxnfield, farrayptr=usxn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( usynfield, farrayptr=usyn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      allocate( fack(1:nk) )
      fack(1:nk) = dden(1:nk) * sig(1:nk)
 
      jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
        isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
        isea = jsea
#endif
        if ( dw(isea).le.zero ) cycle jsea_loop
        depth = max(dmin,dw(isea))
        uzx(:) = zero
        uzy(:) = zero
#ifdef USE_W3OUTG_FOR_EXPORT
        ik_loop: do ik = us3df(2),us3df(3)
          fac1 = tpiinv*dsii(ik)
          kd = max(kdmin_us3d,wn(ik,isea)*dw(isea))
          iz_loop: do iz = 1,nz
            if ( dw(isea)+zl(iz).le.zero ) cycle iz_loop
            kz = wn(ik,isea)*zl(iz)
            if ( kd .lt. kdmax_us3d ) then
              fac2 = fac1*cosh(two*max(zero,kd+kz))/cosh(two*kd)
            else
              fac2 = fac1*exp(two*kz)
            endif
            uzx(iz) = uzx(iz) + us3d(jsea,ik   )*fac2
            uzy(iz) = uzy(iz) + us3d(jsea,nk+ik)*fac2
          enddo iz_loop
        enddo ik_loop
#else
        ik_loop: do ik = 1,nk
          akx = zero
          aky = zero
          ith_loop: do ith = 1,nth
            akx = akx + a(ith,ik,jsea)*ecos(ith)
            aky = aky + a(ith,ik,jsea)*esin(ith)
          enddo ith_loop
          fac1 = fack(ik)*wn(ik,isea)/cg(ik,isea)
          kd = max(kdmin,wn(ik,isea)*depth)
          if ( kd .lt. kdmax ) then
            fac2 = fac1/sinh(kd)**2
          else
            fac2 = fac1*two
          endif
          akx = akx*fac2
          aky = aky*fac2
          iz_loop: do iz = 1,nz
            if ( depth+zl(iz).le.zero ) cycle iz_loop
            kz = wn(ik,isea)*zl(iz)
            if ( kd .lt. kdmax ) then
              fac3 = cosh(two*max(zero,kd+kz))
            else
              fac3 = exp(two*kz)
            endif
            uzx(iz) = uzx(iz) + akx*fac3
            uzy(iz) = uzy(iz) + aky*fac3
          enddo iz_loop
        enddo ik_loop
#endif
        usxn(:,jsea) = uzx(:)
        usyn(:,jsea) = uzy(:)
      enddo jsea_loop
 
      deallocate( fack )
 
    endif !natGridIsLocal
 
    call esmf_fieldredist( usxnfield, usxfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( usynfield, usyfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
#ifdef ESMF_ARBSEQ_WORKAROUND
    call esmf_arraydestroy( usxnarray, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_arraydestroy( usynarray, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
#endif
    call esmf_fielddestroy( usxnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( usynfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
#ifdef TEST_WMESMFMD_STOKES3D
    call esmf_fieldwrite( usxfield, "wmesmfmd_stokes3d_usx.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( usyfield, "wmesmfmd_stokes3d_usy.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcStokes3D ----------------------------------------------- /
    !/
  end subroutine calcstokes3d
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "CalcPStokes"
 
  subroutine calcpstokes ( a, p1xField, p1yField, p2xField,   &
       p2yField, p3xField, p3yField, rc )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           J. Meixner              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-Oct-2019 |
    !/                  +-----------------------------------+
    !/
    !/    DD-MMM-YYYY : Origination.                        ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Calculate partitioned Stokes drift for export
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       a         Real   I   Input spectra (in par list to change shape)
    !       p1Field   Type   I/O
    !       p2Field   Type   I/O
    !       p3Field   Type   I/O
    !       rc        Int    O   Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3adatmd, ONLY: ussp
    USE w3iogomd, ONLY: calc_u3stokes
    IMPLICIT NONE
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    real             :: a(nth,nk,0:nseal)
    type(esmf_field) :: p1xField,p2xField,p3xField
    type(esmf_field) :: p1yField,p2yField,p3yField
    integer          :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    !real(8) :: sxxs, sxys, syys
    type(esmf_field) :: p1xnField, p2xnField, p3xnField
    type(esmf_field) :: p1ynField, p2ynField, p3ynField
    real(ESMF_KIND_RX), pointer :: p1xn(:), p2xn(:), p3xn(:)
    real(ESMF_KIND_RX), pointer :: p1yn(:), p2yn(:), p3yn(:)
    integer, save :: timeSlice = 1
    integer :: isea,jsea
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
 
    p1xnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    p1ynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    p2xnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    p2ynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    p3xnfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    p3ynfield = esmf_fieldcreate( natgrid, natarrayspec2d, &
         staggerloc=natstaggerloc, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call fieldfill( p1xnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( p1ynfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( p2xnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( p2ynfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( p3xnfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call fieldfill( p3ynfield, zerovalue, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    if ( natgridislocal ) then
 
      call esmf_fieldget( p1xnfield, farrayptr=p1xn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( p1ynfield, farrayptr=p1yn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( p2xnfield, farrayptr=p2xn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( p2ynfield, farrayptr=p2yn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( p3xnfield, farrayptr=p3xn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
      call esmf_fieldget( p3ynfield, farrayptr=p3yn, rc=rc )
      if (esmf_logfounderror(rc, passthru)) return
 
      call calc_u3stokes ( a , 2 )
 
      jsea_loop: do jsea = 1,nseal
#ifdef W3_DIST
        isea = iaproc + (jsea-1)*naproc
#endif
#ifdef W3_SHRD
        isea = jsea
#endif
 
        p1xn(jsea)=ussp(jsea,1)
        p1yn(jsea)=ussp(jsea,nk+1)
        p2xn(jsea)=ussp(jsea,2)
        p2yn(jsea)=ussp(jsea,nk+2)
        p3xn(jsea)=ussp(jsea,3)
        p3yn(jsea)=ussp(jsea,nk+3)
      enddo jsea_loop
 
    endif !natGridIsLocal
 
    call esmf_fieldredist( p1xnfield, p1xfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( p1ynfield, p1yfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( p2xnfield, p2xfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( p2ynfield, p2yfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( p3xnfield, p3xfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldredist( p3ynfield, p3yfield, n2erh, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
    call esmf_fielddestroy( p1xnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( p2xnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( p3xnfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( p1ynfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( p2ynfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fielddestroy( p3ynfield, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
 
#ifdef TEST_WMESMFMD_PSTOKES
    call esmf_fieldwrite( p1xfield, "wmesmfmd_pstokes_1x.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( p1yfield, "wmesmfmd_pstokes_1y.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( p2xfield, "wmesmfmd_pstokes_2x.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( p2yfield, "wmesmfmd_pstokes_2y.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( p3xfield, "wmesmfmd_pstokes_3x.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( p3yfield, "wmesmfmd_pstokes_3y.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of CalcPStokes ----------------------------------------------- /
    !/
  end subroutine calcpstokes
  !/ ------------------------------------------------------------------- /
  !/
#undef METHOD
#define METHOD "ReadFromFile"
 
  subroutine readfromfile (idfld, fldwx, fldwy, time0, timen, rc)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           U. Turuncoglu           |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         18-May-2021 |
    !/                  +-----------------------------------+
    !/
    !/    18-May-2021 : Origination.                        ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Read input file to fill unmapped point for regional applications
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       idfld     Str    I/O Field name
    !       fldwx     Type   I/O 2D eastward-component of field
    !       fldwy     Type   I/O 2D northward-component of field
    !       time0     Int    I   Time stamp for current time
    !       timen     Int    I   Time stamp for end time
    !       rc        Int    I/O Return code
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      NONE
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE w3fldsmd, ONLY: w3fldo, w3fldg
    USE wmunitmd, ONLY: wmuget, wmuset
    IMPLICIT NONE
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    character(len=3), intent(inout)  :: idfld
    type(esmf_field), intent(inout)  :: fldwx
    type(esmf_field), intent(inout)  :: fldwy
    integer, intent(in)              :: time0(2)
    integer, intent(in)              :: timen(2)
    integer, intent(inout), optional :: rc
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    integer :: ierr, tw0l(2), twnl(2), lb(2), ub(2)
    real :: wx0l(nx,ny), wy0l(nx,ny)
    real :: wxnl(nx,ny), wynl(nx,ny)
    real :: dt0l(nx,ny), dtnl(nx,ny)
    real(ESMF_KIND_RX), pointer :: dptr(:,:)
    integer :: mdse = 6
    integer :: mdst = 10
    integer, save :: mdsf
    character(256) :: logmsg
    logical :: flagsc = .false.
    integer, parameter :: lde = 0
    logical, save :: firstCall = .true.
    character(len=13) :: tsstr
    character(len=3) :: tsfld
    integer :: nxt, nyt, gtypet, filler(3), tideflag
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_READFROMFILE)
    integer, save :: timeSlice = 1
#endif
    !
    ! -------------------------------------------------------------------- /
    !
    rc = esmf_success
 
    if (firstcall) then
      ! assign unit number for input file
      call wmuget(mdse, mdst, mdsf, 'INP')
      call wmuset(mdse, mdst, mdsf, .true., desc='Input data file')
 
      ! open file
      call w3fldo('READ', idfld, mdsf, mdst, mdse, nx, ny, gtype, ierr)
      if (ierr.ne.0) then
        write(logmsg,*) "Error in opening "//idfld//", iostat = ", ierr
        call esmf_logwrite(trim(logmsg), esmf_logmsg_error)
        rc = esmf_failure
        return
      endif
 
      firstcall = .false.
    end if
 
    ! init variables
    wx0l = 0.0
    wy0l = 0.0
    dt0l = 0.0
    wxnl = 0.0
    wynl = 0.0
    dtnl = 0.0
 
    ! need to rewind to the begining of the file to access
    ! data of requested date correctly
    rewind(mdsf)
 
    ! read header information
    ! this was inside of w3fldo call but since we are opening file
    ! once and rewinding, the header need to be read
    read(mdsf, iostat=ierr) tsstr, tsfld, nxt, nyt, &
         gtypet, filler(1:2), tideflag
 
    ! read input
    call w3fldg('READ', idfld, mdsf, mdst, mdse, nx, ny, &
         nx, ny, time0, timen, tw0l, wx0l, wy0l, dt0l, twnl, &
         wxnl, wynl, dtnl, ierr, flagsc)
 
    ! fill fields with data belong to current time
    if ( impgridislocal ) then
      call esmf_fieldget(fldwx, localde=lde, farrayptr=dptr, &
           exclusivelbound=lb, exclusiveubound=ub, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      dptr(lb(1):ub(1),lb(2):ub(2)) = wx0l(lb(1):ub(1),lb(2):ub(2))
      if (associated(dptr)) nullify(dptr)
      call esmf_fieldget(fldwy, localde=lde, farrayptr=dptr, &
           exclusivelbound=lb, exclusiveubound=ub, rc=rc)
      if (esmf_logfounderror(rc, passthru)) return
      dptr(lb(1):ub(1),lb(2):ub(2)) = wy0l(lb(1):ub(1),lb(2):ub(2))
      if (associated(dptr)) nullify(dptr)
    end if
 
#if defined(TEST_WMESMFMD) || defined(TEST_WMESMFMD_READFROMFILE)
    write(logmsg,*) 'time0 = ', time0(1), time0(2)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'timen = ', timen(1), timen(2)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'tw0 = ', tw0l(1), tw0l(2)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'twn = ', twnl(1), twnl(2)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'wx0 min, max = ', minval(wx0l), maxval(wx0l)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'wy0 min, max = ', minval(wy0l), maxval(wy0l)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'wxn min, max = ', minval(wxnl), maxval(wxnl)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    write(logmsg,*) 'wyn min, max = ', minval(wynl), maxval(wynl)
    call esmf_logwrite(trim(logmsg), esmf_logmsg_info, rc=rc)
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( fldwx, &
         "wmesmfmd_read_wx0.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    call esmf_fieldwrite( fldwy, &
         "wmesmfmd_read_wy0.nc", &
         overwrite=.true., timeslice=timeslice, rc=rc )
    if (esmf_logfounderror(rc, passthru)) return
    timeslice = timeslice + 1
#endif
    !/
    !/ End of ReadFromFile ------------------------------------------- /
    !/
  end subroutine readfromfile
  !/ ------------------------------------------------------------------- /
  !/
  !/ End of module WMESMFMD -------------------------------------------- /
  !/
end module wmesmfmd