Updates for wave coupling in NUOPC (NOAA-GFDL#23)

* Adding extra ensemble slot for waves.

* Updates for wave coupling
- Adding wave information to mech_forcing type to pass from ice_ocean_boundary to wave types
- This is only set-up to read surface Stokes drift and guess the wavelength as something reasonable for now to demonstrate that it works.  This needs to be set-up properly before merging this into the main repository.

* Updates to make Stokes drift from multiple bands work with the coupler approach

* Adding wave stokes drift import to nuopc cap

Co-authored-by: Brandon Reichl <Brandon.Reichl@noaa.gov>

config_src/coupled_driver/MOM_surface_forcing_gfdl.F90

@@ -189,13 +189,20 @@ module MOM_surface_forcing_gfdl
                                                             !! ice-shelves, expressed as a coefficient
                                                             !! for divergence damping, as determined
                                                             !! outside of the ocean model [m3 s-1]
+  real, pointer, dimension(:,:)   :: ustk0           => NULL() !<
+  real, pointer, dimension(:,:)   :: vstk0           => NULL() !<
+  real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !<
+  real, pointer, dimension(:,:,:) :: ustkb           => NULL() !<
+  real, pointer, dimension(:,:,:) :: vstkb           => NULL() !<
+
   integer :: xtype                    !< The type of the exchange - REGRID, REDIST or DIRECT
   type(coupler_2d_bc_type) :: fluxes  !< A structure that may contain an array of named fields
                                       !! used for passive tracer fluxes.
   integer :: wind_stagger = -999      !< A flag indicating the spatial discretization of wind stresses.
                                       !! This flag may be set by the flux-exchange code, based on what
                                       !! the sea-ice model is providing.  Otherwise, the value from
                                       !! the surface_forcing_CS is used.
+  integer :: num_stk_bands            !< Number of Stokes drift bands passed through the coupler
 end type ice_ocean_boundary_type
 
 integer :: id_clock_forcing !< A CPU time clock
@@ -275,7 +282,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   ! flux type has been used.
   if (fluxes%dt_buoy_accum < 0) then
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., &
-                               ustar=.true., press=.true.)
+                               ustar=.true., press=.true. )
 
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -669,7 +676,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   real :: mass_eff      ! effective mass of sea ice for rigidity [kg m-2]
   real :: wt1, wt2      ! Relative weights of previous and current values of ustar, ND.
 
-  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
+  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0, istk
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
   integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
 
@@ -710,6 +717,9 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
        (associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
     call allocate_mech_forcing(G, forces, iceberg=.true.)
 
+  if ( associated(IOB%ustk0) ) &
+    call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands)
+
   if (associated(IOB%ice_rigidity)) then
     rigidity_at_h(:,:) = 0.0
     call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
@@ -769,6 +779,19 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
       forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
     enddo ; enddo
   endif
+  forces%stk_wavenumbers(:) = IOB%stk_wavenumbers
+  do j=js,je; do i=is,ie
+     forces%ustk0(i,j) = IOB%ustk0(i-I0,j-J0) ! How to be careful here that the domains are right?
+     forces%vstk0(i,j) = IOB%vstk0(i-I0,j-J0)
+  enddo ; enddo
+  call pass_vector(forces%ustk0,forces%vstk0, G%domain )
+  do istk = 1,IOB%num_stk_bands
+    do j=js,je; do i=is,ie
+      forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
+      forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
+    enddo; enddo
+    call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
+  enddo
 
   ! Find the net mass source in the input forcing without other adjustments.
   if (CS%approx_net_mass_src .and. associated(forces%net_mass_src)) then

config_src/coupled_driver/ocean_model_MOM.F90

@@ -557,7 +557,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     ! For now, the waves are only updated on the thermodynamics steps, because that is where
     ! the wave intensities are actually used to drive mixing.  At some point, the wave updates
     ! might also need to become a part of the ocean dynamics, according to B. Reichl.
-    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves)
+    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
   endif
 
   if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model.

config_src/nuopc_driver/mom_cap.F90

@@ -606,6 +606,20 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   Ice_ocean_boundary%lrunoff         = 0.0
   Ice_ocean_boundary%frunoff         = 0.0
 
+  if (ocean_state%use_waves) then
+    Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands
+    allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec),         &
+               Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec),         &
+               Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
+               Ice_ocean_boundary% vstkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
+               Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands))
+    Ice_ocean_boundary%ustk0           = 0.0
+    Ice_ocean_boundary%vstk0           = 0.0
+    Ice_ocean_boundary%stk_wavenumbers = ocean_state%Waves%WaveNum_Cen
+    Ice_ocean_boundary%ustkb           = 0.0
+    Ice_ocean_boundary%vstkb           = 0.0
+  endif
+
   ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state
   call ESMF_GridCompSetInternalState(gcomp, ocean_internalstate, rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -649,6 +663,17 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   !These are not currently used and changing requires a nuopc dictionary change
   !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx"        , "will provide")
   !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx"       , "will provide")
+  if (ocean_state%use_waves) then
+    if (Ice_ocean_boundary%num_stk_bands > 3) then 
+      call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this")
+    endif      
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_1" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_1", "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_2" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_2", "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_3" , "will provide")
+    call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_3", "will provide")
+  endif 
 
   !--------- export fields -------------
   call fld_list_add(fldsFrOcn_num, fldsFrOcn, "ocean_mask"                 , "will provide")

config_src/nuopc_driver/mom_cap_methods.F90

@@ -73,6 +73,8 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   character(len=128)              :: fldname
   real(ESMF_KIND_R8), allocatable :: taux(:,:)
   real(ESMF_KIND_R8), allocatable :: tauy(:,:)
+  real(ESMF_KIND_R8), allocatable :: stkx1(:,:),stkx2(:,:),stkx3(:,:)
+  real(ESMF_KIND_R8), allocatable :: stky1(:,:),stky2(:,:),stky3(:,:)
   character(len=*)  , parameter   :: subname = '(mom_import)'
 
   rc = ESMF_SUCCESS
@@ -245,6 +247,56 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
        isc, iec, jsc, jec, ice_ocean_boundary%mi, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+
+  !----
+  ! Partitioned Stokes Drift Components 
+  !----
+  if ( associated(ice_ocean_boundary%ustkb) ) then    
+    allocate(stkx1(isc:iec,jsc:jec))
+    allocate(stky1(isc:iec,jsc:jec))
+    allocate(stkx2(isc:iec,jsc:jec))
+    allocate(stky2(isc:iec,jsc:jec))
+    allocate(stkx3(isc:iec,jsc:jec))
+    allocate(stky3(isc:iec,jsc:jec))
+
+    call state_getimport(importState,'eastward_partitioned_stokes_drift_1' , isc, iec, jsc, jec, stkx1,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call state_getimport(importState,'northward_partitioned_stokes_drift_1', isc, iec, jsc, jec, stky1,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call state_getimport(importState,'eastward_partitioned_stokes_drift_2' , isc, iec, jsc, jec, stkx2,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call state_getimport(importState,'northward_partitioned_stokes_drift_2', isc, iec, jsc, jec, stky2,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call state_getimport(importState,'eastward_partitioned_stokes_drift_3' , isc, iec, jsc, jec, stkx3,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call state_getimport(importState,'northward_partitioned_stokes_drift_3', isc, iec, jsc, jec, stky3,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    ! rotate from true zonal/meridional to local coordinates
+    do j = jsc, jec
+       jg = j + ocean_grid%jsc - jsc
+       do i = isc, iec
+          ig = i + ocean_grid%isc - isc
+          ice_ocean_boundary%ustkb(i,j,1) = ocean_grid%cos_rot(ig,jg)*stkx1(i,j) &
+               - ocean_grid%sin_rot(ig,jg)*stky1(i,j)
+          ice_ocean_boundary%vstkb(i,j,1) = ocean_grid%cos_rot(ig,jg)*stky1(i,j) &
+               + ocean_grid%sin_rot(ig,jg)*stkx1(i,j)
+
+          ice_ocean_boundary%ustkb(i,j,2) = ocean_grid%cos_rot(ig,jg)*stkx2(i,j) &
+               - ocean_grid%sin_rot(ig,jg)*stky2(i,j)
+          ice_ocean_boundary%vstkb(i,j,2) = ocean_grid%cos_rot(ig,jg)*stky2(i,j) &
+               + ocean_grid%sin_rot(ig,jg)*stkx2(i,j)
+
+          ice_ocean_boundary%ustkb(i,j,3) = ocean_grid%cos_rot(ig,jg)*stkx3(i,j) &
+               - ocean_grid%sin_rot(ig,jg)*stky3(i,j)
+          ice_ocean_boundary%vstkb(i,j,3) = ocean_grid%cos_rot(ig,jg)*stky3(i,j) &
+               + ocean_grid%sin_rot(ig,jg)*stkx3(i,j)
+       enddo
+    enddo
+
+    deallocate(stkx1,stkx2,stkx3,stky1,stky2,stky3)
+  endif
+
 end subroutine mom_import
 
 !> Maps outgoing ocean data to ESMF State

config_src/nuopc_driver/mom_ocean_model_nuopc.F90

@@ -143,7 +143,7 @@ module MOM_ocean_model_nuopc
 
   integer :: nstep = 0        !< The number of calls to update_ocean.
   logical :: use_ice_shelf    !< If true, the ice shelf model is enabled.
-  logical :: use_waves        !< If true use wave coupling.
+  logical,public :: use_waves !< If true use wave coupling.
 
   logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
                               !! ocean dynamics and forcing fluxes.
@@ -203,7 +203,7 @@ module MOM_ocean_model_nuopc
   type(marine_ice_CS), pointer :: &
     marine_ice_CSp => NULL()  !< A pointer to the control structure for the
                               !! marine ice effects module.
-  type(wave_parameters_cs), pointer :: &
+  type(wave_parameters_cs), pointer, public :: &
     Waves !< A structure containing pointers to the surface wave fields
   type(surface_forcing_CS), pointer :: &
     forcing_CSp => NULL()     !< A pointer to the MOM forcing control structure
@@ -386,6 +386,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
        "If true, enables surface wave modules.", default=.false.)
   if (OS%use_waves) then
     call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag)
+    call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS",OS%Waves%WaveNum_Cen, &
+           "Central wavenumbers for surface Stokes drift bands.",units='rad/m', &
+           default=0.12566)
   else
     call MOM_wave_interface_init_lite(param_file)
   endif
@@ -570,7 +573,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid, OS%US)
 
   if (OS%use_waves) then
-    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves)
+    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
   endif
 
   if (OS%nstep==0) then

config_src/nuopc_driver/mom_surface_forcing_nuopc.F90

@@ -181,9 +181,15 @@ module MOM_surface_forcing_nuopc
                                                               !! ice-shelves, expressed as a coefficient
                                                               !! for divergence damping, as determined
                                                               !! outside of the ocean model in [m3/s]
+  real, pointer, dimension(:,:)   :: ustk0           => NULL() !<
+  real, pointer, dimension(:,:)   :: vstk0           => NULL() !<
+  real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !<
+  real, pointer, dimension(:,:,:) :: ustkb           => NULL() !<
+  real, pointer, dimension(:,:,:) :: vstkb           => NULL() !<
+  integer :: num_stk_bands            !< Number of Stokes drift bands passed through the coupler
   integer :: xtype                                            !< The type of the exchange - REGRID, REDIST or DIRECT
   type(coupler_2d_bc_type)      :: fluxes                     !< A structure that may contain an array of
-                                                              !! named fields used for passive tracer fluxes.
+                                                              !! namedfields used for passive tracer fluxes.
   integer :: wind_stagger = -999                              !< A flag indicating the spatial discretization of
                                                               !! wind stresses.  This flag may be set by the
                                                               !! flux-exchange code, based on what the sea-ice
@@ -619,7 +625,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
   real :: mass_eff      !< effective mass of sea ice for rigidity (kg/m^2)
 
   integer :: wind_stagger  !< AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
-  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
+  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0, istk
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
   integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
 
@@ -658,6 +664,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
   if ( (associated(IOB%area_berg) .and. (.not. associated(forces%area_berg))) .or. &
        (associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
     call allocate_mech_forcing(G, forces, iceberg=.true.)
+
   if (associated(IOB%ice_rigidity)) then
     rigidity_at_h(:,:) = 0.0
     call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
@@ -668,6 +675,9 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
   if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
   if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
 
+  if ( associated(IOB%ustkb) ) &
+    call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands)
+
   ! applied surface pressure from atmosphere and cryosphere
   if (CS%use_limited_P_SSH) then
      forces%p_surf_SSH => forces%p_surf
@@ -825,6 +835,24 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
   endif   ! endif for wind related fields
 
+  ! wave to ocean coupling
+  if ( associated(IOB%ustkb) ) then
+
+    forces%stk_wavenumbers(:) = IOB%stk_wavenumbers
+    do j=js,je; do i=is,ie
+      forces%ustk0(i,j) = IOB%ustk0(i-I0,j-J0) ! How to be careful here that the domains are right?
+      forces%vstk0(i,j) = IOB%vstk0(i-I0,j-J0)
+    enddo ; enddo
+    call pass_vector(forces%ustk0,forces%vstk0, G%domain )
+    do istk = 1,IOB%num_stk_bands
+      do j=js,je; do i=is,ie
+        forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
+        forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
+      enddo; enddo
+      call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
+    enddo
+  endif 
+
   ! sea ice related dynamic fields
   if (associated(IOB%ice_rigidity)) then
     call pass_var(rigidity_at_h, G%Domain, halo=1)

src/core/MOM_forcing_type.F90

@@ -230,6 +230,15 @@ module MOM_forcing_type
                                 !! ice needs to be accumulated, and the rigidity explicitly
                                 !! reset to zero at the driver level when appropriate.
 
+  real, pointer, dimension(:,:) :: &
+       ustk0 => NULL(), &
+       vstk0 => NULL()
+  real, pointer, dimension(:) :: &
+       stk_wavenumbers => NULL()
+  real, pointer, dimension(:,:,:) :: &
+       ustkb => NULL(), &
+       vstkb => NULL()
+
   logical :: initialized = .false. !< This indicates whether the appropriate arrays have been initialized.
 end type mech_forcing
 
@@ -2875,7 +2884,7 @@ subroutine allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, ic
 end subroutine allocate_forcing_type
 
 !> Conditionally allocate fields within the mechanical forcing type
-subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg)
+subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg, waves, num_stk_bands)
   type(ocean_grid_type), intent(in) :: G       !< Ocean grid structure
   type(mech_forcing), intent(inout) :: forces  !< Forcing fields structure
 
@@ -2884,6 +2893,8 @@ subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg
   logical, optional,     intent(in) :: shelf   !< If present and true, allocate forces for ice-shelf
   logical, optional,     intent(in) :: press   !< If present and true, allocate p_surf and related fields
   logical, optional,     intent(in) :: iceberg !< If present and true, allocate forces for icebergs
+  logical, optional,     intent(in) :: waves   !< If present and true, allocate wave fields
+  integer, optional,     intent(in) :: num_stk_bands !< Number of Stokes bands to allocate
 
   ! Local variables
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
@@ -2910,6 +2921,18 @@ subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg
   call myAlloc(forces%area_berg,isd,ied,jsd,jed, iceberg)
   call myAlloc(forces%mass_berg,isd,ied,jsd,jed, iceberg)
 
+  !These fields should only be allocated when waves are being passed through the coupler
+  call myAlloc(forces%ustk0,isd,ied,jsd,jed, waves)
+  call myAlloc(forces%vstk0,isd,ied,jsd,jed, waves)
+  if (present(waves)) then; if (waves) then; if (.not.associated(forces%ustkb)) then
+    if (.not.(present(num_stk_bands))) call MOM_error(FATAL,"Requested to initialize with waves, but no waves are present.")
+    allocate(forces%stk_wavenumbers(num_stk_bands)) ; forces%stk_wavenumbers (:) = 0.0
+    allocate(forces%ustkb(isd:ied,jsd:jed,num_stk_bands)) ; forces%ustkb(isd:ied,jsd:jed,:) = 0.0
+  endif; endif; endif
+  if (present(waves)) then; if (waves) then; if (.not.associated(forces%vstkb)) then
+    allocate(forces%vstkb(isd:ied,jsd:jed,num_stk_bands)) ; forces%vstkb(isd:ied,jsd:jed,:) = 0.0
+  endif; endif; endif
+
 end subroutine allocate_mech_forcing
 
 !> Allocates and zeroes-out array.