introduce a new wave coupling method: VR12-MA (not fully implemented…

… yet)
NOAA-GFDL · Jun 14, 2021 · 0d59156 · 0d59156
1 parent 3812cf6
commit 0d59156
Show file tree

Hide file tree

Showing 3 changed files with 51 additions and 28 deletions.
diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90
@@ -696,17 +696,21 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   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
+    if (ocean_state%wave_method == "VR12-MA") then
+      continue
+    else
+      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
   endif
 
   ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state
@@ -723,9 +727,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes"                , "will provide")
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes"                , "will provide")
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_hstokes"                , "will provide")
-    !call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_melth"                , "will provide")
-    !call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_meltw"                , "will provide")
-    !call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_fswpen"                 , "will provide")
   else
     !call fld_list_add(fldsToOcn_num, fldsToOcn, "mass_of_overlying_sea_ice" , "will provide")
     !call fld_list_add(fldsFrOcn_num, fldsFrOcn, "sea_lev"                   , "will provide")
@@ -753,15 +754,19 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   !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")
+    if (ocean_state%wave_method == "VR12-MA") then
+      continue
+    else
+      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
-    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 -------------

diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
@@ -145,6 +145,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,public :: use_waves !< If true use wave coupling.
+  character(len=40),public :: wave_method !< Wave coupling method.
 
   logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
                               !! ocean dynamics and forcing fluxes.
@@ -387,10 +388,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, &
        "If true, enables surface wave modules.", default=.false.)
   if (OS%use_waves) then
+    call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.)
     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)
+    if (OS%wave_method == "VR12-MA") then
+      continue
+    else
+      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)
+    endif
   else
     call MOM_wave_interface_init_lite(param_file)
   endif
@@ -575,7 +581,9 @@ 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, OS%forces)
+    if (OS%wave_method /= "VR12-MA") then
+      call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
+    endif
   endif
 
   if (OS%nstep==0) then

diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
@@ -141,6 +141,7 @@ module MOM_wave_interface
                           !!   1 - Surface Stokes Drift Bands
                           !!   2 - DHH85
                           !!   3 - LF17
+                          !!   4 - VR12MA
                           !! -99 - No waves computed, but empirical Langmuir number used.
                           !! \todo Module variable! Move into a control structure.
 
@@ -178,7 +179,7 @@ module MOM_wave_interface
 !! into a control structure.
 ! Switches needed in import_stokes_drift
 integer, parameter :: TESTPROF = 0, SURFBANDS = 1, &
-                      DHH85 = 2, LF17 = 3, NULL_WaveMethod=-99, &
+                      DHH85 = 2, LF17 = 3, VR12MA = 4, NULL_WaveMethod=-99, &
                       DATAOVR = 1, COUPLER = 2, INPUT = 3
 
 ! Options For Test Prof
@@ -208,6 +209,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag )
   character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS"
   character*(5), parameter  :: DHH85_STRING     = "DHH85"
   character*(4), parameter  :: LF17_STRING      = "LF17"
+  character*(7), parameter  :: VR12MA_STRING    = "VR12-MA"
   character*(12), parameter :: DATAOVR_STRING   = "DATAOVERRIDE"
   character*(7), parameter  :: COUPLER_STRING   = "COUPLER"
   character*(5), parameter  :: INPUT_STRING     = "INPUT"
@@ -268,7 +270,11 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag )
        " DHH85         - Uses Donelan et al. 1985 empirical \n"//     &
        "                 wave spectrum with prescribed values. \n"//  &
        " LF17          - Infers Stokes drift profile from wind \n"//  &
-       "                 speed following Li and Fox-Kemper 2017.\n",  &
+       "                 speed following Li and Fox-Kemper 2017.\n"// &
+       " VR12-MA       - Applies an enhancement factor to the KPP\n"//&
+       "                 turbulent velocity scale received from \n"// &
+       "                 WW3 and is based on the surface layer \n"//  &
+       "                 and projected Langmuir number. (Li 2016)\n", &
        units='', default=NULL_STRING)
   select case (TRIM(TMPSTRING1))
   case (NULL_STRING)! No Waves
@@ -363,6 +369,8 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag )
           default=.false.)
   case (LF17_STRING)!Li and Fox-Kemper 17 wind-sea Langmuir number
     WaveMethod = LF17
+  case (VR12MA_STRING)!Li and Fox-Kemper 16
+    WaveMethod = VR12MA
   case default
     call MOM_error(FATAL,'Check WAVE_METHOD.')
   end select
@@ -732,6 +740,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar)
       enddo
       DHH85_is_set = .true.
     endif
+  elseif (WaveMethod==VR12MA) then
+    return ! todo
   else! Keep this else, fallback to 0 Stokes drift
     do kk= 1,GV%ke
       do jj = G%jsd,G%jed