Merge pull request mom-ocean#1491 from gustavo-marques/dev-ncar-main-candidate-2021-08-30

Dev ncar main candidate 2021 08 31

Author: marshallward

config_src/drivers/mct_cap/mom_surface_forcing_mct.F90

@@ -1245,7 +1245,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  "If true, use a 2-dimensional gustiness supplied from "//&
                  "an input file", default=.false.)
   call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
-               "The background gustiness in the winds.", units="Pa", default=0.0)
+                 "The background gustiness in the winds.", units="Pa", default=0.0, &
+                 scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z)
   if (CS%read_gust_2d) then
     call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
                  "The file in which the wind gustiness is found in "//&

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -31,7 +31,7 @@ module MOM_cap_mod
 use MOM_grid,                 only: ocean_grid_type, get_global_grid_size
 use MOM_ocean_model_nuopc,    only: ice_ocean_boundary_type
 use MOM_ocean_model_nuopc,    only: ocean_model_restart, ocean_public_type, ocean_state_type
-use MOM_ocean_model_nuopc,    only: ocean_model_init_sfc
+use MOM_ocean_model_nuopc,    only: ocean_model_init_sfc, ocean_model_flux_init
 use MOM_ocean_model_nuopc,    only: ocean_model_init, update_ocean_model, ocean_model_end
 use MOM_ocean_model_nuopc,    only: get_ocean_grid, get_eps_omesh, query_ocean_state
 use MOM_cap_time,             only: AlarmInit
@@ -422,6 +422,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   logical                                :: isPresent, isPresentDiro, isPresentLogfile, isSet
   logical                                :: existflag
   logical                                :: use_waves  ! If true, the wave modules are active.
+  character(len=40)                      :: wave_method ! Wave coupling method.
   integer                                :: userRc
   integer                                :: localPet
   integer                                :: localPeCount
@@ -650,6 +651,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   ocean_public%is_ocean_pe = .true.
   call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles))
 
+  ! GMM, this call is not needed for NCAR. Check with EMC.
+  ! If this can be deleted, perhaps we should also delete ocean_model_flux_init
+  call ocean_model_flux_init(ocean_state)
+
   call ocean_model_init_sfc(ocean_state, ocean_public)
 
   call mpp_get_compute_domain(ocean_public%domain, isc, iec, jsc, jec)
@@ -669,6 +674,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
              Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),&
              Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec),     &
              Ice_ocean_boundary% mi (isc:iec,jsc:jec),              &
+             Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec),    &
+             Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec),         &
              Ice_ocean_boundary% p (isc:iec,jsc:jec),               &
              Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec),    &
              Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec),    &
@@ -690,25 +697,32 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   Ice_ocean_boundary%seaice_melt     = 0.0
   Ice_ocean_boundary%seaice_melt_heat= 0.0
   Ice_ocean_boundary%mi              = 0.0
+  Ice_ocean_boundary%ice_fraction    = 0.0
+  Ice_ocean_boundary%u10_sqr         = 0.0
   Ice_ocean_boundary%p               = 0.0
   Ice_ocean_boundary%lrunoff_hflx    = 0.0
   Ice_ocean_boundary%frunoff_hflx    = 0.0
   Ice_ocean_boundary%lrunoff         = 0.0
   Ice_ocean_boundary%frunoff         = 0.0
 
-  call query_ocean_state(ocean_state, use_waves=use_waves)
+  call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method)
   if (use_waves) 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
-    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
+    if (wave_method == "EFACTOR") then
+      allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) )
+      Ice_ocean_boundary%lamult          = 0.0
+    else
+      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
+      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
+    endif
   endif
   ! Consider adding this:
   ! if (.not.use_waves) Ice_ocean_boundary%num_stk_bands = 0
@@ -722,18 +736,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
    call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide")
   end if
 
-  if (cesm_coupled) then
-    !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult"                 , "will provide")
-    !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")
-  endif
 
   !--------- import fields -------------
   call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_salt_rate"             , "will provide") ! from ice
@@ -751,21 +753,27 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface"   , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"                  , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"                  , "will provide") !-> ice runoff
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac"                   , "will provide") !-> ice fraction
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"                  , "will provide") !-> wind^2 at 10m
   call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn"        , "will provide")
   !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 (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 (wave_method == "EFACTOR") then
+      call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult"                 , "will provide")
+    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 -------------
@@ -1103,7 +1111,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
            k = k + 1 ! Increment position within gindex
            if (mask(k) /= 0) then
               mesh_areas(k) = dataPtr_mesh_areas(k)
-              model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
+              model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
               mod2med_areacor(k) = model_areas(k) / mesh_areas(k)
               med2mod_areacor(k) = mesh_areas(k) / model_areas(k)
            end if

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -250,12 +250,37 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   !----
   ! Note - preset values to 0, if field does not exist in importState, then will simply return
   ! and preset value will be used
-
   ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8
   call state_getimport(importState, 'mass_of_overlying_ice',  &
        isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+  !----
+  ! sea-ice fraction
+  !----
+  ice_ocean_boundary%ice_fraction(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'Si_ifrac',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%ice_fraction, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+  !----
+  ! 10m wind squared
+  !----
+  ice_ocean_boundary%u10_sqr(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'So_duu10n',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%u10_sqr, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+  !----
+  ! Langmuir enhancement factor
+  !----
+  if ( associated(ice_ocean_boundary%lamult) ) then
+   ice_ocean_boundary%lamult (:,:) = 0._ESMF_KIND_R8
+   call state_getimport(importState, 'Sw_lamult',  &
+        isc, iec, jsc, jec, ice_ocean_boundary%lamult, rc=rc)
+   if (ChkErr(rc,__LINE__,u_FILE_u)) return
+   endif
+
   !----
   ! Partitioned Stokes Drift Components
   !----

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -41,7 +41,6 @@ module MOM_ocean_model_nuopc
 use MOM_time_manager,        only : operator(/=), operator(<=), operator(>=)
 use MOM_time_manager,        only : operator(<), real_to_time_type, time_type_to_real
 use time_interp_external_mod,only : time_interp_external_init
-use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init
 use MOM_tracer_flow_control, only : call_tracer_flux_init
 use MOM_unit_scaling,        only : unit_scale_type
 use MOM_variables,           only : surface
@@ -147,6 +146,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) :: wave_method !< Wave coupling method.
 
   logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
                               !! ocean dynamics and forcing fluxes.
@@ -242,14 +242,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
   character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read
+
   ! Local variables
   real :: Rho0        ! The Boussinesq ocean density, in kg m-3.
   real :: G_Earth     ! The gravitational acceleration in m s-2.
   real :: HFrz        !< If HFrz > 0 (m), melt potential will be computed.
                       !! The actual depth over which melt potential is computed will
                       !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
                       !! If HFrz <= 0 (default), melt potential will not be computed.
-  logical :: use_melt_pot!< If true, allocate melt_potential array
+  logical :: use_melt_pot !< If true, allocate melt_potential array
+  logical :: use_CFC  !< If true, allocated arrays for surface CFCs.
+
 
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
@@ -368,13 +371,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
     use_melt_pot=.false.
   endif
 
+  call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, &
+                 default=.false., do_not_log=.true.)
+
   !   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.
   call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
-                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
+                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn, &
+                              use_meltpot=use_melt_pot, use_cfcs=use_CFC)
 
   call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, &
-                            OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)
+                            OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves)
 
   if (OS%use_ice_shelf)  then
     call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, &
@@ -386,8 +393,12 @@ 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
   ! 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)
@@ -413,6 +424,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   endif
 
+  call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
+
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
 
@@ -572,7 +585,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 /= "EFACTOR") 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
@@ -1002,15 +1017,16 @@ end subroutine ocean_model_flux_init
 
 !> This interface allows certain properties that are stored in the ocean_state_type to be
 !! obtained.
-subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale)
+subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale, wave_method)
   type(ocean_state_type),       intent(in)  :: OS      !< The structure with the complete ocean state
   logical,            optional, intent(out) :: use_waves !< Indicates whether surface waves are in use
   integer,            optional, intent(out) :: NumWaveBands !< If present, this gives the number of
                                                        !! wavenumber partitions in the wave discretization
   real, dimension(:), optional, intent(out) :: Wavenumbers !< If present, this gives the characteristic
                                                        !! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1]
-    logical,          optional, intent(in)  :: unscale !< If present and true, undo any dimensional
+  logical,            optional, intent(in)  :: unscale !< If present and true, undo any dimensional
                                                        !! rescaling and return dimensional values in MKS units
+  character(len=40),  optional, intent(out) :: wave_method !< Wave coupling method.
 
   logical :: undo_scaling
   undo_scaling = .false. ; if (present(unscale)) undo_scaling = unscale
@@ -1022,6 +1038,7 @@ subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale)
   elseif (present(Wavenumbers)) then
     call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers)
   endif
+  if (present(wave_method)) wave_method = OS%wave_method
 
 end subroutine query_ocean_state