diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index e2a557daff..7b32270b4c 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/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 "//& diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 394cf05285..ee498f4184 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/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 diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 0a4e12c647..bb9743bb84 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/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 !---- 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 03f2b1ed9d..39e1046c1c 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/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,6 +242,7 @@ 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. @@ -249,7 +250,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! 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 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 e8673da6f8..c704214930 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -21,6 +21,7 @@ module MOM_surface_forcing_nuopc use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS @@ -80,7 +81,7 @@ module MOM_surface_forcing_nuopc !! the correction for the atmospheric (and sea-ice) !! pressure limited by max_p_surf instead of the !! full atmospheric pressure. The default is true. - + logical :: use_CFC !< enables the MOM_CFC_cap tracer package. real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa] logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied !! from an input file. @@ -129,6 +130,7 @@ module MOM_surface_forcing_nuopc type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing character(len=200) :: inputdir !< directory where NetCDF input files are + character(len=200) :: CFC_BC_file !< filename with cfc11 and cfc12 data character(len=200) :: salt_restore_file !< filename for salt restoring data character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file logical :: mask_srestore !< if true, apply a 2-dimensional mask to the surface @@ -142,9 +144,13 @@ module MOM_surface_forcing_nuopc !! temperature restoring fluxes. The masking file should be !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' + character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file + character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_cfc11_atm = -1 !< id number for time_interp_external. + integer :: id_cfc12_atm = -1 !< id number for time_interp_external. ! Diagnostics handles type(forcing_diags), public :: handles @@ -179,11 +185,14 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: frunoff_hflx =>NULL() !< heat content of frozen runoff [W/m2] real, pointer, dimension(:,:) :: p =>NULL() !< pressure of overlying ice and atmosphere !< on ocean surface [Pa] + real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< fractional ice area [nondim] + real, pointer, dimension(:,:) :: u10_sqr =>NULL() !< wind speed squared at 10m [m2/s2] real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] real, pointer, dimension(:,:) :: ice_rigidity =>NULL() !< rigidity of the sea ice, sea-ice and !! 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] @@ -235,6 +244,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! local variables real, dimension(SZI_(G),SZJ_(G)) :: & + cfc11_atm, & !< CFC11 concentration in the atmopshere [???????] + cfc12_atm, & !< CFC11 concentration in the atmopshere [???????] data_restore, & !< The surface value toward which to restore [g/kg or degC] SST_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C] SSS_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg] @@ -294,7 +305,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! flux type has been used. if (fluxes%dt_buoy_accum < 0) then call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & - press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug) + press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & + cfc=CS%use_CFC) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -434,6 +446,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo endif + ! Check that liquid runoff has a place to go if (CS%liquid_runoff_from_data .and. .not. associated(IOB%lrunoff)) then call MOM_error(FATAL, "liquid runoff is being added via data_override but "// & @@ -533,8 +546,27 @@ 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 + ! 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 @@ -551,6 +583,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure. endif + ! CFCs + if (CS%use_CFC) then + call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + endif + if (associated(IOB%salt_flux)) then do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*IOB%salt_flux(i-i0,j-j0)) @@ -689,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 @@ -1042,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 @@ -1055,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]. @@ -1162,6 +1204,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "coupler. This is used for testing and should be =1.0 for any "//& "production runs.", default=1.0) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, & + default=.false., do_not_log=.true.) + if (restore_salt) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//& @@ -1174,7 +1219,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "The name of the surface salinity variable to read from "//& "SALT_RESTORE_FILE for restoring salinity.", & default="salt") - call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & "If true, the restoring of salinity is applied as a salt "//& "flux instead of as a freshwater flux.", default=.false.) @@ -1270,7 +1314,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, enddo ; enddo endif - call time_interp_external_init + ! initialize time interpolator module + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global ! constant. @@ -1321,8 +1366,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, & "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, 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 "//& @@ -1356,6 +1402,29 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, endif endif ; endif + ! Do not log these params here since they are logged in the CFC cap module + if (CS%use_CFC) then + call get_param(param_file, mdl, "CFC_BC_FILE", CS%CFC_BC_file, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ", do_not_log=.true.) + if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then + ! Add the directory if CFC_BC_file is not already a complete path. + CS%CFC_BC_file = trim(CS%inputdir) // trim(CS%CFC_BC_file) + endif + if (len_trim(CS%CFC_BC_file) > 0) then + call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.) + call get_param(param_file, mdl, "CFC12_VARIABLE", CS%cfc12_var_name, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) + + CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) + CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) + endif + endif + ! Set up any restart fields associated with the forcing. call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") call restart_init_end(CS%restart_CSp) @@ -1402,7 +1471,6 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) intent(in) :: iobt !< An ice-ocean boundary type with fluxes to drive the !! ocean in a coupled model whose checksums are reported - ! local variables ! Local variables integer(kind=int64) :: chks ! A checksum for the field logical :: root ! True only on the root PE @@ -1429,6 +1497,12 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) chks = field_chksum( iobt%lrunoff ) ; if (root) write(outunit,100) 'iobt%lrunoff ', chks chks = field_chksum( iobt%frunoff ) ; if (root) write(outunit,100) 'iobt%frunoff ', chks chks = field_chksum( iobt%p ) ; if (root) write(outunit,100) 'iobt%p ', chks + if (associated(iobt%ice_fraction)) then + chks = field_chksum( iobt%ice_fraction ) ; if (root) write(outunit,100) 'iobt%ice_fraction ', chks + endif + if (associated(iobt%u10_sqr)) then + chks = field_chksum( iobt%u10_sqr ) ; if (root) write(outunit,100) 'iobt%u10_sqr ', chks + endif if (associated(iobt%ustar_berg)) then chks = field_chksum( iobt%ustar_berg ) ; if (root) write(outunit,100) 'iobt%ustar_berg ', chks endif @@ -1438,6 +1512,7 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) if (associated(iobt%mass_berg)) then chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks endif + 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 index 679f147797..7c778d9753 100644 --- a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 @@ -243,7 +243,8 @@ subroutine MESO_surface_forcing_init(Time, G, US, param_file, diag, CS) "parameters from vertical units of m to kg m-2.", & units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) 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) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 682ad03397..2429ce9d2d 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -183,6 +183,16 @@ module MOM_forcing_type real :: C_p !< heat capacity of seawater [Q degC-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. + ! CFC-related arrays needed in the MOM_CFC_cap module + real, pointer, dimension(:,:) :: & + cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> mol m-2 s-1]. + cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> mol m-2 s-1]. + 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. @@ -246,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] @@ -348,6 +357,12 @@ module MOM_forcing_type integer :: id_TKE_tidal = -1 integer :: id_buoy = -1 + ! cfc-related diagnostics handles + integer :: id_cfc11 = -1 + integer :: id_cfc12 = -1 + integer :: id_ice_fraction = -1 + integer :: id_u10_sqr = -1 + ! iceberg diagnostic handles integer :: id_ustar_berg = -1 integer :: id_area_berg = -1 @@ -356,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 @@ -1087,6 +1105,14 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%p_surf)) & call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift , scale=US%RL2_T2_to_Pa) + if (associated(fluxes%u10_sqr)) & + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%L_to_m**2*US%s_to_T**2) + if (associated(fluxes%ice_fraction)) & + call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) + if (associated(fluxes%cfc11_flux)) & + call hchksum(fluxes%cfc11_flux, mesg//" fluxes%cfc11_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + if (associated(fluxes%cfc12_flux)) & + call hchksum(fluxes%cfc12_flux, mesg//" fluxes%cfc12_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & @@ -1100,7 +1126,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%heat_content_lrunoff)) & call hchksum(fluxes%heat_content_lrunoff, mesg//" fluxes%heat_content_lrunoff", G%HI, & - haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_frunoff)) & call hchksum(fluxes%heat_content_frunoff, mesg//" fluxes%heat_content_frunoff", G%HI, & haloshift=hshift, scale=US%QRZ_T_to_W_m2) @@ -1239,13 +1265,15 @@ 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, 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 handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1288,6 +1316,34 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, endif endif + ! units for cfc11_flux and cfc12_flux are mol m-2 s-1 + ! See: + ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html + ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html + if (present(use_cfcs)) then + if (use_cfcs) then + handles%id_cfc11 = register_diag_field('ocean_model', 'cfc11_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC11 into the ocean ', 'mol m-2 s-1', & + conversion= US%Z_to_m*US%s_to_T,& + cmor_field_name='fgcfc11', & + cmor_long_name='Surface Downward CFC11 Flux', & + cmor_standard_name='surface_downward_cfc11_flux') + + handles%id_cfc12 = register_diag_field('ocean_model', 'cfc12_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC12 into the ocean ', 'mol m-2 s-1', & + conversion= US%Z_to_m*US%s_to_T,& + cmor_field_name='fgcfc12', & + cmor_long_name='Surface Downward CFC12 Flux', & + cmor_standard_name='surface_downward_cfc12_flux') + + handles%id_ice_fraction = register_diag_field('ocean_model', 'ice_fraction', diag%axesT1, Time, & + 'Fraction of cell area covered by sea ice', 'm2 m-2') + + handles%id_u10_sqr = register_diag_field('ocean_model', 'u10_sqr', diag%axesT1, Time, & + 'Wind magnitude at 10m, squared', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + endif + endif + handles%id_psurf = register_diag_field('ocean_model', 'p_surf', diag%axesT1, Time, & 'Pressure at ice-ocean or atmosphere-ocean interface', & 'Pa', conversion=US%RL2_T2_to_Pa, cmor_field_name='pso', & @@ -1895,6 +1951,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 @@ -2834,6 +2898,19 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_netFWGlobalScl > 0) & call post_data(handles%id_netFWGlobalScl, fluxes%netFWGlobalScl, diag) + ! post diagnostics related to cfcs ==================================== + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc11_flux)) & + call post_data(handles%id_cfc11, fluxes%cfc11_flux, diag) + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc12_flux)) & + call post_data(handles%id_cfc12, fluxes%cfc12_flux, diag) + + if ((handles%id_ice_fraction > 0) .and. associated(fluxes%ice_fraction)) & + call post_data(handles%id_ice_fraction, fluxes%ice_fraction, diag) + + if ((handles%id_u10_sqr > 0) .and. associated(fluxes%u10_sqr)) & + call post_data(handles%id_u10_sqr, fluxes%u10_sqr, diag) ! remaining boundary terms ================================================== @@ -2858,6 +2935,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) @@ -2872,7 +2953,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) + 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 @@ -2884,6 +2965,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: salt !< If present and true, allocate salt fluxes 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 @@ -2939,6 +3022,16 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%area_berg,isd,ied,jsd,jed, iceberg) call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) + !These fields should only on allocated when USE_CFC_CAP is activated. + call myAlloc(fluxes%cfc11_flux,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%cfc12_flux,isd,ied,jsd,jed, cfc) + 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 @@ -3034,14 +3127,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)) @@ -3185,6 +3285,10 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%ustar_berg)) deallocate(fluxes%ustar_berg) if (associated(fluxes%area_berg)) deallocate(fluxes%area_berg) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) + if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) + if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) + if (associated(fluxes%cfc11_flux)) deallocate(fluxes%cfc11_flux) + if (associated(fluxes%cfc12_flux)) deallocate(fluxes%cfc12_flux) call coupler_type_destructor(fluxes%tr_fluxes) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 86493aad93..08f8dea634 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -11,7 +11,7 @@ module MOM_unit_tests use MOM_diag_vkernels, only : diag_vkernels_unit_tests use MOM_random, only : random_unit_tests use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests - +use MOM_CFC_cap, only : CFC_cap_unit_tests implicit none ; private public unit_tests @@ -41,6 +41,8 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: random_unit_tests FAILED") if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: near_boundary_unit_tests FAILED") + if (CFC_cap_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: CFC_cap_unit_tests FAILED") endif end subroutine unit_tests diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f966ab2ad2..86b373b8dd 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -42,6 +42,8 @@ module MOM_variables SST, & !< The sea surface temperature [degC]. SSS, & !< The sea surface salinity [ppt ~> psu or gSalt/kg]. sfc_density, & !< The mixed layer density [R ~> kg m-3]. + sfc_cfc11, & !< Sea surface concentration of CFC11 [mol kg-1]. + sfc_cfc12, & !< Sea surface concentration of CFC12 [mol kg-1]. Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. @@ -300,7 +302,8 @@ module MOM_variables !> Allocates the fields for the surface (return) properties of !! the ocean model. Unused fields are unallocated. subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & - gas_fields_ocn, use_meltpot, use_iceshelves, omit_frazil) + gas_fields_ocn, use_meltpot, use_iceshelves, & + omit_frazil, use_cfcs) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(surface), intent(inout) :: sfc_state !< ocean surface state type to be allocated. logical, optional, intent(in) :: use_temperature !< If true, allocate the space for thermodynamic variables. @@ -313,13 +316,14 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. logical, optional, intent(in) :: use_meltpot !< If true, allocate the space for melt potential + logical, optional, intent(in) :: use_cfcs !< If true, allocate the space for cfcs logical, optional, intent(in) :: use_iceshelves !< If true, allocate the space for the stresses !! under ice shelves. logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to !! pass frazil fluxes to the coupler ! local variables - logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil + logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_cfcs integer :: is, ie, js, je, isd, ied, jsd, jed integer :: isdB, iedB, jsdB, jedB @@ -330,6 +334,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot + alloc_cfcs = .false. ; if (present(use_cfcs)) alloc_cfcs = use_cfcs alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil @@ -353,6 +358,11 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & allocate(sfc_state%melt_potential(isd:ied,jsd:jed)) ; sfc_state%melt_potential(:,:) = 0.0 endif + if (alloc_cfcs) then + allocate(sfc_state%sfc_cfc11(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc11(:,:) = 0.0 + allocate(sfc_state%sfc_cfc12(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc12(:,:) = 0.0 + endif + if (alloc_integ) then ! Allocate structures for the vertically integrated ocean_mass, ocean_heat, and ocean_salt. allocate(sfc_state%ocean_mass(isd:ied,jsd:jed)) ; sfc_state%ocean_mass(:,:) = 0.0 @@ -396,7 +406,8 @@ subroutine deallocate_surface_state(sfc_state) if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat) if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt) if (allocated(sfc_state%salt_deficit)) deallocate(sfc_state%salt_deficit) - + if (allocated(sfc_state%sfc_cfc11)) deallocate(sfc_state%sfc_cfc11) + if (allocated(sfc_state%sfc_cfc12)) deallocate(sfc_state%sfc_cfc12) call coupler_type_destructor(sfc_state%tr_fields) sfc_state%arrays_allocated = .false. diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 762b2edaea..25946eb631 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -81,7 +81,7 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] - real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [T-1 ~> s-1]. logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. @@ -340,6 +340,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif + if (CS%debug) then + call hchksum(src, "MEKE src", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) + endif + ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -534,6 +538,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif endif ! MEKE_KH>=0 + if (CS%debug) then + call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, scale=US%L_T_to_m_s**2) + endif + ! do j=js,je ; do i=is,ie ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) ! enddo ; enddo @@ -688,7 +696,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 else FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points @@ -824,7 +832,7 @@ subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 enddo ; enddo if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) @@ -957,7 +965,7 @@ subroutine MEKE_lengthScales_0d(CS, US, area, beta, depth, Rd_dx, SN, EKE, & ! Z Leady = 0. endif if (CS%use_min_lscale) then - LmixScale = 1.e7 + LmixScale = 1.e7*US%m_to_L if (CS%aDeform*Ldeform > 0.) LmixScale = min(LmixScale,CS%aDeform*Ldeform) if (CS%aFrict *Lfrict > 0.) LmixScale = min(LmixScale,CS%aFrict *Lfrict) if (CS%aRhines*Lrhines > 0.) LmixScale = min(LmixScale,CS%aRhines*Lrhines) @@ -1069,7 +1077,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) if (CS%MEKE_equilibrium_restoring) then call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & - default=1e6, scale=US%T_to_s) + default=1e6, scale=US%s_to_T) CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index fb70f5d679..dabeb27159 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -506,10 +506,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency - !! at u-points [T-2 ~> s-2] + !! at u-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency - !! at v-points [T-2 ~> s-2] + !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), pointer :: CS !< Variable mixing coefficients type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -654,9 +654,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C endif if (CS%debug) then - call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) + call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, & + scale=US%Z_to_L, haloshift=1) call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, & - scale=US%s_to_T**2, scalar_pair=.true.) + scale=US%L_to_Z**2 * US%s_to_T**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & scale=US%s_to_T, scalar_pair=.true.) endif @@ -1484,7 +1485,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB)) ; CS%Depth_fn_v(:,:) = 0.0 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, & "The depth above which KHTH is scaled away.",& - units="m", default=1000.) + units="m", scale=US%m_to_Z, default=1000.) call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, & "The exponent used in the depth dependent scaling function for KHTH.",& units="nondim", default=3.0) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index da62ffc6b7..4d038ed31e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -433,7 +433,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m) if (use_stored_slopes) then call uvchksum("VarMix%slope_[xy]", VarMix%slope_x, VarMix%slope_y, & - G%HI, haloshift=0) + G%HI, haloshift=0, scale=US%Z_to_L) endif if (associated(tv%eqn_of_state)) then call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1) @@ -505,7 +505,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(GV%m_to_H,GV%Z_to_H*G%bathyT(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 083f5ed000..4dcaa70bc2 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 @@ -408,17 +411,27 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) "Unrecognized KPP_LT_K_SHAPE option: "//trim(string)) end select call get_param(paramFile, mdl, "KPP_LT_K_METHOD", string , & - 'Method to enhance mixing coefficient in KPP. '// & + 'Method to enhance mixing coefficient in KPP. '// & 'Valid options are: \n'// & '\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'// & '\t VR12 = Function of Langmuir number based on VR12\n'// & - '\t RW16 = Function of Langmuir number based on RW16', & + '\t (Van Roekel et al. 2012)\n'// & + '\t (Li et al. 2016, OM) \n'// & + '\t RW16 = Function of Langmuir number based on RW16\n'// & + '\t (Reichl et al., 2016, JPO)', & 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 @@ -433,20 +446,33 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) 'in Bulk Richardson Number.', units="", Default=.false.) if (CS%LT_Vt2_Enhancement) then call get_param(paramFile, mdl, "KPP_LT_VT2_METHOD",string , & - 'Method to enhance Vt2 in KPP. '// & + 'Method to enhance Vt2 in KPP. '// & 'Valid options are: \n'// & '\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'// & '\t VR12 = Function of Langmuir number based on VR12\n'// & + '\t (Van Roekel et al., 2012) \n'// & + '\t (Li et al. 2016, OM) \n'// & '\t RW16 = Function of Langmuir number based on RW16\n'// & - '\t LF17 = Function of Langmuir number based on LF17', & + '\t (Reichl et al., 2016, JPO) \n'// & + '\t LF17 = Function of Langmuir number based on LF17\n'// & + '\t (Li and Fox-Kemper, 2017, JPO)', & 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 +482,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 +498,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 +629,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 +650,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 @@ -644,7 +674,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%debug) then call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m) call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T) - call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0) + call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) endif @@ -661,7 +691,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 +744,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 +802,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 +815,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 +865,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 +912,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 +925,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 +939,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 +979,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 @@ -967,8 +992,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%debug) then call hchksum(Salt, "KPP in: S",G%HI,haloshift=0) call hchksum(Temp, "KPP in: T",G%HI,haloshift=0) - call hchksum(u, "KPP in: u",G%HI,haloshift=0) - call hchksum(v, "KPP in: v",G%HI,haloshift=0) + call hchksum(u, "KPP in: u",G%HI,haloshift=0,scale=US%L_T_to_m_s) + call hchksum(v, "KPP in: v",G%HI,haloshift=0,scale=US%L_T_to_m_s) endif call cpu_clock_begin(id_clock_KPP_compute_BLD) @@ -987,7 +1012,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 +1171,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 +1232,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 +1286,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_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 8b007a0b11..a615c9f40b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -288,8 +288,8 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) ! if (CS%id_kv_conv > 0) & ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kd)) call hchksum(Kv, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) + if (present(Kd)) call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) endif ! send diagnostics to post_data diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index fb759b13b2..526eb569af 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -628,11 +628,19 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim 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, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, 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, 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) @@ -1208,11 +1216,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, & + 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, & + 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) @@ -1794,11 +1810,19 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e enddo ; enddo ; enddo endif - 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, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, 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, 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/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 new file mode 100644 index 0000000000..e1770b0d52 --- /dev/null +++ b/src/tracer/MOM_CFC_cap.F90 @@ -0,0 +1,691 @@ +!> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover +!! provided via cap (only NUOPC cap is implemented so far). +module MOM_CFC_cap + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_hor_index, only : hor_index_type +use MOM_grid, only : ocean_grid_type +use MOM_io, only : file_exists, MOM_read_data, slasher +use MOM_io, only : vardesc, var_desc, query_vardesc, stdout +use MOM_open_boundary, only : ocean_OBC_type +use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_time_manager, only : time_type +use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public register_CFC_cap, initialize_CFC_cap, CFC_cap_unit_tests +public CFC_cap_column_physics, CFC_cap_surface_state, CFC_cap_fluxes +public CFC_cap_stock, CFC_cap_end + +integer, parameter :: NTR = 2 !< the number of tracers in this module. + +!> The control structure for the CFC_cap tracer package +type, public :: CFC_cap_CS ; private + character(len=200) :: IC_file !< The file in which the CFC initial values can + !! be found, or an empty string for internal initilaization. + logical :: Z_IC_file !< If true, the IC_file is in Z-space. The default is false. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the MOM6 tracer registry + real, pointer, dimension(:,:,:) :: & + CFC11 => NULL(), & !< The CFC11 concentration [mol kg-1]. + CFC12 => NULL() !< The CFC12 concentration [mol kg-1]. + ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12. + real :: CFC11_IC_val = 0.0 !< The initial value assigned to CFC11 [mol kg-1]. + real :: CFC12_IC_val = 0.0 !< The initial value assigned to CFC12 [mol kg-1]. + real :: CFC11_land_val = -1.0 !< The value of CFC11 used where land is masked out [mol kg-1]. + real :: CFC12_land_val = -1.0 !< The value of CFC12 used where land is masked out [mol kg-1]. + logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code + !! if they are not found in the restart files. + character(len=16) :: CFC11_name !< CFC11 variable name + character(len=16) :: CFC12_name !< CFC12 variable name + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate + !! the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Model restart control structure + + ! The following vardesc types contain a package of metadata about each tracer. + type(vardesc) :: CFC11_desc !< A set of metadata for the CFC11 tracer + type(vardesc) :: CFC12_desc !< A set of metadata for the CFC12 tracer + !>@{ Diagnostic IDs + integer :: id_cfc11_cmor = -1, id_cfc12_cmor = -1 + !>@} +end type CFC_cap_CS + +contains + +!> Register the CFCs to be used with MOM and read the parameters +!! that are used with this tracer package +function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(CFC_cap_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module. + type(tracer_registry_type), & + pointer :: tr_Reg !< A pointer to the tracer registry. + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. + + ! Local variables + character(len=40) :: mdl = "MOM_CFC_cap" ! This module's name. + character(len=200) :: inputdir ! The directory where NetCDF input files are. + ! This include declares and sets the variable "version". +#include "version_variable.h" + real, dimension(:,:,:), pointer :: tr_ptr => NULL() + character(len=200) :: dummy ! Dummy variable to store params that need to be logged here. + logical :: register_CFC_cap + integer :: isd, ied, jsd, jed, nz, m + + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(WARNING, "register_CFC_cap called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "CFC_IC_FILE", CS%IC_file, & + "The file in which the CFC initial values can be "//& + "found, or an empty string for internal initialization.", & + default=" ") + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + ! Add the directory if CS%IC_file is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file) + endif + call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", CS%Z_IC_file, & + "If true, CFC_IC_FILE is in depth space, not layer space", & + default=.false.) + call get_param(param_file, mdl, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + call get_param(param_file, mdl, "CFC11_IC_VAL", CS%CFC11_IC_val, & + "Value that CFC_11 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) + call get_param(param_file, mdl, "CFC12_IC_VAL", CS%CFC12_IC_val, & + "Value that CFC_12 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) + + ! the following params are not used in this module. Instead, they are used in + ! the cap but are logged here to keep all the CFC cap params together. + call get_param(param_file, mdl, "CFC_BC_FILE", dummy, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ") + if ((len_trim(dummy) > 0) .and. (scan(dummy,'/') == 0)) then + ! Add the directory if dummy is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + dummy = trim(slasher(inputdir))//trim(dummy) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", dummy) + endif + if (len_trim(dummy) > 0) then + call get_param(param_file, mdl, "CFC11_VARIABLE", dummy, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11") + call get_param(param_file, mdl, "CFC12_VARIABLE", dummy, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12") + endif + + ! The following vardesc types contain a package of metadata about each tracer, + ! including, the name; units; longname; and grid information. + CS%CFC11_name = "CFC_11" ; CS%CFC12_name = "CFC_12" + CS%CFC11_desc = var_desc(CS%CFC11_name,"mol kg-1","Moles Per Unit Mass of CFC-11 in sea water", caller=mdl) + CS%CFC12_desc = var_desc(CS%CFC12_name,"mol kg-1","Moles Per Unit Mass of CFC-12 in sea water", caller=mdl) + + allocate(CS%CFC11(isd:ied,jsd:jed,nz)) ; CS%CFC11(:,:,:) = 0.0 + allocate(CS%CFC12(isd:ied,jsd:jed,nz)) ; CS%CFC12(:,:,:) = 0.0 + + ! This pointer assignment is needed to force the compiler not to do a copy in + ! the registration calls. Curses on the designers and implementers of F90. + tr_ptr => CS%CFC11 + ! Register CFC11 for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC11_desc, registry_diags=.true., & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + ! Do the same for CFC12 + tr_ptr => CS%CFC12 + call register_tracer(tr_ptr, Tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC12_desc, registry_diags=.true., & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + register_CFC_cap = .true. + +end function register_CFC_cap + +!> Initialize the CFC tracer fields and set up the tracer output. +subroutine initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate + !! diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type + !! specifies whether, where, and what + !! open boundary conditions are used. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: from_file = .false. + + if (.not.associated(CS)) return + + CS%Time => day + CS%diag => diag + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC11, CS%CFC11_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC11, CS%CFC11_name, CS%CFC11_land_val, & + CS%CFC11_IC_val, G, GV, US, CS) + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC12, CS%CFC12_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC12, CS%CFC12_name, CS%CFC12_land_val, & + CS%CFC12_IC_val, G, GV, US, CS) + + + ! cmor diagnostics + ! CFC11 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/42625c97b8fe75124a345962c4430982.html + CS%id_cfc11_cmor = register_diag_field('ocean_model', 'cfc11', diag%axesTL, day, & + 'Mole Concentration of CFC11 in Sea Water', 'mol m-3') + ! CFC12 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/3ab8e10027d7014f18f9391890369235.html + CS%id_cfc12_cmor = register_diag_field('ocean_model', 'cfc12', diag%axesTL, day, & + 'Mole Concentration of CFC12 in Sea Water', 'mol m-3') + + if (associated(OBC)) then + ! Steal from updated DOME in the fullness of time. + ! GMM: TODO this must be coded if we intend to use this module in regional applications + endif + +end subroutine initialize_CFC_cap + +!>This subroutine initializes a tracer array. +subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, GV, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: tr !< The tracer concentration array + character(len=*), intent(in) :: name !< The tracer name + real, intent(in) :: land_val !< A value the tracer takes over land + real, intent(in) :: IC_val !< The initial condition value for the tracer + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: OK + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file, G%Domain)) & + call MOM_error(FATAL, "initialize_CFC_cap: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr, h, CS%IC_file, name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr, h, CS%IC_file, trim(name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_CFC_cap: "//& + "Unable to read "//trim(name)//" from "//& + trim(CS%IC_file)//".") + endif + else + call MOM_read_data(CS%IC_file, trim(name), tr, G%Domain) + endif + else + do k=1,nz ; do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) < 0.5) then + tr(i,j,k) = land_val + else + tr(i,j,k) = IC_val + endif + enddo ; enddo ; enddo + endif + +end subroutine init_tracer_CFC + +!> Applies diapycnal diffusion, souces and sinks and any other column +!! tracer physics to the CFC cap tracers. CFCs are relatively simple, +!! as they are passive tracers with only a surface flux as a source. +subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, & + evap_CFL_limit, minimum_forcing_depth) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes!< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [H ~> m or kg m-2] + + ! The arguments to this subroutine are redundant in that + ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + + ! Local variables + real, pointer, dimension(:,:,:) :: CFC11 => NULL(), CFC12 => NULL() + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] + real :: scale_factor ! convert from cfc1[12]_flux to units of sfc_flux in tracer_vertdiff + integer :: i, j, k, m, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (.not.associated(CS)) return + + ! set factor to convert from cfc1[12]_flux units to tracer_vertdiff argument sfc_flux units + ! cfc1[12]_flux units are CU Z T-1 kg m-3 + ! tracer_vertdiff argument sfc_flux units are CU kg m-2 T-1 + scale_factor = US%Z_to_m + + CFC11 => CS%CFC11 ; CFC12 => CS%CFC12 + + ! Use a tridiagonal solver to determine the concentrations after the + ! surface source is applied and diapycnal advection and diffusion occurs. + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC11, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC12, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + else + call tracer_vertdiff(h_old, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + call tracer_vertdiff(h_old, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + endif + + ! If needed, write out any desired diagnostics from tracer sources & sinks here. + if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC11, CS%diag) + if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC12, CS%diag) + +end subroutine CFC_cap_column_physics + +!> Calculates the mass-weighted integral of all tracer stocks, +!! returning the number of stocks it has calculated. If the stock_index +!! is present, only the stock corresponding to that coded index is returned. +function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of each + !! tracer, in kg times concentration units [kg conc]. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< The coded index of a specific + !! stock being sought. + integer :: CFC_cap_stock !< The number of stocks calculated here. + + ! Local variables + real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or nondim] + real :: mass ! The cell volume or mass [H L2 ~> m3 or kg] + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + CFC_cap_stock = 0 + if (.not.associated(CS)) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + call query_vardesc(CS%CFC11_desc, name=names(1), units=units(1), caller="CFC_cap_stock") + call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="CFC_cap_stock") + units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg" + + stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2 + stocks(1) = 0.0 ; stocks(2) = 0.0 + do k=1,nz ; do j=js,je ; do i=is,ie + mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k) + stocks(1) = stocks(1) + CS%CFC11(i,j,k) * mass + stocks(2) = stocks(2) + CS%CFC12(i,j,k) * mass + enddo ; enddo ; enddo + stocks(1) = stock_scale * stocks(1) + stocks(2) = stock_scale * stocks(2) + + CFC_cap_stock = 2 + +end function CFC_cap_stock + +!> Extracts the ocean surface CFC concentrations and copies them to sfc_state. +subroutine CFC_cap_surface_state(sfc_state, G, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + type(CFC_cap_CS), pointer :: CS!< The control structure returned by a previous + !! call to register_CFC_cap. + + ! Local variables + integer :: i, j, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + if (.not.associated(CS)) return + + do j=js,je ; do i=is,ie + sfc_state%sfc_cfc11(i,j) = CS%CFC11(i,j,1) + sfc_state%sfc_cfc12(i,j) = CS%CFC12(i,j,1) + enddo ; enddo + +end subroutine CFC_cap_surface_state + +!> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM +!! concentration, and calculating the solubility, Schmidt number, and gas exchange. +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) + type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type + type(surface), intent(in ) :: sfc_state !< A structure containing fields + !! that describe the surface state of the ocean. + type(forcing), intent(inout) :: fluxes !< A structure containing pointers + !! to thermodynamic and tracer forcing fields. Unused fields + !! have NULL ptrs. + real, intent(in ) :: Rho0 !< The mean ocean density [R ~> kg m-3] + type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the + !! CFC's concentration in the atmosphere. + integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. + integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: & + kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1]. + kw, & ! gas transfer velocity [Z T-1 ~> m s-1]. + cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) + ! [mol kg-1]. + cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] + cfc12_atm !< CFC11 concentration in the atmopshere [pico mol/mol] + real :: ta ! Absolute sea surface temperature [hectoKelvin] + real :: sal ! Surface salinity [PSU]. + real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. + real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. + real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12. + real :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + real :: Rho0_kg_m3 ! Rho0 in non-scaled units [kg m-3] + real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. + real :: press_to_atm ! converts from model pressure units to atm + integer :: i, j, m, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! CFC11 ATM concentration + if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then + call time_interp_external(id_cfc11_atm, Time, cfc11_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc11_atm = cfc11_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc11_atm internally" //& + "has not been implemented yet.") + endif + + ! CFC12 ATM concentration + if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then + call time_interp_external(id_cfc12_atm, Time, cfc12_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc12_atm = cfc12_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc12_atm internally" //& + "has not been implemented yet.") + endif + + !--------------------------------------------------------------------- + ! Gas exchange/piston velocity parameter + !--------------------------------------------------------------------- + ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 + ! = 6.97e-7 m/s s^2/m^2 [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + kw_coeff = (US%m_to_Z*US%s_to_T*US%L_to_m**2) * 6.97e-7 + + ! set unit conversion factors + Rho0_kg_m3 = Rho0 * US%R_to_kg_m3 + press_to_atm = US%R_to_kg_m3*US%L_T_to_m_s**2 * pa_to_atm + + do j=js,je ; do i=is,ie + ! ta in hectoKelvin + ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) + sal = sfc_state%SSS(i,j) + + ! Calculate solubilities + call get_solubility(alpha_11, alpha_12, ta, sal , G%mask2dT(i,j)) + + ! Calculate Schmidt numbers using coefficients given by + ! Wanninkhof (2014); doi:10.4319/lom.2014.12.351. + call comp_CFC_schmidt(sfc_state%SST(i,j), sc_11, sc_12) + + kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) + + ! air concentrations and cfcs BC's fluxes + ! CFC flux units: CU Z T-1 kg m-3 = mol kg-1 Z T-1 kg m-3 ~> mol m-2 s-1 + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_11) + cair(i,j) = press_to_atm * alpha_11 * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC11(i,j)) * Rho0_kg_m3 + + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_12) + cair(i,j) = press_to_atm * alpha_12 * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC12(i,j)) * Rho0_kg_m3 + enddo ; enddo + +end subroutine CFC_cap_fluxes + +!> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32. +subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask) + real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1] + real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1] + real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin] + real, intent(in ) :: sal !< Surface salinity [PSU]. + real, intent(in ) :: mask !< ocean mask + + ! Local variables + + ! Coefficients for calculating CFC11 solubilities + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + + real, parameter :: d1_11 = -232.0411 ! [nondim] + real, parameter :: d2_11 = 322.5546 ! [hectoKelvin-1] + real, parameter :: d3_11 = 120.4956 ! [log(hectoKelvin)-1] + real, parameter :: d4_11 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_11 = -0.146531 ! [PSU-1] + real, parameter :: e2_11 = 0.093621 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_11 = -0.0160693 ! [PSU-2 hectoKelvin-2] + + ! Coefficients for calculating CFC12 solubilities + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + + real, parameter :: d1_12 = -220.2120 ! [nondim] + real, parameter :: d2_12 = 301.8695 ! [hectoKelvin-1] + real, parameter :: d3_12 = 114.8533 ! [log(hectoKelvin)-1] + real, parameter :: d4_12 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_12 = -0.147718 ! [PSU-1] + real, parameter :: e2_12 = 0.093175 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_12 = -0.0157340 ! [PSU-2 hectoKelvin-2] + + real :: factor ! introduce units to result [mol kg-1 atm-1] + + ! Eq. 9 from Warner and Weiss (1985) DSR, vol 32. + factor = 1.0 + alpha_11 = exp(d1_11 + d2_11/ta + d3_11*log(ta) + d4_11*ta**2 +& + sal * ((e3_11 * ta + e2_11) * ta + e1_11)) * & + factor * mask + alpha_12 = exp(d1_12 + d2_12/ta + d3_12*log(ta) + d4_12*ta**2 +& + sal * ((e3_12 * ta + e2_12) * ta + e1_12)) * & + factor * mask + +end subroutine get_solubility + + +!> Compute Schmidt numbers of CFCs following Wanninkhof (2014); doi:10.4319/lom.2014.12.351 +!! Range of validity of fit is -2:40. +subroutine comp_CFC_schmidt(sst_in, cfc11_sc, cfc12_sc) + real, intent(in) :: sst_in !< The sea surface temperature [degC]. + real, intent(inout) :: cfc11_sc !< Schmidt number of CFC11 [nondim]. + real, intent(inout) :: cfc12_sc !< Schmidt number of CFC12 [nondim]. + + !local variables + real , parameter :: a_11 = 3579.2 + real , parameter :: b_11 = -222.63 + real , parameter :: c_11 = 7.5749 + real , parameter :: d_11 = -0.14595 + real , parameter :: e_11 = 0.0011874 + real , parameter :: a_12 = 3828.1 + real , parameter :: b_12 = -249.86 + real , parameter :: c_12 = 8.7603 + real , parameter :: d_12 = -0.1716 + real , parameter :: e_12 = 0.001408 + real :: sst + + + ! clip SST to avoid bad values + sst = MAX(-2.0, MIN(40.0, sst_in)) + cfc11_sc = a_11 + sst * (b_11 + sst * (c_11 + sst * (d_11 + sst * e_11))) + cfc12_sc = a_12 + sst * (b_12 + sst * (c_12 + sst * (d_12 + sst * e_12))) + +end subroutine comp_CFC_schmidt + +!> Deallocate any memory associated with the CFC cap tracer package +subroutine CFC_cap_end(CS) + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + ! local variables + integer :: m + + if (associated(CS)) then + if (associated(CS%CFC11)) deallocate(CS%CFC11) + if (associated(CS%CFC12)) deallocate(CS%CFC12) + + deallocate(CS) + endif +end subroutine CFC_cap_end + +!> Unit tests for the CFC cap module. +logical function CFC_cap_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, output additional + !! information for debugging unit tests + + ! Local variables + real :: dummy1, dummy2, ta, sal + character(len=120) :: test_name ! Title of the unit test + + CFC_cap_unit_tests = .false. + write(stdout,*) '==== MOM_CFC_cap =======================' + + ! test comp_CFC_schmidt, Table 1 in Wanninkhof (2014); doi:10.4319/lom.2014.12.351 + test_name = 'Schmidt number calculation' + call comp_CFC_schmidt(20.0, dummy1, dummy2) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 1179.0, 0.5) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 1188.0, 0.5) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)') "Passed "//test_name + + test_name = 'Solubility function, SST = 1.0 C, and SSS = 10 psu' + ta = max(0.01, (1.0 + 273.15) * 0.01); sal = 10. + ! cfc1 = 3.238 10-2 mol kg-1 atm-1 + ! cfc2 = 7.943 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 3.238e-2, 5.0e-6) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 7.943e-3, 5.0e-6) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + + test_name = 'Solubility function, SST = 20.0 C, and SSS = 35 psu' + ta = max(0.01, (20.0 + 273.15) * 0.01); sal = 35. + ! cfc1 = 0.881 10-2 mol kg-1 atm-1 + ! cfc2 = 2.446 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 8.8145e-3, 5.0e-8) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 2.4462e-3, 5.0e-8) + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + +end function CFC_cap_unit_tests + +!> Test that ans and calc are approximately equal by computing the difference +!! and comparing it against limit. +logical function compare_values(verbose, test_name, calc, ans, limit) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + real, intent(in) :: calc !< computed value + real, intent(in) :: ans !< correct value + real, intent(in) :: limit !< value above which test fails + + ! Local variables + real :: diff + + diff = ans - calc + + compare_values = .false. + if (diff > limit ) then + compare_values = .true. + write(stdout,*) "CFC_cap_unit_tests, UNIT TEST FAILED: ", test_name + write(stdout,10) calc, ans + elseif (verbose) then + write(stdout,10) calc, ans + endif + +10 format("calc=",f20.16," ans",f20.16) +end function compare_values + +!> \namespace mom_CFC_cap +!! +!! This module contains the code that is needed to simulate +!! CFC-11 and CFC-12 using atmospheric and sea ice variables +!! provided via cap (only NUOPC cap is implemented so far). + +end module MOM_CFC_cap diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index 9be4af08dc..c1e39598cc 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -152,7 +152,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) @@ -185,14 +185,14 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ea(i,j,1) b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ea(i,j,k) @@ -200,7 +200,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = eb(i,j,nz-1) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ea(i,j,nz) @@ -208,7 +208,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ea(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo @@ -267,8 +267,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & !! ensure positive definiteness [H ~> m or kg m-2]. real :: h_neglect !< A thickness that is so small it is usually lost !! in roundoff and can be neglected [H ~> m or kg m-2]. - logical :: convert_flux = .true. - + logical :: convert_flux integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -279,6 +278,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & return endif + convert_flux = .true. if (present(convert_flux_in)) convert_flux = convert_flux_in h_neglect = GV%H_subroundoff sink_dist = 0.0 @@ -351,7 +351,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) @@ -384,14 +384,14 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ent(i,j,1) b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ent(i,j,K) @@ -399,7 +399,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = ent(i,j,nz) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ent(i,j,nz) @@ -407,7 +407,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ent(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 26ef197ae2..439e9b5396 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -42,6 +42,9 @@ module MOM_tracer_flow_control use MOM_OCMIP2_CFC, only : register_OCMIP2_CFC, initialize_OCMIP2_CFC, flux_init_OCMIP2_CFC use MOM_OCMIP2_CFC, only : OCMIP2_CFC_column_physics, OCMIP2_CFC_surface_state use MOM_OCMIP2_CFC, only : OCMIP2_CFC_stock, OCMIP2_CFC_end, OCMIP2_CFC_CS +use MOM_CFC_cap, only : register_CFC_cap, initialize_CFC_cap +use MOM_CFC_cap, only : CFC_cap_column_physics, CFC_cap_surface_state +use MOM_CFC_cap, only : CFC_cap_stock, CFC_cap_end, CFC_cap_CS use oil_tracer, only : register_oil_tracer, initialize_oil_tracer use oil_tracer, only : oil_tracer_column_physics, oil_tracer_surface_state use oil_tracer, only : oil_stock, oil_tracer_end, oil_tracer_CS @@ -80,6 +83,7 @@ module MOM_tracer_flow_control logical :: use_oil = .false. !< If true, use the oil tracer package logical :: use_advection_test_tracer = .false. !< If true, use the advection_test_tracer package logical :: use_OCMIP2_CFC = .false. !< If true, use the OCMIP2_CFC tracer package + logical :: use_CFC_cap = .false. !< If true, use the CFC_cap tracer package logical :: use_MOM_generic_tracer = .false. !< If true, use the MOM_generic_tracer packages logical :: use_pseudo_salt_tracer = .false. !< If true, use the psuedo_salt tracer package logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package @@ -94,6 +98,7 @@ module MOM_tracer_flow_control type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL() type(advection_test_tracer_CS), pointer :: advection_test_tracer_CSp => NULL() type(OCMIP2_CFC_CS), pointer :: OCMIP2_CFC_CSp => NULL() + type(CFC_cap_CS), pointer :: CFC_cap_CSp => NULL() type(MOM_generic_tracer_CS), pointer :: MOM_generic_tracer_CSp => NULL() type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL() type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL() @@ -114,7 +119,7 @@ subroutine call_tracer_flux_init(verbosity) type(param_file_type) :: param_file ! A structure to parse for run-time parameters character(len=40) :: mdl = "call_tracer_flux_init" ! This module's name. - logical :: use_OCMIP_CFCs, use_MOM_generic_tracer + logical :: use_OCMIP_CFCs, use_MOM_generic_tracer, use_CFC_caps ! Determine which tracer routines with tracer fluxes are to be called. Note ! that not every tracer package is required to have a flux_init call. @@ -193,6 +198,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_OCMIP2_CFC", CS%use_OCMIP2_CFC, & "If true, use the MOM_OCMIP2_CFC tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC_cap, & + "If true, use the MOM_CFC_cap tracer package.", & + default=.false.) call get_param(param_file, mdl, "USE_generic_tracer", CS%use_MOM_generic_tracer, & "If true and _USE_GENERIC_TRACER is defined as a "//& "preprocessor macro, use the MOM_generic_tracer packages.", & @@ -237,6 +245,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_OCMIP2_CFC) CS%use_OCMIP2_CFC = & register_OCMIP2_CFC(HI, GV, param_file, CS%OCMIP2_CFC_CSp, & tr_Reg, restart_CS) + if (CS%use_CFC_cap) CS%use_CFC_cap = & + register_CFC_cap(HI, GV, param_file, CS%CFC_cap_CSp, & + tr_Reg, restart_CS) if (CS%use_MOM_generic_tracer) CS%use_MOM_generic_tracer = & register_MOM_generic_tracer(HI, GV, param_file, CS%MOM_generic_tracer_CSp, & tr_Reg, restart_CS) @@ -317,6 +328,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag if (CS%use_OCMIP2_CFC) & call initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS%OCMIP2_CFC_CSp, & sponge_CSp) + if (CS%use_CFC_cap) & + call initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS%CFC_cap_CSp) + if (CS%use_MOM_generic_tracer) & call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp) @@ -462,6 +476,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, G, GV, US, CS%OCMIP2_CFC_CSp, & evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -516,6 +535,9 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -624,6 +646,12 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni set_pkg_name, max_ns, ns_tot, stock_names, stock_units) endif + if (CS%use_CFC_cap) then + ns = CFC_cap_stock(h, values, G, GV, CS%CFC_cap_CSp, names, units, stock_index) + call store_stocks("MOM_CFC_cap", ns, names, units, values, index, stock_values, & + set_pkg_name, max_ns, ns_tot, stock_names, stock_units) + endif + if (CS%use_advection_test_tracer) then ns = advection_test_stock( h, values, G, GV, CS%advection_test_tracer_CSp, & names, units, stock_index ) @@ -753,6 +781,8 @@ subroutine call_tracer_surface_state(sfc_state, h, G, GV, CS) call advection_test_tracer_surface_state(sfc_state, h, G, GV, CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_surface_state(sfc_state, h, G, GV, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_surface_state(sfc_state, G, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) & call MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS%MOM_generic_tracer_CSp) @@ -772,6 +802,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) call advection_test_tracer_end(CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) call OCMIP2_CFC_end(CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) call CFC_cap_end(CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) call end_MOM_generic_tracer(CS%MOM_generic_tracer_CSp) if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp) if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp) 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)