diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 0090468a05..67d650aded 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -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 @@ -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 @@ -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 @@ -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 ------------- diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index b7e3b13c1e..bb9743bb84 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -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 !---- diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index e193c9ee38..39e1046c1c 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -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. @@ -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, & @@ -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) @@ -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 @@ -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 @@ -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 diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index c83e32711c..823d79c951 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -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] @@ -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 @@ -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 @@ -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 @@ -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]. @@ -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 "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7d8375a8af..7022c8412c 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -190,6 +190,9 @@ module MOM_forcing_type ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] + real, pointer, dimension(:,:) :: & + lamult => NULL() !< Langmuir enhancement factor [nondim] + ! passive tracer surface fluxes type(coupler_2d_bc_type) :: tr_fluxes !< This structure contains arrays of !! of named fields used for passive tracer fluxes. @@ -253,7 +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] @@ -369,6 +371,9 @@ module MOM_forcing_type ! Iceberg + Ice shelf diagnostic handles integer :: id_ustar_ice_cover = -1 integer :: id_frac_ice_cover = -1 + + ! wave forcing diagnostics handles. + integer :: id_lamult = -1 !>@} integer :: id_clock_forcing = -1 !< CPU clock id @@ -1256,13 +1261,14 @@ end subroutine forcing_SinglePointPrint !> Register members of the forcing type for diagnostics -subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_cfcs) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_waves, use_cfcs) type(time_type), intent(in) :: Time !< time type type(diag_ctrl), intent(inout) :: diag !< diagnostic control type type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< True if T/S are in use type(forcing_diags), intent(inout) :: handles !< handles for diagnostics logical, optional, intent(in) :: use_berg_fluxes !< If true, allow iceberg flux diagnostics + logical, optional, intent(in) :: use_waves !< If true, allow wave forcing diagnostics logical, optional, intent(in) :: use_cfcs !< If true, allow cfc related diagnostics ! Clock for forcing diagnostics @@ -1941,6 +1947,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_total_saltFluxAdded = register_scalar_field('ocean_model', 'total_salt_Flux_Added', & Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg s-1') + !=============================================================== + ! wave forcing diagnostics + if (present(use_waves)) then + if (use_waves) then + handles%id_lamult = register_diag_field('ocean_model', 'lamult', & + diag%axesT1, Time, long_name='Langmuir enhancement factor received from WW3', units="nondim") + endif + endif end subroutine register_forcing_type_diags @@ -2917,6 +2931,10 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_ustar_ice_cover > 0) .and. associated(fluxes%ustar_shelf)) & call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) + ! wave forcing =============================================================== + if (handles%id_lamult > 0) & + call post_data(handles%id_lamult, fluxes%lamult, diag) + ! endif ! query_averaging_enabled call disable_averaging(diag) @@ -2931,7 +2949,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) + shelf, iceberg, salt, fix_accum_bug, cfc, waves) 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 @@ -2944,6 +2962,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: fix_accum_bug !< If present and true, avoid using a bug in !! 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 ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -3005,6 +3024,10 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) + !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) + if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group @@ -3100,14 +3123,21 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & !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.associated(forces%ustkb)) then - if (.not.(present(num_stk_bands))) call MOM_error(FATAL,"Requested to & + 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.") - 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 + 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 + 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)) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 083f5ed000..a26f251a2b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -15,7 +15,7 @@ module MOM_CVMix_KPP use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number +use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number, get_wave_method use MOM_domains, only : pass_var use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE @@ -198,6 +198,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) # include "version_variable.h" character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module character(len=20) :: string !< local temporary string + character(len=20) :: langmuir_mixing_opt = 'NONE' !< langmuir mixing opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_entrainment_opt = 'NONE' !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: wave_method logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function) !! False => compute G'(1) as in LMD94 @@ -415,10 +418,17 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) '\t RW16 = Function of Langmuir number based on RW16', & default='CONSTANT') select case ( trim(string)) - case ("CONSTANT") ; CS%LT_K_METHOD = LT_K_MODE_CONSTANT - case ("VR12") ; CS%LT_K_METHOD = LT_K_MODE_VR12 - case ("RW16") ; CS%LT_K_METHOD = LT_K_MODE_RW16 - case default ; call MOM_error(FATAL,"KPP_init: "//& + case ("CONSTANT") + CS%LT_K_METHOD = LT_K_MODE_CONSTANT + langmuir_mixing_opt = 'LWF16' + case ("VR12") + CS%LT_K_METHOD = LT_K_MODE_VR12 + langmuir_mixing_opt = 'LWF16' + case ("RW16") + CS%LT_K_METHOD = LT_K_MODE_RW16 + langmuir_mixing_opt = 'RWHGK16' + case default + call MOM_error(FATAL,"KPP_init: "//& "Unrecognized KPP_LT_K_METHOD option: "//trim(string)) end select if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then @@ -441,12 +451,21 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) '\t LF17 = Function of Langmuir number based on LF17', & default='CONSTANT') select case ( trim(string)) - case ("CONSTANT") ; CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT - case ("VR12") ; CS%LT_VT2_METHOD = LT_VT2_MODE_VR12 - case ("RW16") ; CS%LT_VT2_METHOD = LT_VT2_MODE_RW16 - case ("LF17") ; CS%LT_VT2_METHOD = LT_VT2_MODE_LF17 - case default ; call MOM_error(FATAL,"KPP_init: "//& - "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string)) + case ("CONSTANT") + CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT + langmuir_entrainment_opt = 'LWF16' + case ("VR12") + CS%LT_VT2_METHOD = LT_VT2_MODE_VR12 + langmuir_entrainment_opt = 'LWF16' + case ("RW16") + CS%LT_VT2_METHOD = LT_VT2_MODE_RW16 + langmuir_entrainment_opt = 'RWHGK16' + case ("LF17") + CS%LT_VT2_METHOD = LT_VT2_MODE_LF17 + langmuir_entrainment_opt = 'LF17' + case default + call MOM_error(FATAL,"KPP_init: "//& + "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string)) end select if (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then call get_param(paramFile, mdl, "KPP_VT2_ENH_FAC",CS%KPP_VT2_ENH_FAC , & @@ -456,6 +475,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) endif call closeParameterBlock(paramFile) + call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) call CVMix_init_kpp( Ri_crit=CS%Ri_crit, & @@ -471,6 +491,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) lenhanced_diff=CS%enhance_diffusion,& lnonzero_surf_nonlocal=Cs_is_one ,& lnoDGat1=lnoDGat1 ,& + langmuir_mixing_str=langmuir_mixing_opt,& + langmuir_entrainment_str=langmuir_entrainment_opt,& CVMix_kpp_params_user=CS%KPP_params ) ! Register diagnostics @@ -600,7 +622,7 @@ end function KPP_init !> KPP vertical diffusivity/viscosity and non-local tracer transport subroutine KPP_calculate(CS, G, GV, US, h, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& - nonLocalTransScalar, waves) + nonLocalTransScalar, waves, lamult) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -621,7 +643,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & !! [Z2 T-1 ~> m2 s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport [m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local trans. [m s-1] - type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement multiplier ! Local variables integer :: i, j, k ! Loop indices @@ -661,7 +684,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & !$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, & !$OMP sigmaRatio) & !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, & - !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves) + !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves, lamult) ! loop over horizontal points on processor do j = G%jsc, G%jec do i = G%isc, G%iec @@ -714,6 +737,48 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & Kviscosity(:) = US%Z2_T_to_m2_s * Kv(i,j,:) endif + IF (CS%LT_K_ENHANCEMENT) then + if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then + LangEnhK = CS%KPP_K_ENH_FAC + elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then + if (present(lamult)) then + LangEnhK = lamult(i,j) + else + LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & + (5.4*CS%La_SL(i,j))**(-4)) + endif + elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then + !This maximum value is proposed in Reichl et al., 2016 JPO formula + LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j)) + else + !This shouldn't be reached. + !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT") + LangEnhK = 1.0 + endif + + ! diffusivities don't need to be enhanced below anymore since LangEnhK is applied within CVMix. + ! todo: need to double check if the distinction between the two different options of LT_K_SHAPE may need to be + ! treated specially. + do k=1,GV%ke + if (CS%LT_K_SHAPE== LT_K_CONSTANT) then + if (CS%id_EnhK > 0) CS%EnhK(i,j,:) = LangEnhK + !Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK + !Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK + !Kviscosity(k) = Kviscosity(k) * LangEnhK + elseif (CS%LT_K_SHAPE == LT_K_SCALED) then + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) + SigmaRatio = sigma * (1. - sigma)**2 / 0.148148037 + if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio) + !Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) + !Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) + !Kviscosity(k) = Kviscosity(k) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) + endif + enddo + endif + call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] @@ -730,6 +795,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] GV%ke, & ! (in) Number of levels to compute coeffs for GV%ke, & ! (in) Number of levels in array shape + Langmuir_EFactor=LangEnhK,& ! Langmuir enhancement multiplier CVMix_kpp_params_user=CS%KPP_params ) ! safety check, Kviscosity and Kdiffusivity must be >= 0 @@ -742,41 +808,6 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & endif enddo - IF (CS%LT_K_ENHANCEMENT) then - if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then - LangEnhK = CS%KPP_K_ENH_FAC - elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then - ! Added minimum value for La_SL, so removed maximum value for LangEnhK. - LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & - (5.4*CS%La_SL(i,j))**(-4)) - elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then - !This maximum value is proposed in Reichl et al., 2016 JPO formula - LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j)) - else - !This shouldn't be reached. - !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT") - LangEnhK = 1.0 - endif - do k=1,GV%ke - if (CS%LT_K_SHAPE== LT_K_CONSTANT) then - if (CS%id_EnhK > 0) CS%EnhK(i,j,:) = LangEnhK - Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK - Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK - Kviscosity(k) = Kviscosity(k) * LangEnhK - elseif (CS%LT_K_SHAPE == LT_K_SCALED) then - sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) - SigmaRatio = sigma * (1. - sigma)**2 / 0.148148037 - if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio) - Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) - Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) - Kviscosity(k) = Kviscosity(k) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) - endif - enddo - endif - ! Over-write CVMix NLT shape function with one of the following choices. ! The CVMix code has yet to update for thse options, so we compute in MOM6. ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with @@ -827,18 +858,6 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & Kviscosity(:) = 0.0 endif - - ! compute unresolved squared velocity for diagnostics - if (CS%id_Vt2 > 0) then -!BGR Now computing VT2 above so can modify for LT -! therefore, don't repeat this operation here -! CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & -! cellHeight(1:GV%ke), & ! Depth of cell center [m] -! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] -! N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] -! CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - endif - ! Copy 1d data into 3d diagnostic arrays !/ grabbing obldepth_0d for next time step. CS%OBLdepthprev(i,j)=CS%OBLdepth(i,j) @@ -886,7 +905,6 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%id_OBLdepth_original > 0) call post_data(CS%id_OBLdepth_original,CS%OBLdepth_original,CS%diag) if (CS%id_sigma > 0) call post_data(CS%id_sigma, CS%sigma, CS%diag) if (CS%id_Ws > 0) call post_data(CS%id_Ws, CS%Ws, CS%diag) - if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) if (CS%id_uStar > 0) call post_data(CS%id_uStar, uStar, CS%diag) if (CS%id_buoyFlux > 0) call post_data(CS%id_buoyFlux, buoyFlux, CS%diag) if (CS%id_Kt_KPP > 0) call post_data(CS%id_Kt_KPP, CS%Kt_KPP, CS%diag) @@ -900,7 +918,7 @@ end subroutine KPP_calculate !> Compute OBL depth -subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFlux, Waves) +subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFlux, Waves, lamult) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -914,8 +932,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Velocity j-component [L T-1 ~> m s-1] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult!< Langmuir enhancement factor ! Local variables integer :: i, j, k, km1 ! Loop indices @@ -953,8 +972,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl integer :: kk, ksfc, ktmp ! For Langmuir Calculations - real :: LangEnhW ! Langmuir enhancement for turbulent velocity scale - real, dimension(GV%ke) :: LangEnhVt2 ! Langmuir enhancement for unresolved shear + real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear real, dimension(GV%ke) :: U_H, V_H real :: MLD_GUESS, LA real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir @@ -987,7 +1005,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, & !$OMP BulkRi_1d, zBottomMinusOffset) & !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & - !$OMP Temp, Salt, waves, tv, GoRho, u, v) + !$OMP Temp, Salt, waves, tv, GoRho, u, v, lamult) do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1146,67 +1164,43 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params ) - !Compute CVMix VT2 - CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & - zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] - N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] - CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - - !Modify CVMix VT2 + ! Determine the enhancement factor for unresolved shear IF (CS%LT_VT2_ENHANCEMENT) then IF (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then - do k=1,GV%ke - LangEnhVT2(k) = CS%KPP_VT2_ENH_FAC - enddo + LangEnhVT2 = CS%KPP_VT2_ENH_FAC elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed. - enhvt2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & - (5.4*CS%La_SL(i,j))**(-4)) - do k=1,GV%ke - LangEnhVT2(k) = enhvt2 - enddo - elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_RW16) then - !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed. - enhvt2 = 1. + 2.3*CS%La_SL(i,j)**(-0.5) - do k=1,GV%ke - LangEnhVT2(k) = enhvt2 - enddo - elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_LF17) then - CS%CS=cvmix_get_kpp_real('c_s',CS%KPP_params) - do k=1,GV%ke - WST = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellHeight(k)))**(1./3.) - LangEnhVT2(k) = sqrt((0.15*WST**3. + 0.17*surfFricVel**3.* & - (1.+0.49*CS%La_SL(i,j)**(-2.))) / & - (0.2*ws_1d(k)**3/(CS%cs*CS%surf_layer_ext*CS%vonKarman**4.))) - enddo + if (present(lamult)) then + LangEnhVT2 = lamult(i,j) + else + LangEnhVT2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & + (5.4*CS%La_SL(i,j))**(-4)) + endif else - !This shouldn't be reached. - !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2") - LangEnhVT2(:) = 1.0 + ! for other methods (e.g., LT_VT2_MODE_RW16, LT_VT2_MODE_LF17), the enhancement factor is + ! computed internally within CVMix using LaSL, bfsfc, and ustar to be passed to CVMix. + LangEnhVT2 = 1.0 endif else - LangEnhVT2(:) = 1.0 + LangEnhVT2 = 1.0 endif - do k=1,GV%ke - CS%Vt2(i,j,k)=CS%Vt2(i,j,k)*LangEnhVT2(k) - if (CS%id_EnhVt2 > 0) CS%EnhVt2(i,j,k)=LangEnhVT2(k) - enddo + surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! Calculate Bulk Richardson number from eq (21) of LMD94 BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] - Vt_sqr_cntr=CS%Vt2(i,j,:), & ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] - N_iface=CS%N(i,j,:)) ! Buoyancy frequency [s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1] + EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] + LaSL = CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] + bfsfc = surfBuoyFlux, & ! surface buoyancy flux [m2 s-3] + uStar = uStar(i,j), & ! surface friction velocity [m s-1] + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters - surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit - ! h to Monin-Obukov (default is false, ie. not used) - call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number iFaceHeight, & ! (in) Height of interfaces [m] @@ -1231,6 +1225,18 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + ! compute unresolved squared velocity for diagnostics + if (CS%id_Vt2 > 0) then + CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & + cellHeight(1:GV%ke), & ! Depth of cell center [m] + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] + EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] + LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] + bfsfc=surfBuoyFlux, & ! surface buoyancy flux [m2 s-3] + uStar=uStar(i,j), & ! surface friction velocity [m s-1] + CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + endif ! recompute wscale for diagnostics, now that we in fact know boundary layer depth !BGR consider if LTEnhancement is wanted for diagnostics @@ -1273,6 +1279,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%id_EnhK > 0) call post_data(CS%id_EnhK, CS%EnhK, CS%diag) if (CS%id_EnhVt2 > 0) call post_data(CS%id_EnhVt2, CS%EnhVt2, CS%diag) if (CS%id_La_SL > 0) call post_data(CS%id_La_SL, CS%La_SL, CS%diag) + if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) ! BLD smoothing: if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 7a75802a84..b44d4c84b4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1190,11 +1190,19 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) ! The KPP scheme calculates boundary layer diffusivities and non-local transport. - call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + if ( associated(fluxes%lamult) ) then + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + else + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + endif if (associated(Hml)) then call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 38aa6b13a5..2f0d95a62d 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -39,6 +39,7 @@ module MOM_wave_interface ! serious consideration of the full 3d wave-averaged Navier-Stokes ! CL2 effects. public Waves_end ! public interface to deallocate and free wave related memory. +public get_wave_method ! public interface to obtain the wave method string ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -173,12 +174,20 @@ module MOM_wave_interface ! Switches needed in import_stokes_drift !>@{ Enumeration values for the wave method -integer, parameter :: TESTPROF = 0, SURFBANDS = 1, DHH85 = 2, LF17 = 3, NULL_WaveMethod = -99 +integer, parameter :: TESTPROF = 0, SURFBANDS = 1, DHH85 = 2, LF17 = 3, EFACTOR = 4, NULL_WaveMethod = -99 !>@} !>@{ Enumeration values for the wave data source integer, parameter :: DATAOVR = 1, COUPLER = 2, INPUT = 3 !>@} +! Strings for the wave method +character*(5), parameter :: NULL_STRING = "EMPTY" !< null wave method string +character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" !< test profile string +character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" !< surface bands string +character*(5), parameter :: DHH85_STRING = "DHH85" !< DHH85 wave method string +character*(4), parameter :: LF17_STRING = "LF17" !< LF17 wave method string +character*(7), parameter :: EFACTOR_STRING = "EFACTOR" !< EFACTOR (based on vr12-ma) wave method string + contains !> Initializes parameters related to MOM_wave_interface @@ -196,11 +205,6 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) ! This include declares and sets the variable "version". # include "version_variable.h" character*(13) :: TMPSTRING1, TMPSTRING2 - character*(5), parameter :: NULL_STRING = "EMPTY" - character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" - character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" - character*(5), parameter :: DHH85_STRING = "DHH85" - character*(4), parameter :: LF17_STRING = "LF17" character*(12), parameter :: DATAOVR_STRING = "DATAOVERRIDE" character*(7), parameter :: COUPLER_STRING = "COUPLER" character*(5), parameter :: INPUT_STRING = "INPUT" @@ -283,7 +287,12 @@ 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"// & + " EFACTOR - Applies an enhancement factor to the KPP\n"//& + " turbulent velocity scale received \n"// & + " directly from WW3 and is based on the \n"// & + " surface layer and projected Langmuir \n"// & + " number (Li 2016)\n", & units='', default=NULL_STRING) select case (TRIM(TMPSTRING1)) case (NULL_STRING)! No Waves @@ -382,6 +391,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 CS%WaveMethod = LF17 + case (EFACTOR_STRING)!Li and Fox-Kemper 16 + CS%WaveMethod = EFACTOR case default call MOM_error(FATAL,'Check WAVE_METHOD.') end select @@ -760,6 +771,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo CS%DHH85_is_set = .true. endif + elseif (CS%WaveMethod==EFACTOR) then + return ! pass else! Keep this else, fallback to 0 Stokes drift do kk= 1,GV%ke do jj = G%jsd,G%jed @@ -1046,6 +1059,31 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, & end subroutine get_Langmuir_Number +!> function to return the wave method string set in the param file +function get_wave_method(CS) + character(:), allocatable :: get_wave_method + type(wave_parameters_CS), pointer :: CS !< Control structure + + if (associated(CS)) then + select case(CS%WaveMethod) + case (Null_WaveMethod) + get_wave_method = NULL_STRING + case (TESTPROF) + get_wave_method = TESTPROF_STRING + case (SURFBANDS) + get_wave_method = SURFBANDS_STRING + case (DHH85) + get_wave_method = DHH85_STRING + case (LF17) + get_wave_method = LF17_STRING + case (EFACTOR) + get_wave_method = EFACTOR_STRING + end select + else + get_wave_method = NULL_STRING + endif +end function get_wave_method + !> Get SL averaged Stokes drift from Li/FK 17 method !! !! Original description: diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index b3c8f45843..136e5f9eee 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -49,7 +49,7 @@ function register_tidal_bay_OBC(param_file, CS, US, OBC_Reg) call get_param(param_file, mdl, "TIDAL_BAY_FLOW", CS%tide_flow, & "Maximum total tidal volume flux.", & - units="m3 s-1", default=3.0d6, scale=US%m_s_to_L_T*US%m_to_L*US%m_to_Z) + units="m3 s-1", default=3.0e6, scale=US%m_s_to_L_T*US%m_to_L*US%m_to_Z) ! Register the open boundaries. call register_OBC(casename, param_file, OBC_Reg)