From 0d591560846c96bb3b667e2f59bc7324bfe106ed Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Sun, 13 Jun 2021 22:56:47 -0600 Subject: [PATCH 1/5] introduce a new wave coupling method: VR12-MA (not fully implemented yet) --- config_src/drivers/nuopc_cap/mom_cap.F90 | 49 ++++++++++--------- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 16 ++++-- src/user/MOM_wave_interface.F90 | 14 +++++- 3 files changed, 51 insertions(+), 28 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..8d383c3d3f 100644 --- 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 index 493762f4bc..7e7716d386 100644 --- 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 index a8e3e207a6..1ec1e9314c 100644 --- 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 From f45659d47d48ef82fd2c26d7eb94ad49928d037f Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Mon, 14 Jun 2021 20:22:57 -0600 Subject: [PATCH 2/5] introduce the lamult import from ww3 in nuopc cap and add a diag --- config_src/drivers/nuopc_cap/mom_cap.F90 | 17 +++---- .../drivers/nuopc_cap/mom_cap_methods.F90 | 8 ++++ .../nuopc_cap/mom_ocean_model_nuopc.F90 | 7 +-- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 19 ++++++-- src/core/MOM_forcing_type.F90 | 46 +++++++++++++++---- 5 files changed, 71 insertions(+), 26 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 8d383c3d3f..e0aaff8534 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -697,7 +697,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (ocean_state%use_waves) then if (ocean_state%wave_method == "VR12-MA") then - continue + allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) ) + Ice_ocean_boundary%lamult = 0.0 else Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & @@ -722,15 +723,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") - 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 @@ -755,7 +747,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") if (ocean_state%use_waves) then if (ocean_state%wave_method == "VR12-MA") then - continue + 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") 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") diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 0a4e12c647..085dbe818c 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -256,6 +256,14 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- + ! Langmuir enhancement factor + !---- + 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 + !---- ! 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 7e7716d386..a07b24e8d8 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -367,13 +367,16 @@ 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_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. 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) 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, & @@ -385,8 +388,6 @@ 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 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) 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 7168823fbc..61d900f0b6 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -183,6 +183,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] @@ -688,8 +689,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 @@ -849,6 +853,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) endif ! endif for wind related fields ! wave to ocean coupling + if ( associated(IOB%lamult)) then + do j=js,je; do i=is,ie + forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + enddo ; enddo + call pass_var(forces%lamult, G%domain, halo=1 ) + endif + if ( associated(IOB%ustkb) ) then forces%stk_wavenumbers(:) = IOB%stk_wavenumbers @@ -1041,7 +1052,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 @@ -1054,6 +1065,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]. @@ -1321,7 +1334,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "If true, makes available diagnostics of fluxes from icebergs "//& "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_berg_fluxes=iceberg_flux_diags, use_waves=use_waves) 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 682ad03397..6e2c0040c8 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -246,7 +246,8 @@ 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(:,:) :: & + lamult => NULL() !< Langmuir enhancement factor [nondim] real, pointer, dimension(:,:) :: & ustk0 => NULL(), & !< Surface Stokes drift, zonal [m/s] vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] @@ -356,6 +357,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 @@ -1239,13 +1243,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) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_waves) 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 ! Clock for forcing diagnostics handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1895,6 +1900,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 @@ -2266,6 +2279,10 @@ subroutine mech_forcing_diags(forces_in, dt, G, time_end, diag, handles) ! endif + ! wave forcing =============================================================== + if (handles%id_lamult > 0) & + call post_data(handles%id_lamult, forces%lamult, diag) + call disable_averaging(diag) if (turns /= 0) then @@ -3034,14 +3051,25 @@ 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%lamult)) then + allocate(forces%lamult(isd:ied,jsd:jed)) + endif + else !num_stk_bands > 0 + 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)) From 5f756ab7e2312039686fa97e2a642db819c4a2e1 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 18:48:39 -0600 Subject: [PATCH 3/5] Enable the export of ice fraction when wave coupling is on. The export of ice fraction was previously added for the newly added CFC module. --- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 1 + .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 18 +++++++----------- src/core/MOM_forcing_type.F90 | 6 +++++- 3 files changed, 13 insertions(+), 12 deletions(-) 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 1af4185a78..d9e5e53fdd 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -395,6 +395,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif if (OS%use_waves) then + call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true.) 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) if (OS%wave_method == "VR12-MA") then 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 f1429bf867..4b03ccf743 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -545,18 +545,14 @@ 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) - enddo ; enddo + ! 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) - 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) - enddo ; enddo - endif + enddo ; enddo ! applied surface pressure from atmosphere and cryosphere if (associated(IOB%p)) then diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 69114d9339..cd4aa06edd 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -2948,7 +2948,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 @@ -2961,6 +2961,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 @@ -3022,6 +3023,9 @@ 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) + if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group From 7d85ab2912ced10af8bcb10e18df3e0c6aec903a Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 19:49:09 -0600 Subject: [PATCH 4/5] set lamult to 1 below ice --- config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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 4b03ccf743..228b47dc69 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -879,7 +879,11 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) ! wave to ocean coupling if ( associated(IOB%lamult)) then do j=js,je; do i=is,ie - forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + if (IOB%ice_fraction(i-i0,j-j0) <= 0.05 ) then + forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + else + forces%lamult(i,j) = 1.0 + endif enddo ; enddo call pass_var(forces%lamult, G%domain, halo=1 ) endif From 997636e3cd42c1ccc2a7e088d46bce61339dc9db Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 20:54:27 -0600 Subject: [PATCH 5/5] move lamult from forces to fluxes --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 24 +++++++++---------- src/core/MOM_forcing_type.F90 | 20 +++++++--------- 2 files changed, 21 insertions(+), 23 deletions(-) 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 228b47dc69..4743e0d268 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -554,6 +554,18 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo + ! 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 if (associated(IOB%p)) then if (CS%max_p_surf >= 0.0) then @@ -876,18 +888,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) endif ! endif for wind related fields - ! 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 - forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) - else - forces%lamult(i,j) = 1.0 - endif - enddo ; enddo - call pass_var(forces%lamult, G%domain, halo=1 ) - endif - if ( associated(IOB%ustkb) ) then forces%stk_wavenumbers(:) = IOB%stk_wavenumbers diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index cd4aa06edd..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,8 +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(:,:) :: & - lamult => NULL() !< Langmuir enhancement factor [nondim] real, pointer, dimension(:,:) :: & ustk0 => NULL(), & !< Surface Stokes drift, zonal [m/s] vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] @@ -2325,10 +2326,6 @@ subroutine mech_forcing_diags(forces_in, dt, G, time_end, diag, handles) ! endif - ! wave forcing =============================================================== - if (handles%id_lamult > 0) & - call post_data(handles%id_lamult, forces%lamult, diag) - call disable_averaging(diag) if (turns /= 0) then @@ -2934,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) @@ -3025,6 +3026,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) if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group @@ -3126,11 +3128,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & 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%lamult)) then - allocate(forces%lamult(isd:ied,jsd:jed)) - endif - else !num_stk_bands > 0 + if (num_stk_bands > 0) then if (.not.associated(forces%ustkb)) then allocate(forces%stk_wavenumbers(num_stk_bands)) forces%stk_wavenumbers(:) = 0.0