Merge remote-tracking branch 'NCAR/surfbands_refactor' into BGR_Stoke…

…s_DDT_Fixes

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -714,16 +714,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
       Ice_ocean_boundary%lamult          = 0.0
     else if (wave_method == "SURFACE_BANDS") then
       call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands)
-      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
+      allocate(Ice_ocean_boundary%ustkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), source=0.0)
+      allocate(Ice_ocean_boundary%vstkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), source=0.0)
+      allocate(Ice_ocean_boundary%stk_wavenumbers(Ice_ocean_boundary%num_stk_bands), source=0.0)
       call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.)
-      Ice_ocean_boundary%ustkb           = 0.0
-      Ice_ocean_boundary%vstkb           = 0.0
     else
       call MOM_error(FATAL, "Unsupported WAVE_METHOD encountered in NUOPC cap.")
     endif

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -373,6 +373,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, &
                  default=.false., do_not_log=.true.)
+  call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, &
+       "If true, enables surface wave modules.", default=.false.)
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
@@ -393,12 +395,11 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
       call allocate_forcing_type(OS%grid, OS%fluxes, shelf=.true.)
   endif
 
-  call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true.)
-  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.)
   endif
+  call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true., lamult=(trim(OS%wave_method)=="EFACTOR"))
+
   ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because
   ! it also initializes statistical waves.
   call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag, OS%restart_CSp)

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -193,8 +193,6 @@ module MOM_surface_forcing_nuopc
                                                               !! for divergence damping, as determined
                                                               !! outside of the ocean model in [m3/s]
   real, pointer, dimension(:,:)   :: lamult          => NULL() !< Langmuir enhancement factor [nondim]
-  real, pointer, dimension(:,:)   :: ustk0           => NULL() !< Surface Stokes drift, zonal [m/s]
-  real, pointer, dimension(:,:)   :: vstk0           => NULL() !< Surface Stokes drift, meridional [m/s]
   real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m]
   real, pointer, dimension(:,:,:) :: ustkb           => NULL() !< Stokes Drift spectrum, zonal [m/s]
                                                                !! Horizontal  - u points
@@ -893,17 +891,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
   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 )
+      call pass_var(forces%ustkb(:,:,istk), G%domain )
+      call pass_var(forces%vstkb(:,:,istk), G%domain )
     enddo
   endif
 

src/core/MOM_forcing_type.F90

@@ -256,9 +256,6 @@ module MOM_forcing_type
   logical :: accumulate_rigidity = .false. !< If true, the rigidity due to various types of
                                 !! ice needs to be accumulated, and the rigidity explicitly
                                 !! reset to zero at the driver level when appropriate.
-  real, pointer, dimension(:,:) :: &
-       ustk0 => NULL(), &       !< Surface Stokes drift, zonal [m/s]
-       vstk0 => NULL()          !< Surface Stokes drift, meridional [m/s]
   real, pointer, dimension(:) :: &
        stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m]
   real, pointer, dimension(:,:,:) :: &
@@ -2953,7 +2950,7 @@ end subroutine forcing_diagnostics
 
 !> Conditionally allocate fields within the forcing type
 subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
-                                     shelf, iceberg, salt, fix_accum_bug, cfc, waves)
+                                     shelf, iceberg, salt, fix_accum_bug, cfc, waves, lamult)
   type(ocean_grid_type), intent(in) :: G       !< Ocean grid structure
   type(forcing),      intent(inout) :: fluxes  !< A structure containing thermodynamic forcing fields
   logical, optional,     intent(in) :: water   !< If present and true, allocate water fluxes
@@ -2967,6 +2964,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
                                                !! accumulation of ustar_gustless
   logical, optional,     intent(in) :: cfc     !< If present and true, allocate cfc fluxes
   logical, optional,     intent(in) :: waves   !< If present and true, allocate wave fields
+  logical, optional,     intent(in) :: lamult  !< If present and true, allocate langmuir enhancement factor
 
   ! Local variables
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
@@ -3030,7 +3028,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, &
 
   !These fields should only on allocated when wave coupling is activated.
   call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, waves)
-  call myAlloc(fluxes%lamult,isd,ied,jsd,jed, waves)
+  call myAlloc(fluxes%lamult,isd,ied,jsd,jed, lamult)
 
   if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug
 end subroutine allocate_forcing_by_group
@@ -3125,29 +3123,20 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, &
   call myAlloc(forces%mass_berg,isd,ied,jsd,jed, iceberg)
 
   !These fields should only be allocated when waves
-  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. present(num_stk_bands)) then
       call MOM_error(FATAL,"Requested to &
       initialize with waves, but no waves are present.")
     endif
     if (num_stk_bands > 0) then
       if (.not.associated(forces%ustkb)) then
-        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
+        allocate(forces%stk_wavenumbers(num_stk_bands), source=0.0)
+        allocate(forces%ustkb(isd:ied,jsd:jed,num_stk_bands), source=0.0)
+        allocate(forces%vstkb(isd:ied,jsd:jed,num_stk_bands), source=0.0)
       endif
     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_by_group
 
 

src/user/MOM_wave_interface.F90

@@ -465,12 +465,9 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag, restar
     !allocate(CS%Us_x_prev(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
     !allocate(CS%Us_y_prev(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
     allocate(CS%ddt_Us_x(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
-    CS%ddt_Us_x(:,:,:) = 0.0
     allocate(CS%ddt_Us_y(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
-    CS%ddt_Us_y(:,:,:) = 0.0
     allocate(CS%Us_x_from_ddt(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
     allocate(CS%Us_y_from_ddt(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
-
   endif
   ! b. Surface Values
   allocate(CS%US0_x(G%isdB:G%iedB,G%jsd:G%jed), source=0.0)
@@ -1542,7 +1539,7 @@ subroutine CoriolisStokes(G, GV, dt, h, u, v, WAVES, US)
   integer :: i,j,k
 
   do k = 1, GV%ke
-do j = G%jsc, G%jec
+    do j = G%jsc, G%jec
       do I = G%iscB, G%iecB
         DVel = 0.25*(WAVES%us_y(i,j+1,k)+WAVES%us_y(i-1,j+1,k))*G%CoriolisBu(i,j+1) + &
                0.25*(WAVES%us_y(i,j,k)+WAVES%us_y(i-1,j,k))*G%CoriolisBu(i,j)