Skip to content

Commit

Permalink
Merge pull request #194 from NCAR/merge_ww3_coupling_vr12ma
Browse files Browse the repository at this point in the history 
Introduce the legacy CESM WW3 coupling option (VR12-MA)
  • Loading branch information
@gustavo-marques
gustavo-marques authored Aug 27, 2021
2 parents 3e1b0cd + 88b3642 commit 380615d
Show file tree
Hide file tree
Showing 9 changed files with 291 additions and 177 deletions.
60 changes: 29 additions & 31 deletions config_src/drivers/nuopc_cap/mom_cap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -704,19 +705,24 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
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
Expand All @@ -730,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
Expand All @@ -767,15 +761,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 (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 -------------
Expand Down
10 changes: 10 additions & 0 deletions config_src/drivers/nuopc_cap/mom_cap_methods.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,6 +271,16 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
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
!----
Expand Down
17 changes: 13 additions & 4 deletions config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,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.
Expand Down Expand Up @@ -380,7 +381,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
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, &
Expand All @@ -392,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)
Expand Down Expand Up @@ -580,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
Expand Down Expand Up @@ -1010,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
Expand All @@ -1030,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

Expand Down
36 changes: 25 additions & 11 deletions config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,7 @@ 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(:,:) :: 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]
Expand Down Expand Up @@ -545,17 +546,25 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)

! sea ice fraction [nondim]
if (associated(IOB%ice_fraction) .and. associated(fluxes%ice_fraction)) &
fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0)
! 10-m wind speed squared [m2/s2]
if (associated(IOB%u10_sqr) .and. associated(fluxes%u10_sqr)) &
fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0)

enddo ; enddo

if (CS%use_CFC) then
do j=js,je ; do i=is,ie
! sea ice fraction [nondim]
if (associated(IOB%ice_fraction)) &
fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0)
! 10-m wind speed squared [m2/s2]
if (associated(IOB%u10_sqr)) &
fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0)
! wave to ocean coupling
if ( associated(IOB%lamult)) then
do j=js,je; do i=is,ie
if (IOB%ice_fraction(i-i0,j-j0) <= 0.05 ) then
fluxes%lamult(i,j) = IOB%lamult(i-i0,j-j0)
else
fluxes%lamult(i,j) = 1.0
endif
enddo ; enddo
call pass_var(fluxes%lamult, G%domain, halo=1 )
endif

! applied surface pressure from atmosphere and cryosphere
Expand Down Expand Up @@ -717,8 +726,11 @@ 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) ) &
if ( associated(IOB%lamult) ) then
call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=0)
elseif ( associated(IOB%ustkb) ) then
call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands)
endif

! applied surface pressure from atmosphere and cryosphere
if (CS%use_limited_P_SSH) then
Expand Down Expand Up @@ -1070,7 +1082,7 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
end subroutine forcing_save_restart

!> Initialize the surface forcing, including setting parameters and allocating permanent memory.
subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp)
subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp, use_waves)
type(time_type), intent(in) :: Time !< The current model time
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
Expand All @@ -1083,6 +1095,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
!! restoring will be applied in this model.
logical, optional, intent(in) :: restore_temp !< If present and true surface temperature
!! restoring will be applied in this model.
logical, optional, intent(in) :: use_waves !< If present and true, use waves and activate
!! the corresponding wave forcing diagnostics

! Local variables
real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1].
Expand Down Expand Up @@ -1354,7 +1368,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
"as seen by MOM6.", default=.false.)

call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, &
use_berg_fluxes=iceberg_flux_diags, use_cfcs=CS%use_CFC)
use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves, use_cfcs=CS%use_CFC)

call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, &
"If true, allows flux adjustments to specified via the "//&
Expand Down
Loading

0 comments on commit 380615d

Please sign in to comment.