diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index d0ee3aad87..5f4b2e19ca 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -737,6 +737,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary% hrofi (isc:iec,jsc:jec), & Ice_ocean_boundary% hevap (isc:iec,jsc:jec), & Ice_ocean_boundary% hcond (isc:iec,jsc:jec), & + Ice_ocean_boundary% lrunoff_glc (isc:iec,jsc:jec), & + Ice_ocean_boundary% frunoff_glc (isc:iec,jsc:jec), & + Ice_ocean_boundary% hrofl_glc (isc:iec,jsc:jec), & + Ice_ocean_boundary% hrofi_glc (isc:iec,jsc:jec), & source=0.0) ! Needed for MARBL @@ -797,6 +801,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_pslv" , "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 + if (cesm_coupled) then + call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofl_glc" , "will provide") !-> liquid glc runoff + call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofi_glc" , "will provide") !-> frozen glc runoff + endif 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, "Fioi_meltw" , "will provide") @@ -808,6 +816,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hcond" , "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofl" , "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofi" , "will provide") + if (cesm_coupled) then + call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofl_glc" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofi_glc" , "will provide") + endif if (Ice_ocean_boundary%ice_ncat > 0) then call fld_list_add(fldsToOcn_num, fldsToOcn, "Sf_afracr", "will provide") @@ -2855,6 +2867,34 @@ end subroutine shr_log_setLogUnit !! !! !! +!! Forr_rofl_glc +!! kg m-2 s-1 +!! runoff +!! mass flux of liquid glc runoff +!! +!! +!! +!! Forr_rofi_glc +!! kg m-2 s-1 +!! runoff +!! mass flux of frozen glc runoff +!! +!! +!! +!! Foxx_hrofi_glc +!! W m-2 +!! hrofi_glc +!! heat content (enthalpy) of frozen glc runoff +!! +!! +!! +!! Foxx_hrofl_glc +!! W m-2 +!! hrofl_glc +!! heat content (enthalpy) of liquid glc runoff +!! +!! +!! !! Fioi_salt !! kg m-2 s-1 !! salt_flux diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index d5ec9dc259..bb12dc6092 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -216,6 +216,22 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! liquid glc runoff + if ( associated(ice_ocean_boundary%lrunoff_glc) ) then + ice_ocean_boundary%lrunoff_glc (:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Forr_rofl_glc', & + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_glc, areacor=med2mod_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + ! frozen glc runoff + if ( associated(ice_ocean_boundary%frunoff_glc) ) then + ice_ocean_boundary%frunoff_glc (:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Forr_rofi_glc', & + isc, iec, jsc, jec, ice_ocean_boundary%frunoff_glc, areacor=med2mod_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + !---- ! Enthalpy terms !---- @@ -256,6 +272,23 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, if (ChkErr(rc,__LINE__,u_FILE_u)) return end if + !---- + ! enthalpy from liquid glc runoff (hrofl_glc) + !---- + if ( associated(ice_ocean_boundary%hrofl_glc) ) then + call state_getimport(importState, 'Foxx_hrofl_glc', isc, iec, jsc, jec, & + ice_ocean_boundary%hrofl_glc, areacor=med2mod_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + !---- + ! enthalpy from frozen glc runoff (hrofi_glc) + !---- + if ( associated(ice_ocean_boundary%hrofi_glc) ) then + call state_getimport(importState, 'Foxx_hrofi_glc', isc, iec, jsc, jec, & + ice_ocean_boundary%hrofi_glc, areacor=med2mod_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if !---- ! enthalpy from evaporation (hevap) !---- 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 122a9d00ca..897491711f 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -164,6 +164,8 @@ module MOM_surface_forcing_nuopc type, public :: ice_ocean_boundary_type real, pointer, dimension(:,:) :: lrunoff =>NULL() !< liquid runoff [kg/m2/s] real, pointer, dimension(:,:) :: frunoff =>NULL() !< ice runoff [kg/m2/s] + real, pointer, dimension(:,:) :: lrunoff_glc =>NULL() !< liquid glc runoff via rof [kg/m2/s] + real, pointer, dimension(:,:) :: frunoff_glc =>NULL() !< frozen glc runoff via rof [kg/m2/s] real, pointer, dimension(:,:) :: u_flux =>NULL() !< i-direction wind stress [Pa] real, pointer, dimension(:,:) :: v_flux =>NULL() !< j-direction wind stress [Pa] real, pointer, dimension(:,:) :: t_flux =>NULL() !< sensible heat flux [W/m2] @@ -183,6 +185,8 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: mass_berg =>NULL() !< mass of icebergs(kg/m2) real, pointer, dimension(:,:) :: hrofl =>NULL() !< heat content from liquid runoff [W/m2] real, pointer, dimension(:,:) :: hrofi =>NULL() !< heat content from frozen runoff [W/m2] + real, pointer, dimension(:,:) :: hrofl_glc =>NULL() !< heat content from liquid glc runoff [W/m2] + real, pointer, dimension(:,:) :: hrofi_glc =>NULL() !< heat content from frozen glc runoff [W/m2] real, pointer, dimension(:,:) :: hrain =>NULL() !< heat content from liquid precipitation [W/m2] real, pointer, dimension(:,:) :: hsnow =>NULL() !< heat content from frozen precipitation [W/m2] real, pointer, dimension(:,:) :: hevap =>NULL() !< heat content from evaporation [W/m2] @@ -321,7 +325,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, press=.true., fix_accum_bug=.not.CS%ustar_gustless_bug, & cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, & tau_mag=.true., ice_ncat=IOB%ice_ncat) - !call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed) @@ -494,6 +498,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j) endif + ! add liquid glc runoff flux via rof + if (associated(IOB%lrunoff_glc)) then + fluxes%lrunoff_glc(i,j) = kg_m2_s_conversion * IOB%lrunoff_glc(i-i0,j-j0) * G%mask2dT(i,j) + endif + + ! ice glc runoff flux via rof + if (associated(IOB%frunoff_glc)) then + fluxes%frunoff_glc(i,j) = kg_m2_s_conversion * IOB%frunoff_glc(i-i0,j-j0) * G%mask2dT(i,j) + endif + if (associated(IOB%ustar_berg)) & fluxes%ustar_berg(i,j) = US%m_to_Z*US%T_to_s * IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) @@ -531,6 +545,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * & IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion endif + ! notice minus sign since frunoff_glc is positive into the ocean + if (associated(IOB%frunoff_glc)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - & + IOB%frunoff_glc(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion + fluxes%latent_frunoff_glc_diag(i,j) = fluxes%latent_frunoff_glc_diag(i,j) - G%mask2dT(i,j) * & + IOB%frunoff_glc(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion + endif if (associated(IOB%q_flux)) then fluxes%latent(i,j) = fluxes%latent(i,j) + & IOB%q_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_vapor @@ -572,6 +593,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%hcond)) & fluxes%heat_content_cond(i,j) = US%W_m2_to_QRZ_T * IOB%hcond(i-i0,j-j0) * G%mask2dT(i,j) + + if (associated(IOB%hrofl_glc)) & + fluxes%heat_content_lrunoff_glc(i,j) = US%W_m2_to_QRZ_T * IOB%hrofl_glc(i-i0,j-j0) * G%mask2dT(i,j) + + if (associated(IOB%hrofi_glc)) & + fluxes%heat_content_frunoff_glc(i,j) = US%W_m2_to_QRZ_T * IOB%hrofi_glc(i-i0,j-j0) * G%mask2dT(i,j) endif ! sea ice fraction [nondim] @@ -633,7 +660,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, do j=js,je ; do i=is,ie net_FW(i,j) = US%RZ_T_to_kg_m2s * & (((fluxes%lprec(i,j) + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + & - (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + & + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) + & + fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) + & (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * US%L_to_m**2*G%areaT(i,j) net_FW2(i,j) = net_FW(i,j) / (US%L_to_m**2*G%areaT(i,j)) enddo ; enddo @@ -734,7 +762,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed) call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed) - !call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed) + call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed) if (CS%rigid_sea_ice) then call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) @@ -895,7 +923,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) - !forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j)) + forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j)) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) else ! C-grid wind stresses. @@ -1133,7 +1161,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, ! Local variables real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1]. type(directories) :: dirs - logical :: new_sim, iceberg_flux_diags, fix_ustar_gustless_bug + logical :: new_sim, iceberg_flux_diags, glc_runoff_diags, fix_ustar_gustless_bug logical :: test_value ! This is used to determine whether a logical parameter is being set explicitly. logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly. type(time_type) :: Time_frc @@ -1431,8 +1459,13 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "If true, makes available diagnostics of fluxes from icebergs "//& "as seen by MOM6.", default=.false.) + call get_param(param_file, mdl, "ALLOW_GLC_RUNOFF_DIAGNOSTICS", glc_runoff_diags, & + "If true, makes available diagnostics of separate glacier runoff fluxes"//& + "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_waves=use_waves, use_cfcs=CS%use_CFC) + use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves, & + use_cfcs=CS%use_CFC, use_glc_runoff=glc_runoff_diags) call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & "If true, allows flux adjustments to specified via the "//& @@ -1541,6 +1574,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) chks = field_chksum( iobt%fprec ) ; if (root) write(outunit,100) 'iobt%fprec ', chks 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%lrunoff_glc ) ; if (root) write(outunit,100) 'iobt%lrunoff_glc ', chks + chks = field_chksum( iobt%frunoff_glc ) ; if (root) write(outunit,100) 'iobt%frunoff_glc ', 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 @@ -1631,6 +1666,12 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) if (associated(iobt%hcond)) then chks = field_chksum( iobt%hcond ) ; if (root) write(outunit,100) 'iobt%hcond ', chks endif + if (associated(iobt%hrofl_glc)) then + chks = field_chksum( iobt%hrofl_glc ) ; if (root) write(outunit,100) 'iobt%hrofl_glc ', chks + endif + if (associated(iobt%hrofl_glc)) then + chks = field_chksum( iobt%hrofl_glc ) ; if (root) write(outunit,100) 'iobt%hrofl_glc ', chks + endif 100 FORMAT(" CHECKSUM::",A20," = ",Z20) 110 FORMAT(" CHECKSUM::",A30," = ",Z20) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 52aac958e0..3ec78bac83 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 52aac958e05cdb2471dc73f9ef7fb4e816c550f2 +Subproject commit 3ec78bac8306ef2e61a33e0c4beafa0875a2c787 diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 8faec6c495..58bb35d5a7 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -200,7 +200,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m character(len=80) :: string, string2, varName ! Temporary strings character(len=40) :: coord_units, coord_res_param ! Temporary strings character(len=MAX_PARAM_LENGTH) :: param_name - character(len=200) :: inputdir, fileName + character(len=200) :: inputdir, fileName, longString character(len=320) :: message ! Temporary strings character(len=12) :: expected_units, alt_units ! Temporary strings logical :: tmpLogical, do_sum, main_parameters @@ -680,7 +680,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! Optionally specify maximum thicknesses for each layer, enforced by moving ! the interface below a layer downward. - call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, & + call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", longString, & "Determines how to specify the maximum layer thicknesses.\n"//& "Valid options are:\n"//& " NONE - there are no maximum layer thicknesses\n"//& @@ -692,26 +692,26 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m default='NONE') message = "The list of maximum thickness for each layer." allocate(h_max(ke)) - if ( trim(string) == "NONE") then + if ( trim(longString) == "NONE") then ! Do nothing. - elseif ( trim(string) == "PARAM") then + elseif ( trim(longString) == "PARAM") then call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, & trim(message), units="m", fail_if_missing=.true., scale=GV%m_to_H) call set_regrid_max_thickness(CS, h_max) - elseif (index(trim(string),'FILE:')==1) then - if (string(6:6)=='.' .or. string(6:6)=='/') then + elseif (index(trim(longString),'FILE:')==1) then + if (longString(6:6)=='.' .or. longString(6:6)=='/') then ! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path - fileName = trim( extractWord(trim(string(6:80)), 1) ) + fileName = trim( extractWord(trim(longString(6:200)), 1) ) else ! Otherwise assume we should look for the file in INPUTDIR - fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) + fileName = trim(inputdir) // trim( extractWord(trim(longString(6:200)), 1) ) endif if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(longString)//")") - varName = trim( extractWord(trim(string(6:)), 2) ) + varName = trim( extractWord(trim(longString(6:)), 2) ) if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(longString)//")") if (len_trim(varName)==0) then if (field_exists(fileName,'h_max')) then; varName = 'h_max' elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' @@ -723,14 +723,14 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, & trim(message), units=coordinateUnits(coord_mode)) call set_regrid_max_thickness(CS, h_max, GV%m_to_H) - elseif (index(trim(string),'FNC1:')==1) then - call dz_function1( trim(string(6:)), h_max ) + elseif (index(trim(longString),'FNC1:')==1) then + call dz_function1( trim(longString(6:)), h_max ) call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, & trim(message), units=coordinateUnits(coord_mode)) call set_regrid_max_thickness(CS, h_max, GV%m_to_H) else call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string)) + "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longString)) endif deallocate(h_max) endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2cbbf6bcc7..1c94cf18fb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -78,6 +78,7 @@ module MOM use MOM_dynamics_split_RK2, only : step_MOM_dyn_split_RK2, register_restarts_dyn_split_RK2 use MOM_dynamics_split_RK2, only : initialize_dyn_split_RK2, end_dyn_split_RK2 use MOM_dynamics_split_RK2, only : MOM_dyn_split_RK2_CS, remap_dyn_split_rk2_aux_vars +use MOM_dynamics_split_RK2, only : init_dyn_split_RK2_diabatic use MOM_dynamics_split_RK2b, only : step_MOM_dyn_split_RK2b, register_restarts_dyn_split_RK2b use MOM_dynamics_split_RK2b, only : initialize_dyn_split_RK2b, end_dyn_split_RK2b use MOM_dynamics_split_RK2b, only : MOM_dyn_split_RK2b_CS, remap_dyn_split_RK2b_aux_vars @@ -2102,6 +2103,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & logical :: symmetric ! If true, use symmetric memory allocation. logical :: save_IC ! If true, save the initial conditions. logical :: do_unit_tests ! If true, call unit tests. + logical :: fpmix ! Needed to decide if BLD should be passed to RK2. logical :: test_grid_copy = .false. logical :: bulkmixedlayer ! If true, a refined bulk mixed layer scheme is used @@ -2206,6 +2208,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & default=.false.) endif + ! FPMIX is needed to decide if boundary layer depth should be passed to RK2 + call get_param(param_file, '', "FPMIX", fpmix, & + "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", & + default=.false., do_not_log=.true.) + + if (fpmix .and. .not. CS%split) then + call MOM_error(FATAL, "initialize_MOM: "//& + "FPMIX=True only works when SPLIT=True.") + endif + call get_param(param_file, "MOM", "BOUSSINESQ", Boussinesq, & "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.) call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_Boussinesq, & @@ -3334,6 +3346,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp) endif + ! GMM, the following is needed to get BLDs into the dynamics module + if (CS%split .and. fpmix) then + call init_dyn_split_RK2_diabatic(CS%diabatic_CSp, CS%dyn_split_RK2_CSp) + endif + if (associated(CS%sponge_CSp)) & call init_sponge_diags(Time, G, GV, US, diag, CS%sponge_CSp) @@ -3609,7 +3626,7 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) ! hML is needed when using the ice shelf module call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., & do_not_log=.true.) - if (use_ice_shelf) then + if (use_ice_shelf .and. associated(CS%Hml)) then call register_restart_field(CS%Hml, "hML", .false., restart_CSp, & "Mixed layer thickness", "m", conversion=US%Z_to_m) endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 35e9c2722e..217ec42c20 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -12,6 +12,7 @@ module MOM_dynamics_split_RK2 use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member use MOM_diag_mediator, only : diag_mediator_init, enable_averages use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : post_product_u, post_product_sum_u @@ -45,7 +46,9 @@ module MOM_dynamics_split_RK2 use MOM_continuity, only : continuity_init, continuity_stencil use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_CS use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_debugging, only : check_redundant +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil @@ -137,14 +140,16 @@ module MOM_dynamics_split_RK2 real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure !! anomaly in each layer due to free surface height !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure - real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean - !! to the seafloor [R L Z T-2 ~> Pa] - real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean - !! to the seafloor [R L Z T-2 ~> Pa] - type(BT_cont_type), pointer :: BT_cont => NULL() !< A structure with elements that describe the - !! effective summed open face areas as a function - !! of barotropic flow. + real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean + !! to the seafloor [R L Z T-2 ~> Pa] + real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean + !! to the seafloor [R L Z T-2 ~> Pa] + type(BT_cont_type), pointer :: BT_cont => NULL() !< A structure with elements that describe the + !! effective summed open face areas as a function + !! of barotropic flow. ! This is to allow the previous, velocity-based coupling with between the ! baroclinic and barotropic modes. @@ -174,13 +179,13 @@ module MOM_dynamics_split_RK2 !! the extent to which the treatment of gravity waves !! is forward-backward (0) or simulated backward !! Euler (1) [nondim]. 0 is often used. - logical :: debug !< If true, write verbose checksums for debugging purposes. + real :: Cemp_NL !< Empirical coefficient of non-local momentum mixing [nondim] + logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do debugging calls for open boundary conditions. - logical :: fpmix = .false. !< If true, applies profiles of momentum flux magnitude and direction. + logical :: fpmix !< If true, add non-local momentum flux increments and diffuse down the Eulerian gradient. logical :: module_is_initialized = .false. !< Record whether this module has been initialized. !>@{ Diagnostic IDs - integer :: id_uold = -1, id_vold = -1 integer :: id_uh = -1, id_vh = -1 integer :: id_umo = -1, id_vmo = -1 integer :: id_umo_2d = -1, id_vmo_2d = -1 @@ -267,6 +272,7 @@ module MOM_dynamics_split_RK2 public register_restarts_dyn_split_RK2 public initialize_dyn_split_RK2 public remap_dyn_split_RK2_aux_vars +public init_dyn_split_RK2_diabatic public end_dyn_split_RK2 !>@{ CPU time clock IDs @@ -394,7 +400,8 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f logical :: Use_Stokes_PGF ! If true, add Stokes PGF to hydrostatic PGF !---For group halo pass logical :: showCallTree, sym - + logical :: lFPpost ! Used to only post diagnostics in vertFPmix when fpmix=true and + ! in the corrector step (not the predict) integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: cont_stencil, obc_stencil, vel_stencil @@ -710,16 +717,22 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, G, GV, US, CS%vertvisc_CSp, & CS%OBC, VarMix) - call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, & - GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) if (CS%fpmix) then hbl(:,:) = 0.0 - if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:) - call vertFPmix(up, vp, uold, vold, hbl, h, forces, & - dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) - call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & - GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) + if (ASSOCIATED(CS%energetic_PBL_CSp)) & + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) + + ! lFPpost must be false in the predictor step to avoid averaging into the diagnostics + lFPpost = .false. + call vertFPmix(up, vp, uold, vold, hbl, h, forces, dt_pred, lFPpost, CS%Cemp_NL, & + G, GV, US, CS%vertvisc_CSp, CS%OBC, waves=waves) + call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, & + GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves) + else + call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, & + GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) endif if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)") @@ -960,12 +973,15 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix) - call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & - CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves) if (CS%fpmix) then - call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, & - G, GV, US, CS%vertvisc_CSp, CS%OBC) + lFPpost = .true. + call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, lFPpost, CS%Cemp_NL, & + G, GV, US, CS%vertvisc_CSp, CS%OBC, Waves=Waves) + call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & + CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves) + + else call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) endif @@ -1052,11 +1068,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f CS%CAu_pred_stored = .false. endif - if (CS%fpmix) then - if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag) - if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag) - endif - ! The time-averaged free surface height has already been set by the last call to btstep. ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH @@ -1290,6 +1301,17 @@ subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_ end subroutine remap_dyn_split_RK2_aux_vars +!> Initializes aspects of the dyn_split_RK2 that depend on diabatic processes. +!! Needed when BLDs are used in the dynamics. +subroutine init_dyn_split_RK2_diabatic(diabatic_CSp, CS) + type(diabatic_CS), intent(in) :: diabatic_CSp !< diabatic structure + type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure + + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + +end subroutine init_dyn_split_RK2_diabatic + !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & @@ -1402,8 +1424,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p "timestep for use in the predictor step of the next split RK2 timestep.", & default=.true.) call get_param(param_file, mdl, "FPMIX", CS%fpmix, & - "If true, apply profiles of momentum flux magnitude and "//& - " direction", default=.false.) + "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", & + default=.false.) + if (CS%fpmix) then + call get_param(param_file, "MOM", "CEMP_NL", CS%Cemp_NL, & + "Empirical coefficient of non-local momentum mixing.", & + units="nondim", default=3.6) + endif call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, & "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//& "variables that are needed to reproduce across restarts, similarly to "//& @@ -1480,7 +1507,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) - call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, ntrunc, CS%vertvisc_CSp) + call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, & + ntrunc, CS%vertvisc_CSp, CS%fpmix) CS%set_visc_CSp => set_visc call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) ) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index adf9d29ea3..d25df710ce 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -81,7 +81,7 @@ module MOM_forcing_type ! surface stress components and turbulent velocity scale real, pointer, dimension(:,:) :: & - !omega_w2x => NULL(), & !< the counter-clockwise angle of the wind stress with respect + omega_w2x => NULL(), & !< the counter-clockwise angle of the wind stress with respect ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. tau_mag => NULL(), & !< Magnitude of the wind stress averaged over tracer cells, !! including any contributions from sub-gridscale variability @@ -114,9 +114,10 @@ module MOM_forcing_type ! components of latent heat fluxes used for diagnostic purposes real, pointer, dimension(:,:) :: & - latent_evap_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from evaporating liquid water (typically < 0) - latent_fprec_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting fprec (typically < 0) - latent_frunoff_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from melting frunoff (calving) (typically < 0) + latent_evap_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from evaporating liquid water (typically < 0) + latent_fprec_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting fprec (typically < 0) + latent_frunoff_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting frunoff (calving) (typically < 0) + latent_frunoff_glc_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from melting glacier frunoff (typically < 0) ! water mass fluxes into the ocean [R Z T-1 ~> kg m-2 s-1]; these fluxes impact the ocean mass real, pointer, dimension(:,:) :: & @@ -126,6 +127,8 @@ module MOM_forcing_type vprec => NULL(), & !< virtual liquid precip associated w/ SSS restoring [R Z T-1 ~> kg m-2 s-1] lrunoff => NULL(), & !< liquid river runoff entering ocean [R Z T-1 ~> kg m-2 s-1] frunoff => NULL(), & !< frozen river runoff (calving) entering ocean [R Z T-1 ~> kg m-2 s-1] + lrunoff_glc => NULL(), & !< liquid river glacier runoff entering ocean [R Z T-1 ~> kg m-2 s-1] + frunoff_glc => NULL(), & !< frozen river glacier runoff entering ocean [R Z T-1 ~> kg m-2 s-1] seaice_melt => NULL() !< snow/seaice melt (positive) or formation (negative) [R Z T-1 ~> kg m-2 s-1] ! Integrated water mass fluxes into the ocean, used for passive tracer sources [H ~> m or kg m-2] @@ -137,15 +140,17 @@ module MOM_forcing_type ! heat associated with water crossing ocean surface real, pointer, dimension(:,:) :: & - heat_content_cond => NULL(), & !< heat content associated with condensating water [Q R Z T-1 ~> W m-2] - heat_content_evap => NULL(), & !< heat content associated with evaporating water [Q R Z T-1 ~> W m-2] - heat_content_lprec => NULL(), & !< heat content associated with liquid >0 precip [Q R Z T-1 ~> W m-2] - heat_content_fprec => NULL(), & !< heat content associated with frozen precip [Q R Z T-1 ~> W m-2] - heat_content_vprec => NULL(), & !< heat content associated with virtual >0 precip [Q R Z T-1 ~> W m-2] - heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff [Q R Z T-1 ~> W m-2] - heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff [Q R Z T-1 ~> W m-2] - heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean [Q R Z T-1 ~> W m-2] - heat_content_massin => NULL() !< heat content associated with mass entering ocean [Q R Z T-1 ~> W m-2] + heat_content_cond => NULL(), & !< heat content associated with condensating water [Q R Z T-1 ~> W m-2] + heat_content_evap => NULL(), & !< heat content associated with evaporating water [Q R Z T-1 ~> W m-2] + heat_content_lprec => NULL(), & !< heat content associated with liquid >0 precip [Q R Z T-1 ~> W m-2] + heat_content_fprec => NULL(), & !< heat content associated with frozen precip [Q R Z T-1 ~> W m-2] + heat_content_vprec => NULL(), & !< heat content associated with virtual >0 precip [Q R Z T-1 ~> W m-2] + heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff [Q R Z T-1 ~> W m-2] + heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff [Q R Z T-1 ~> W m-2] + heat_content_lrunoff_glc => NULL(), & !< heat content associated with liquid runoff [Q R Z T-1 ~> W m-2] + heat_content_frunoff_glc => NULL(), & !< heat content associated with frozen runoff [Q R Z T-1 ~> W m-2] + heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean [Q R Z T-1 ~> W m-2] + heat_content_massin => NULL() !< heat content associated with mass entering ocean [Q R Z T-1 ~> W m-2] ! salt mass flux (contributes to ocean mass only if non-Bouss ) real, pointer, dimension(:,:) :: & @@ -258,8 +263,8 @@ module MOM_forcing_type tau_mag => NULL(), & !< Magnitude of the wind stress averaged over tracer cells, including any !! contributions from sub-gridscale variability or gustiness [R L Z T-2 ~> Pa] ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. - net_mass_src => NULL() !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1] - !omega_w2x => NULL() !< the counter-clockwise angle of the wind stress with respect + net_mass_src => NULL(), & !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1] + omega_w2x => NULL() !< the counter-clockwise angle of the wind stress with respect !! to the horizontal abscissa (x-coordinate) at tracer points [rad]. ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice) @@ -323,6 +328,7 @@ module MOM_forcing_type integer :: id_precip = -1, id_vprec = -1 integer :: id_lprec = -1, id_fprec = -1 integer :: id_lrunoff = -1, id_frunoff = -1 + integer :: id_lrunoff_glc = -1, id_frunoff_glc = -1 integer :: id_net_massout = -1, id_net_massin = -1 integer :: id_massout_flux = -1, id_massin_flux = -1 integer :: id_seaice_melt = -1 @@ -332,6 +338,7 @@ module MOM_forcing_type integer :: id_total_precip = -1, id_total_vprec = -1 integer :: id_total_lprec = -1, id_total_fprec = -1 integer :: id_total_lrunoff = -1, id_total_frunoff = -1 + integer :: id_total_lrunoff_glc = -1, id_total_frunoff_glc = -1 integer :: id_total_net_massout = -1, id_total_net_massin = -1 integer :: id_total_seaice_melt = -1 @@ -341,34 +348,38 @@ module MOM_forcing_type integer :: id_precip_ga = -1, id_vprec_ga= -1 ! heat flux diagnostic handles - integer :: id_net_heat_coupler = -1, id_net_heat_surface = -1 - integer :: id_sens = -1, id_LwLatSens = -1 - integer :: id_sw = -1, id_lw = -1 - integer :: id_sw_vis = -1, id_sw_nir = -1 - integer :: id_lat_evap = -1, id_lat_frunoff = -1 - integer :: id_lat = -1, id_lat_fprec = -1 - integer :: id_heat_content_lrunoff= -1, id_heat_content_frunoff = -1 - integer :: id_heat_content_lprec = -1, id_heat_content_fprec = -1 - integer :: id_heat_content_cond = -1, id_heat_content_surfwater= -1 - integer :: id_heat_content_evap = -1 - integer :: id_heat_content_vprec = -1, id_heat_content_massout = -1 - integer :: id_heat_added = -1, id_heat_content_massin = -1 - integer :: id_hfrainds = -1, id_hfrunoffds = -1 - integer :: id_seaice_melt_heat = -1 + integer :: id_net_heat_coupler = -1, id_net_heat_surface = -1 + integer :: id_sens = -1, id_LwLatSens = -1 + integer :: id_sw = -1, id_lw = -1 + integer :: id_sw_vis = -1, id_sw_nir = -1 + integer :: id_lat_evap = -1, id_lat_frunoff = -1 + integer :: id_lat_frunoff_glc = -1 + integer :: id_lat = -1, id_lat_fprec = -1 + integer :: id_heat_content_lrunoff = -1, id_heat_content_frunoff = -1 + integer :: id_heat_content_lrunoff_glc= -1, id_heat_content_frunoff_glc= -1 + integer :: id_heat_content_lprec = -1, id_heat_content_fprec = -1 + integer :: id_heat_content_cond = -1, id_heat_content_surfwater = -1 + integer :: id_heat_content_evap = -1 + integer :: id_heat_content_vprec = -1, id_heat_content_massout = -1 + integer :: id_heat_added = -1, id_heat_content_massin = -1 + integer :: id_hfrainds = -1, id_hfrunoffds = -1 + integer :: id_seaice_melt_heat = -1 ! global area integrated heat flux diagnostic handles - integer :: id_total_net_heat_coupler = -1, id_total_net_heat_surface = -1 - integer :: id_total_sens = -1, id_total_LwLatSens = -1 - integer :: id_total_sw = -1, id_total_lw = -1 - integer :: id_total_lat_evap = -1, id_total_lat_frunoff = -1 - integer :: id_total_lat = -1, id_total_lat_fprec = -1 - integer :: id_total_heat_content_lrunoff= -1, id_total_heat_content_frunoff = -1 - integer :: id_total_heat_content_lprec = -1, id_total_heat_content_fprec = -1 - integer :: id_total_heat_content_cond = -1, id_total_heat_content_surfwater= -1 - integer :: id_total_heat_content_evap = -1 - integer :: id_total_heat_content_vprec = -1, id_total_heat_content_massout = -1 - integer :: id_total_heat_added = -1, id_total_heat_content_massin = -1 - integer :: id_total_seaice_melt_heat = -1 + integer :: id_total_net_heat_coupler = -1, id_total_net_heat_surface = -1 + integer :: id_total_sens = -1, id_total_LwLatSens = -1 + integer :: id_total_sw = -1, id_total_lw = -1 + integer :: id_total_lat_evap = -1, id_total_lat_frunoff = -1 + integer :: id_total_lat_frunoff_glc = -1 + integer :: id_total_lat = -1, id_total_lat_fprec = -1 + integer :: id_total_heat_content_lrunoff = -1, id_total_heat_content_frunoff = -1 + integer :: id_total_heat_content_lrunoff_glc= -1, id_total_heat_content_frunoff_glc=-1 + integer :: id_total_heat_content_lprec = -1, id_total_heat_content_fprec = -1 + integer :: id_total_heat_content_cond = -1, id_total_heat_content_surfwater = -1 + integer :: id_total_heat_content_evap = -1 + integer :: id_total_heat_content_vprec = -1, id_total_heat_content_massout = -1 + integer :: id_total_heat_added = -1, id_total_heat_content_massin = -1 + integer :: id_total_seaice_melt_heat = -1 ! global area averaged heat flux diagnostic handles integer :: id_net_heat_coupler_ga = -1, id_net_heat_surface_ga = -1 @@ -397,7 +408,7 @@ module MOM_forcing_type integer :: id_taux = -1 integer :: id_tauy = -1 integer :: id_ustar = -1 - !integer :: id_omega_w2x = -1 + integer :: id_omega_w2x = -1 integer :: id_tau_mag = -1 integer :: id_psurf = -1 integer :: id_TKE_tidal = -1 @@ -609,23 +620,27 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & ! net volume/mass of liquid and solid passing through surface boundary fluxes netMassInOut(i) = dt * (scale * & - (((((( fluxes%lprec(i,j) & + (((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & + + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & - + fluxes%frunoff(i,j) )) + + fluxes%frunoff(i,j) ) & + + fluxes%frunoff_glc(i,j))) if (do_NMIOr) then ! Repeat the above code without multiplying by a timestep for legacy reasons netMassInOut_rate(i) = (scale * & - (((((( fluxes%lprec(i,j) & + (((((((( fluxes%lprec(i,j) & + fluxes%fprec(i,j) ) & + fluxes%evap(i,j) ) & + fluxes%lrunoff(i,j) ) & + + fluxes%lrunoff_glc(i,j)) & + fluxes%vprec(i,j) ) & + fluxes%seaice_melt(i,j)) & - + fluxes%frunoff(i,j) )) + + fluxes%frunoff(i,j) ) & + + fluxes%frunoff_glc(i,j))) endif ! smg: @@ -700,6 +715,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & ! remove lrunoff*SST here, to counteract its addition elsewhere net_heat(i) = (net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_lrunoff(i,j)) - & (GV%RZ_to_H * (scale * dt)) * fluxes%lrunoff(i,j) * T(i,1) + net_heat(i) = (net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_lrunoff_glc(i,j)) - & + (GV%RZ_to_H * (scale * dt)) * fluxes%lrunoff_glc(i,j) * T(i,1) !BGR-Jul 5, 2017{ !Intentionally neglect the following contribution to rate for legacy reasons. !if (do_NHR) net_heat_rate(i) = (net_heat_rate(i) + (scale*I_Cp_Hconvert) * fluxes%heat_content_lrunoff(i,j)) - & @@ -708,6 +725,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (calculate_diags .and. associated(tv%TempxPmE)) then tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & (I_Cp*fluxes%heat_content_lrunoff(i,j) - fluxes%lrunoff(i,j)*T(i,1)) + tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & + (I_Cp*fluxes%heat_content_lrunoff_glc(i,j) - fluxes%lrunoff_glc(i,j)*T(i,1)) endif endif @@ -717,6 +736,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & ! remove frunoff*SST here, to counteract its addition elsewhere net_heat(i) = net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_frunoff(i,j) - & (GV%RZ_to_H * (scale * dt)) * fluxes%frunoff(i,j) * T(i,1) + net_heat(i) = net_heat(i) + (scale*(dt * I_Cp_Hconvert)) * fluxes%heat_content_frunoff_glc(i,j) - & + (GV%RZ_to_H * (scale * dt)) * fluxes%frunoff_glc(i,j) * T(i,1) !BGR-Jul 5, 2017{ !Intentionally neglect the following contribution to rate for legacy reasons. ! if (do_NHR) net_heat_rate(i) = net_heat_rate(i) + (scale*I_Cp_Hconvert) * fluxes%heat_content_frunoff(i,j) - & @@ -725,6 +746,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (calculate_diags .and. associated(tv%TempxPmE)) then tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & (I_Cp*fluxes%heat_content_frunoff(i,j) - fluxes%frunoff(i,j)*T(i,1)) + tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * & + (I_Cp*fluxes%heat_content_frunoff_glc(i,j) - fluxes%frunoff_glc(i,j)*T(i,1)) endif endif @@ -748,6 +771,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (.not. do_enthalpy) then net_heat(i) = net_heat(i) + (scale * dt * I_Cp_Hconvert * & (fluxes%heat_content_lrunoff(i,j) + fluxes%heat_content_frunoff(i,j) + & + fluxes%heat_content_lrunoff_glc(i,j) + fluxes%heat_content_frunoff_glc(i,j) + & fluxes%heat_content_lprec(i,j) + fluxes%heat_content_fprec(i,j) + & fluxes%heat_content_evap(i,j) + fluxes%heat_content_cond(i,j))) endif @@ -876,6 +900,9 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (associated(fluxes%lrunoff) .and. associated(fluxes%heat_content_lrunoff)) then fluxes%heat_content_lrunoff(i,j) = tv%C_p*fluxes%lrunoff(i,j)*T(i,1) endif + if (associated(fluxes%lrunoff_glc) .and. associated(fluxes%heat_content_lrunoff_glc)) then + fluxes%heat_content_lrunoff_glc(i,j) = tv%C_p*fluxes%lrunoff_glc(i,j)*T(i,1) + endif endif ! Icebergs enter ocean at SST if land model does not provide calving heat content. @@ -883,6 +910,9 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & if (associated(fluxes%frunoff) .and. associated(fluxes%heat_content_frunoff)) then fluxes%heat_content_frunoff(i,j) = tv%C_p*fluxes%frunoff(i,j)*T(i,1) endif + if (associated(fluxes%frunoff_glc) .and. associated(fluxes%heat_content_frunoff_glc)) then + fluxes%heat_content_frunoff_glc(i,j) = tv%C_p*fluxes%frunoff_glc(i,j)*T(i,1) + endif endif elseif (.not. do_enthalpy) then @@ -905,6 +935,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & fluxes%heat_content_fprec(i,j) + & fluxes%heat_content_lrunoff(i,j) + & fluxes%heat_content_frunoff(i,j) + & + fluxes%heat_content_lrunoff_glc(i,j) + & + fluxes%heat_content_frunoff_glc(i,j) + & fluxes%heat_content_evap(i,j) + & fluxes%heat_content_cond(i,j)) endif @@ -1309,6 +1341,9 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) if (associated(fluxes%latent_frunoff_diag)) & call hchksum(fluxes%latent_frunoff_diag, mesg//" fluxes%latent_frunoff_diag", G%HI, & haloshift=hshift, scale=US%QRZ_T_to_W_m2) + if (associated(fluxes%latent_frunoff_glc_diag)) & + call hchksum(fluxes%latent_frunoff_glc_diag, mesg//" fluxes%latent_frunoff_glc_diag", G%HI, & + haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%sens)) & call hchksum(fluxes%sens, mesg//" fluxes%sens", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%evap)) & @@ -1338,14 +1373,24 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%ustar_tidal, mesg//" fluxes%ustar_tidal", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) if (associated(fluxes%lrunoff)) & call hchksum(fluxes%lrunoff, mesg//" fluxes%lrunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + if (associated(fluxes%lrunoff_glc)) & + call hchksum(fluxes%lrunoff_glc, mesg//" fluxes%lrunoff_glc", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%frunoff)) & call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + if (associated(fluxes%frunoff_glc)) & + call hchksum(fluxes%frunoff_glc, mesg//" fluxes%frunoff_glc", 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%QRZ_T_to_W_m2) + if (associated(fluxes%heat_content_lrunoff_glc)) & + call hchksum(fluxes%heat_content_lrunoff_glc, mesg//" fluxes%heat_content_lrunoff_glc", G%HI, & + 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) + if (associated(fluxes%heat_content_frunoff_glc)) & + call hchksum(fluxes%heat_content_frunoff_glc, mesg//" fluxes%heat_content_frunoff_glc", G%HI, & + haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_lprec)) & call hchksum(fluxes%heat_content_lprec, mesg//" fluxes%heat_content_lprec", G%HI, & haloshift=hshift, scale=US%QRZ_T_to_W_m2) @@ -1448,6 +1493,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) call locMsg(fluxes%latent_evap_diag,'latent_evap_diag') call locMsg(fluxes%latent_fprec_diag,'latent_fprec_diag') call locMsg(fluxes%latent_frunoff_diag,'latent_frunoff_diag') + call locMsg(fluxes%latent_frunoff_glc_diag,'latent_frunoff_glc_diag') call locMsg(fluxes%sens,'sens') call locMsg(fluxes%evap,'evap') call locMsg(fluxes%lprec,'lprec') @@ -1460,9 +1506,13 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) call locMsg(fluxes%TKE_tidal,'TKE_tidal') call locMsg(fluxes%ustar_tidal,'ustar_tidal') call locMsg(fluxes%lrunoff,'lrunoff') + call locMsg(fluxes%lrunoff_glc,'lrunoff_glc') call locMsg(fluxes%frunoff,'frunoff') + call locMsg(fluxes%frunoff_glc,'frunoff_glc') call locMsg(fluxes%heat_content_lrunoff,'heat_content_lrunoff') + call locMsg(fluxes%heat_content_lrunoff_glc,'heat_content_lrunoff_glc') call locMsg(fluxes%heat_content_frunoff,'heat_content_frunoff') + call locMsg(fluxes%heat_content_frunoff_glc,'heat_content_frunoff_glc') call locMsg(fluxes%heat_content_lprec,'heat_content_lprec') call locMsg(fluxes%heat_content_fprec,'heat_content_fprec') call locMsg(fluxes%heat_content_vprec,'heat_content_vprec') @@ -1489,7 +1539,8 @@ end subroutine forcing_SinglePointPrint !> Register members of the forcing type for diagnostics -subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_waves, use_cfcs) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_waves, & + use_cfcs, use_glc_runoff) 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 @@ -1498,6 +1549,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 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 + logical, optional, intent(in) :: use_glc_runoff !< If true, allow separate glacial runoff diagnostics ! Clock for forcing diagnostics handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1525,8 +1577,8 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', & 'm s-1', conversion=US%Z_to_m*US%s_to_T) - !handles%id_omega_w2x = register_diag_field('ocean_model', 'omega_w2x', diag%axesT1, Time, & - ! 'Counter-clockwise angle of the wind stress from the horizontal axis.', 'rad') + handles%id_omega_w2x = register_diag_field('ocean_model', 'omega_w2x', diag%axesT1, Time, & + 'Counter-clockwise angle of the wind stress from the horizontal axis.', 'rad') if (present(use_berg_fluxes)) then if (use_berg_fluxes) then @@ -1634,6 +1686,18 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_standard_name='water_flux_into_sea_water_from_rivers', & cmor_long_name='Water Flux into Sea Water From Rivers') + if (present(use_glc_runoff)) then + handles%id_frunoff_glc = register_diag_field('ocean_model', 'frunoff_glc', diag%axesT1, Time, & + 'Frozen glacier runoff (calving) and iceberg melt into ocean', & + units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + standard_name='glc_water_flux_into_sea_water_from_icebergs') ! todo: update cmor names + + handles%id_lrunoff_glc = register_diag_field('ocean_model', 'lrunoff_glc', diag%axesT1, Time, & + 'Liquid runoff (glaciers) into ocean', & + units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + standard_name='water_flux_into_sea_water_from_glaciers') ! todo: update cmor names + endif + handles%id_net_massout = register_diag_field('ocean_model', 'net_massout', diag%axesT1, Time, & 'Net mass leaving the ocean due to evaporation, seaice formation', & 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) @@ -1707,6 +1771,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_standard_name='water_flux_into_sea_water_from_rivers_area_integrated', & cmor_long_name='Water Flux into Sea Water From Rivers Area Integrated') + if (present(use_glc_runoff)) then + handles%id_total_frunoff_glc = register_scalar_field('ocean_model', 'total_frunoff_glc', Time, diag, & + long_name='Area integrated frozen glacier runoff (calving) & iceberg melt into ocean', units='kg s-1') + + handles%id_total_lrunoff_glc = register_scalar_field('ocean_model', 'total_lrunoff_glc', Time, diag,& + long_name='Area integrated liquid glacier runoff into ocean', units='kg s-1') + endif + handles%id_total_net_massout = register_scalar_field('ocean_model', 'total_net_massout', Time, diag, & long_name='Area integrated mass leaving ocean due to evap and seaice form', units='kg s-1') @@ -1765,6 +1837,16 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 'W m-2', conversion=US%QRZ_T_to_W_m2, & standard_name='temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water') + if (present(use_glc_runoff)) then + handles%id_heat_content_frunoff_glc = register_diag_field('ocean_model', 'heat_content_frunoff_glc', & + diag%axesT1, Time, 'Heat content (relative to 0C) of solid glacier runoff into ocean', & + 'W m-2', conversion=US%QRZ_T_to_W_m2) + + handles%id_heat_content_lrunoff_glc = register_diag_field('ocean_model', 'heat_content_lrunoff_glc', & + diag%axesT1, Time, 'Heat content (relative to 0C) of liquid glacier runoff into ocean', & + 'W m-2', conversion=US%QRZ_T_to_W_m2) + endif + handles%id_hfrunoffds = register_diag_field('ocean_model', 'hfrunoffds', & diag%axesT1, Time, 'Heat content (relative to 0C) of liquid+solid runoff into ocean', & 'W m-2', conversion=US%QRZ_T_to_W_m2, & @@ -1868,6 +1950,11 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_standard_name='heat_flux_into_sea_water_due_to_iceberg_thermodynamics', & cmor_long_name='Latent Heat to Melt Frozen Runoff/Iceberg') + if (present(use_glc_runoff)) then + handles%id_lat_frunoff_glc = register_diag_field('ocean_model', 'latent_frunoff_glc', diag%axesT1, Time, & + 'Latent heat flux into ocean due to melting of frozen glacier runoff', 'W m-2', conversion=US%QRZ_T_to_W_m2) + endif + handles%id_sens = register_diag_field('ocean_model', 'sensible', diag%axesT1, Time, & 'Sensible heat flux into ocean', 'W m-2', conversion=US%QRZ_T_to_W_m2, & standard_name='surface_downward_sensible_heat_flux', & @@ -1907,6 +1994,18 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_long_name= & 'Temperature Flux due to Runoff Expressed as Heat Flux into Sea Water Area Integrated') + if (present(use_glc_runoff)) then + handles%id_total_heat_content_frunoff_glc = register_scalar_field('ocean_model', & + 'total_heat_content_frunoff_glc', Time, diag, & + long_name='Area integrated heat content (relative to 0C) of solid glacier runoff', & + units='W') ! todo: update cmor names + + handles%id_total_heat_content_lrunoff_glc = register_scalar_field('ocean_model', & + 'total_heat_content_lrunoff_glc', Time, diag, & + long_name='Area integrated heat content (relative to 0C) of liquid glacier runoff', & + units='W') ! todo: update cmor names + endif + handles%id_total_heat_content_lprec = register_scalar_field('ocean_model', & 'total_heat_content_lprec', Time, diag, & long_name='Area integrated heat content (relative to 0C) of liquid precip', & @@ -2024,6 +2123,13 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_long_name= & 'Heat Flux into Sea Water due to Iceberg Thermodynamics Area Integrated') + if (present(use_glc_runoff)) then + handles%id_total_lat_frunoff_glc = register_scalar_field('ocean_model', & + 'total_lat_frunoff_glc', Time, diag, & + long_name='Area integrated latent heat flux due to melting frozen glacier runoff', & + units='W') ! todo: update cmor names + endif + handles%id_total_sens = register_scalar_field('ocean_model', & 'total_sens', Time, diag, & long_name='Area integrated downward sensible heat flux', & @@ -2287,6 +2393,8 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%vprec(i,j) = wt1*fluxes%vprec(i,j) + wt2*flux_tmp%vprec(i,j) fluxes%lrunoff(i,j) = wt1*fluxes%lrunoff(i,j) + wt2*flux_tmp%lrunoff(i,j) fluxes%frunoff(i,j) = wt1*fluxes%frunoff(i,j) + wt2*flux_tmp%frunoff(i,j) + fluxes%lrunoff_glc(i,j) = wt1*fluxes%lrunoff_glc(i,j) + wt2*flux_tmp%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j) = wt1*fluxes%frunoff_glc(i,j) + wt2*flux_tmp%frunoff_glc(i,j) fluxes%seaice_melt(i,j) = wt1*fluxes%seaice_melt(i,j) + wt2*flux_tmp%seaice_melt(i,j) fluxes%sw(i,j) = wt1*fluxes%sw(i,j) + wt2*flux_tmp%sw(i,j) fluxes%sw_vis_dir(i,j) = wt1*fluxes%sw_vis_dir(i,j) + wt2*flux_tmp%sw_vis_dir(i,j) @@ -2340,6 +2448,18 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%heat_content_frunoff(i,j) = wt1*fluxes%heat_content_frunoff(i,j) + wt2*flux_tmp%heat_content_frunoff(i,j) enddo ; enddo endif + if (associated(fluxes%heat_content_lrunoff_glc) .and. associated(flux_tmp%heat_content_lrunoff_glc)) then + do j=js,je ; do i=is,ie + fluxes%heat_content_lrunoff_glc(i,j) = wt1*fluxes%heat_content_lrunoff_glc(i,j) + & + wt2*flux_tmp%heat_content_lrunoff_glc(i,j) + enddo ; enddo + endif + if (associated(fluxes%heat_content_frunoff_glc) .and. associated(flux_tmp%heat_content_frunoff_glc)) then + do j=js,je ; do i=is,ie + fluxes%heat_content_frunoff_glc(i,j) = wt1*fluxes%heat_content_frunoff_glc(i,j) + & + wt2*flux_tmp%heat_content_frunoff_glc(i,j) + enddo ; enddo + endif if (associated(fluxes%ustar_shelf) .and. associated(flux_tmp%ustar_shelf)) then do i=isd,ied ; do j=jsd,jed @@ -2389,11 +2509,11 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) fluxes%ustar(i,j) = forces%ustar(i,j) enddo ; enddo endif - !if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then - ! do j=js,je ; do i=is,ie - ! fluxes%omega_w2x(i,j) = forces%omega_w2x(i,j) - ! enddo ; enddo - !endif + if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then + do j=js,je ; do i=is,ie + fluxes%omega_w2x(i,j) = forces%omega_w2x(i,j) + enddo ; enddo + endif if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then do j=js,je ; do i=is,ie fluxes%tau_mag(i,j) = forces%tau_mag(i,j) @@ -2511,6 +2631,12 @@ subroutine get_net_mass_forcing(fluxes, G, US, net_mass_src) if (associated(fluxes%frunoff)) then ; do j=js,je ; do i=is,ie net_mass_src(i,j) = net_mass_src(i,j) + fluxes%frunoff(i,j) enddo ; enddo ; endif + if (associated(fluxes%lrunoff_glc)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%lrunoff_glc(i,j) + enddo ; enddo ; endif + if (associated(fluxes%frunoff_glc)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%frunoff_glc(i,j) + enddo ; enddo ; endif if (associated(fluxes%evap)) then ; do j=js,je ; do i=is,ie net_mass_src(i,j) = net_mass_src(i,j) + fluxes%evap(i,j) enddo ; enddo ; endif @@ -2535,11 +2661,11 @@ subroutine copy_back_forcing_fields(fluxes, forces, G) forces%ustar(i,j) = fluxes%ustar(i,j) enddo ; enddo endif - !if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then - ! do j=js,je ; do i=is,ie - ! forces%omega_w2x(i,j) = fluxes%omega_w2x(i,j) - ! enddo ; enddo - !endif + if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then + do j=js,je ; do i=is,ie + forces%omega_w2x(i,j) = fluxes%omega_w2x(i,j) + enddo ; enddo + endif if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then do j=js,je ; do i=is,ie forces%tau_mag(i,j) = fluxes%tau_mag(i,j) @@ -2671,6 +2797,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%evap)) res(i,j) = res(i,j) + fluxes%evap(i,j) if (associated(fluxes%lrunoff)) res(i,j) = res(i,j) + fluxes%lrunoff(i,j) if (associated(fluxes%frunoff)) res(i,j) = res(i,j) + fluxes%frunoff(i,j) + if (associated(fluxes%lrunoff_glc)) res(i,j) = res(i,j) + fluxes%lrunoff_glc(i,j) + if (associated(fluxes%frunoff_glc)) res(i,j) = res(i,j) + fluxes%frunoff_glc(i,j) if (associated(fluxes%vprec)) res(i,j) = res(i,j) + fluxes%vprec(i,j) if (associated(fluxes%seaice_melt)) res(i,j) = res(i,j) + fluxes%seaice_melt(i,j) enddo ; enddo @@ -2718,6 +2846,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%fprec)) res(i,j) = res(i,j) + fluxes%fprec(i,j) if (associated(fluxes%lrunoff)) res(i,j) = res(i,j) + fluxes%lrunoff(i,j) if (associated(fluxes%frunoff)) res(i,j) = res(i,j) + fluxes%frunoff(i,j) + if (associated(fluxes%lrunoff_glc)) res(i,j) = res(i,j) + fluxes%lrunoff_glc(i,j) + if (associated(fluxes%frunoff_glc)) res(i,j) = res(i,j) + fluxes%frunoff_glc(i,j) if (associated(fluxes%lprec)) then if (fluxes%lprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j) @@ -2813,6 +2943,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h endif endif + if (associated(fluxes%lrunoff_glc)) then + if (handles%id_lrunoff_glc > 0) call post_data(handles%id_lrunoff_glc, fluxes%lrunoff_glc, diag) + if (handles%id_total_lrunoff_glc > 0) then + total_transport = global_area_integral(fluxes%lrunoff_glc, G, scale=US%RZ_T_to_kg_m2s) + call post_data(handles%id_total_lrunoff_glc, total_transport, diag) + endif + endif + if (associated(fluxes%frunoff)) then if (handles%id_frunoff > 0) call post_data(handles%id_frunoff, fluxes%frunoff, diag) if (handles%id_total_frunoff > 0) then @@ -2821,6 +2959,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h endif endif + if (associated(fluxes%frunoff_glc)) then + if (handles%id_frunoff_glc > 0) call post_data(handles%id_frunoff_glc, fluxes%frunoff_glc, diag) + if (handles%id_total_frunoff_glc > 0) then + total_transport = global_area_integral(fluxes%frunoff_glc, G, scale=US%RZ_T_to_kg_m2s) + call post_data(handles%id_total_frunoff_glc, total_transport, diag) + endif + endif + if (associated(fluxes%seaice_melt)) then if (handles%id_seaice_melt > 0) call post_data(handles%id_seaice_melt, fluxes%seaice_melt, diag) if (handles%id_total_seaice_melt > 0) then @@ -2838,12 +2984,26 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_total_heat_content_lrunoff, total_transport, diag) endif + + if ((handles%id_heat_content_lrunoff_glc > 0) .and. associated(fluxes%heat_content_lrunoff_glc)) & + call post_data(handles%id_heat_content_lrunoff_glc, fluxes%heat_content_lrunoff_glc, diag) + if ((handles%id_total_heat_content_lrunoff_glc > 0) .and. associated(fluxes%heat_content_lrunoff_glc)) then + total_transport = global_area_integral(fluxes%heat_content_lrunoff_glc, G, scale=US%QRZ_T_to_W_m2) + call post_data(handles%id_total_heat_content_lrunoff_glc, total_transport, diag) + endif + if ((handles%id_heat_content_frunoff > 0) .and. associated(fluxes%heat_content_frunoff)) & call post_data(handles%id_heat_content_frunoff, fluxes%heat_content_frunoff, diag) if ((handles%id_total_heat_content_frunoff > 0) .and. associated(fluxes%heat_content_frunoff)) then total_transport = global_area_integral(fluxes%heat_content_frunoff, G, scale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_frunoff, total_transport, diag) endif + if ((handles%id_heat_content_frunoff_glc > 0) .and. associated(fluxes%heat_content_frunoff_glc)) & + call post_data(handles%id_heat_content_frunoff_glc, fluxes%heat_content_frunoff_glc, diag) + if ((handles%id_total_heat_content_frunoff_glc > 0) .and. associated(fluxes%heat_content_frunoff_glc)) then + total_transport = global_area_integral(fluxes%heat_content_frunoff_glc, G, scale=US%QRZ_T_to_W_m2) + call post_data(handles%id_total_heat_content_frunoff_glc, total_transport, diag) + endif if ((handles%id_heat_content_lprec > 0) .and. associated(fluxes%heat_content_lprec)) & call post_data(handles%id_heat_content_lprec, fluxes%heat_content_lprec, diag) @@ -2929,6 +3089,10 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) if (associated(fluxes%heat_content_frunoff)) & res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j) + if (associated(fluxes%heat_content_lrunoff_glc)) & + res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j) + if (associated(fluxes%heat_content_frunoff_glc)) & + res(i,j) = res(i,j) + fluxes%heat_content_frunoff_glc(i,j) if (associated(fluxes%heat_content_lprec)) & res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j) if (associated(fluxes%heat_content_fprec)) & @@ -2961,12 +3125,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_heat_content_surfwater > 0 .or. handles%id_total_heat_content_surfwater > 0) then do j=js,je ; do i=is,ie res(i,j) = 0.0 - if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) - if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j) - if (associated(fluxes%heat_content_lprec)) res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j) - if (associated(fluxes%heat_content_fprec)) res(i,j) = res(i,j) + fluxes%heat_content_fprec(i,j) - if (associated(fluxes%heat_content_vprec)) res(i,j) = res(i,j) + fluxes%heat_content_vprec(i,j) - if (associated(fluxes%heat_content_cond)) res(i,j) = res(i,j) + fluxes%heat_content_cond(i,j) + if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) + if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j) + if (associated(fluxes%heat_content_lrunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j) + if (associated(fluxes%heat_content_frunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff_glc(i,j) + if (associated(fluxes%heat_content_lprec)) res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j) + if (associated(fluxes%heat_content_fprec)) res(i,j) = res(i,j) + fluxes%heat_content_fprec(i,j) + if (associated(fluxes%heat_content_vprec)) res(i,j) = res(i,j) + fluxes%heat_content_vprec(i,j) + if (associated(fluxes%heat_content_cond)) res(i,j) = res(i,j) + fluxes%heat_content_cond(i,j) if (mom_enthalpy) then if (associated(fluxes%heat_content_massout)) res(i,j) = res(i,j) + fluxes%heat_content_massout(i,j) else @@ -2986,6 +3152,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h res(i,j) = 0.0 if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j) + if (associated(fluxes%heat_content_lrunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff_glc(i,j) + if (associated(fluxes%heat_content_frunoff_glc)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff_glc(i,j) enddo ; enddo call post_data(handles%id_hfrunoffds, res, diag) endif @@ -3095,6 +3263,14 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_total_lat_frunoff, total_transport, diag) endif + if ((handles%id_lat_frunoff_glc > 0) .and. associated(fluxes%latent_frunoff_glc_diag)) then + call post_data(handles%id_lat_frunoff_glc, fluxes%latent_frunoff_glc_diag, diag) + endif + if (handles%id_total_lat_frunoff_glc > 0 .and. associated(fluxes%latent_frunoff_glc_diag)) then + total_transport = global_area_integral(fluxes%latent_frunoff_glc_diag, G, scale=US%QRZ_T_to_W_m2) + call post_data(handles%id_total_lat_frunoff_glc, total_transport, diag) + endif + if ((handles%id_sens > 0) .and. associated(fluxes%sens)) then call post_data(handles%id_sens, fluxes%sens, diag) endif @@ -3191,8 +3367,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) & call post_data(handles%id_ustar, fluxes%ustar, diag) - !if ((handles%id_omega_w2x > 0) .and. associated(fluxes%omega_w2x)) & - ! call post_data(handles%id_omega_w2x, fluxes%omega_w2x, diag) + if ((handles%id_omega_w2x > 0) .and. associated(fluxes%omega_w2x)) & + call post_data(handles%id_omega_w2x, fluxes%omega_w2x, diag) if ((handles%id_ustar_berg > 0) .and. associated(fluxes%ustar_berg)) & call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) @@ -3278,6 +3454,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%vprec,isd,ied,jsd,jed, water) call myAlloc(fluxes%lrunoff,isd,ied,jsd,jed, water) call myAlloc(fluxes%frunoff,isd,ied,jsd,jed, water) + call myAlloc(fluxes%lrunoff_glc,isd,ied,jsd,jed, water) + call myAlloc(fluxes%frunoff_glc,isd,ied,jsd,jed, water) call myAlloc(fluxes%seaice_melt,isd,ied,jsd,jed, water) call myAlloc(fluxes%netMassOut,isd,ied,jsd,jed, water) call myAlloc(fluxes%netMassIn,isd,ied,jsd,jed, water) @@ -3289,6 +3467,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%latent_evap_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_fprec_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%latent_frunoff_diag,isd,ied,jsd,jed, heat) + call myAlloc(fluxes%latent_frunoff_glc_diag,isd,ied,jsd,jed, heat) call myAlloc(fluxes%salt_flux,isd,ied,jsd,jed, salt) @@ -3300,6 +3479,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%heat_content_vprec,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_lrunoff,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_frunoff,isd,ied,jsd,jed, .true.) + call myAlloc(fluxes%heat_content_lrunoff_glc,isd,ied,jsd,jed, .true.) + call myAlloc(fluxes%heat_content_frunoff_glc,isd,ied,jsd,jed, .true.) call myAlloc(fluxes%heat_content_massout,isd,ied,jsd,jed, enthalpy_mom) call myAlloc(fluxes%heat_content_massin,isd,ied,jsd,jed, enthalpy_mom) endif ; endif @@ -3567,7 +3748,7 @@ end subroutine myAlloc_3d subroutine deallocate_forcing_type(fluxes) type(forcing), intent(inout) :: fluxes !< Forcing fields structure - !if (associated(fluxes%omega_w2x)) deallocate(fluxes%omega_w2x) + if (associated(fluxes%omega_w2x)) deallocate(fluxes%omega_w2x) if (associated(fluxes%ustar)) deallocate(fluxes%ustar) if (associated(fluxes%ustar_gustless)) deallocate(fluxes%ustar_gustless) if (associated(fluxes%tau_mag)) deallocate(fluxes%tau_mag) @@ -3583,10 +3764,13 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%latent_evap_diag)) deallocate(fluxes%latent_evap_diag) if (associated(fluxes%latent_fprec_diag)) deallocate(fluxes%latent_fprec_diag) if (associated(fluxes%latent_frunoff_diag)) deallocate(fluxes%latent_frunoff_diag) + if (associated(fluxes%latent_frunoff_glc_diag)) deallocate(fluxes%latent_frunoff_glc_diag) if (associated(fluxes%sens)) deallocate(fluxes%sens) if (associated(fluxes%heat_added)) deallocate(fluxes%heat_added) if (associated(fluxes%heat_content_lrunoff)) deallocate(fluxes%heat_content_lrunoff) if (associated(fluxes%heat_content_frunoff)) deallocate(fluxes%heat_content_frunoff) + if (associated(fluxes%heat_content_lrunoff_glc)) deallocate(fluxes%heat_content_lrunoff_glc) + if (associated(fluxes%heat_content_frunoff_glc)) deallocate(fluxes%heat_content_frunoff_glc) if (associated(fluxes%heat_content_lprec)) deallocate(fluxes%heat_content_lprec) if (associated(fluxes%heat_content_fprec)) deallocate(fluxes%heat_content_fprec) if (associated(fluxes%heat_content_cond)) deallocate(fluxes%heat_content_cond) @@ -3599,6 +3783,8 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%vprec)) deallocate(fluxes%vprec) if (associated(fluxes%lrunoff)) deallocate(fluxes%lrunoff) if (associated(fluxes%frunoff)) deallocate(fluxes%frunoff) + if (associated(fluxes%lrunoff_glc)) deallocate(fluxes%lrunoff_glc) + if (associated(fluxes%frunoff_glc)) deallocate(fluxes%frunoff_glc) if (associated(fluxes%seaice_melt)) deallocate(fluxes%seaice_melt) if (associated(fluxes%netMassOut)) deallocate(fluxes%netMassOut) if (associated(fluxes%netMassIn)) deallocate(fluxes%netMassIn) @@ -3635,7 +3821,7 @@ end subroutine deallocate_forcing_type subroutine deallocate_mech_forcing(forces) type(mech_forcing), intent(inout) :: forces !< Forcing fields structure - !if (associated(forces%omega_w2x)) deallocate(forces%omega_w2x) + if (associated(forces%omega_w2x)) deallocate(forces%omega_w2x) if (associated(forces%taux)) deallocate(forces%taux) if (associated(forces%tauy)) deallocate(forces%tauy) if (associated(forces%ustar)) deallocate(forces%ustar) @@ -3682,6 +3868,8 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%vprec, turns, fluxes%vprec) call rotate_array(fluxes_in%lrunoff, turns, fluxes%lrunoff) call rotate_array(fluxes_in%frunoff, turns, fluxes%frunoff) + call rotate_array(fluxes_in%lrunoff_glc, turns, fluxes%lrunoff_glc) + call rotate_array(fluxes_in%frunoff_glc, turns, fluxes%frunoff_glc) call rotate_array(fluxes_in%seaice_melt, turns, fluxes%seaice_melt) call rotate_array(fluxes_in%netMassOut, turns, fluxes%netMassOut) call rotate_array(fluxes_in%netMassIn, turns, fluxes%netMassIn) @@ -3696,6 +3884,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%latent_evap_diag, turns, fluxes%latent_evap_diag) call rotate_array(fluxes_in%latent_fprec_diag, turns, fluxes%latent_fprec_diag) call rotate_array(fluxes_in%latent_frunoff_diag, turns, fluxes%latent_frunoff_diag) + call rotate_array(fluxes_in%latent_frunoff_glc_diag, turns, fluxes%latent_frunoff_glc_diag) endif if (do_salt) then @@ -3708,7 +3897,9 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) call rotate_array(fluxes_in%heat_content_fprec, turns, fluxes%heat_content_fprec) call rotate_array(fluxes_in%heat_content_vprec, turns, fluxes%heat_content_vprec) call rotate_array(fluxes_in%heat_content_lrunoff, turns, fluxes%heat_content_lrunoff) + call rotate_array(fluxes_in%heat_content_lrunoff_glc, turns, fluxes%heat_content_lrunoff_glc) call rotate_array(fluxes_in%heat_content_frunoff, turns, fluxes%heat_content_frunoff) + call rotate_array(fluxes_in%heat_content_frunoff_glc, turns, fluxes%heat_content_frunoff_glc) if (associated (fluxes_in%heat_content_evap)) then call rotate_array(fluxes_in%heat_content_evap, turns, fluxes%heat_content_evap) else @@ -3956,6 +4147,8 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%vprec, G, tmp_scale=US%RZ_T_to_kg_m2s) call homogenize_field_t(fluxes%lrunoff, G, tmp_scale=US%RZ_T_to_kg_m2s) call homogenize_field_t(fluxes%frunoff, G, tmp_scale=US%RZ_T_to_kg_m2s) + call homogenize_field_t(fluxes%lrunoff_glc, G, tmp_scale=US%RZ_T_to_kg_m2s) + call homogenize_field_t(fluxes%frunoff_glc, G, tmp_scale=US%RZ_T_to_kg_m2s) call homogenize_field_t(fluxes%seaice_melt, G, tmp_scale=US%RZ_T_to_kg_m2s) ! These two calls might not be needed. call homogenize_field_t(fluxes%netMassOut, G, tmp_scale=GV%H_to_mks) @@ -3974,6 +4167,7 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%latent_evap_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_fprec_diag, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%latent_frunoff_diag, G, tmp_scale=US%QRZ_T_to_W_m2) + call homogenize_field_t(fluxes%latent_frunoff_glc_diag, G, tmp_scale=US%QRZ_T_to_W_m2) endif if (do_salt) call homogenize_field_t(fluxes%salt_flux, G, tmp_scale=US%RZ_T_to_kg_m2s) @@ -3985,6 +4179,8 @@ subroutine homogenize_forcing(fluxes, G, GV, US) call homogenize_field_t(fluxes%heat_content_vprec, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%heat_content_lrunoff, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%heat_content_frunoff, G, tmp_scale=US%QRZ_T_to_W_m2) + call homogenize_field_t(fluxes%heat_content_lrunoff_glc, G, tmp_scale=US%QRZ_T_to_W_m2) + call homogenize_field_t(fluxes%heat_content_frunoff_glc, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%heat_content_massout, G, tmp_scale=US%QRZ_T_to_W_m2) call homogenize_field_t(fluxes%heat_content_massin, G, tmp_scale=US%QRZ_T_to_W_m2) endif diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index fb95b79a91..fdcee8107d 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -966,8 +966,8 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then do j=js,je ; do i=is,ie FW_in(i,j) = RZL2_to_kg * dt*G%areaT(i,j)*(fluxes%evap(i,j) + & - (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + & - (fluxes%fprec(i,j) + fluxes%frunoff(i,j)))) + (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j) + fluxes%lrunoff_glc(i,j)) + & + (fluxes%fprec(i,j) + fluxes%frunoff(i,j) + fluxes%frunoff_glc(i,j)))) enddo ; enddo else call MOM_error(WARNING, & diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index fec7219e7d..816d5d7498 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -30,6 +30,7 @@ module MOM_CVMix_KPP use CVMix_kpp, only : CVMix_kpp_compute_unresolved_shear use CVMix_kpp, only : CVMix_kpp_params_type use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth +use CVMix_kpp, only : CVMix_kpp_compute_StokesXi implicit none ; private @@ -82,6 +83,7 @@ module MOM_CVMix_KPP logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer. character(len=32) :: interpType !< Type of interpolation to compute bulk Richardson number character(len=32) :: interpType2 !< Type of interpolation to compute diff and visc at OBL_depth + logical :: StokesMOST !< If True, use Stokes similarity package logical :: computeEkman !< If True, compute Ekman depth limit for OBLdepth logical :: computeMoninObukhov !< If True, compute Monin-Obukhov limit for OBLdepth logical :: passiveMode !< If True, makes KPP passive meaning it does NOT alter the diffusivity @@ -145,11 +147,15 @@ module MOM_CVMix_KPP integer :: id_EnhW = -1 integer :: id_La_SL = -1 integer :: id_OBLdepth_original = -1 + integer :: id_StokesXI = -1 + integer :: id_Lam2 = -1 !>@} ! Diagnostics arrays real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of ocean boundary layer (OBL) [Z ~> m] real, allocatable, dimension(:,:) :: OBLdepth_original !< Depth (positive) of OBL [Z ~> m] without smoothing + real, allocatable, dimension(:,:) :: StokesParXI !< Stokes similarity parameter + real, allocatable, dimension(:,:) :: Lam2 !< La^(-2) = Ustk0/u* real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent [nondim] real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL [Z ~> m] real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP [nondim] @@ -272,6 +278,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 'Type of interpolation to compute diff and visc at OBL_depth.\n'// & 'Allowed types are: linear, quadratic, cubic or LMD94.', & default='LMD94') + call get_param(paramFile, mdl, 'STOKES_MOST', CS%StokesMOST, & + 'If True, use Stokes Similarity package.', & + default=.False.) call get_param(paramFile, mdl, 'COMPUTE_EKMAN', CS%computeEkman, & 'If True, limit OBL depth to be no deeper than Ekman depth.', & default=.False.) @@ -498,6 +507,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) interp_type=CS%interpType, & interp_type2=CS%interpType2, & lEkman=CS%computeEkman, & + lStokesMOST=CS%StokesMOST, & lMonOb=CS%computeMoninObukhov, & MatchTechnique=CS%MatchTechnique, & lenhanced_diff=CS%enhance_diffusion,& @@ -524,6 +534,12 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', & cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') endif + if( CS%StokesMOST ) then + CS%id_StokesXI = register_diag_field('ocean_model', 'StokesXI', diag%axesT1, Time, & + 'Stokes Similarity Parameter', 'nondim') + CS%id_Lam2 = register_diag_field('ocean_model', 'Lam2', diag%axesT1, Time, & + 'Ustk0_ustar', 'nondim') + endif CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, & 'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', & 'kg/m3', conversion=US%R_to_kg_m3) @@ -584,6 +600,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) allocate( CS%N( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. ) allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ), source=0. ) + allocate( CS%StokesParXI( SZI_(G), SZJ_(G) ), source=0. ) + allocate( CS%Lam2 ( SZI_(G), SZJ_(G) ), source=0. ) allocate( CS%kOBL( SZI_(G), SZJ_(G) ), source=0. ) allocate( CS%La_SL( SZI_(G), SZJ_(G) ), source=0. ) allocate( CS%Vt2( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. ) @@ -804,6 +822,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, & 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 + StokesXi = CS%StokesParXI(i,j), & ! Stokes forcing parameter CVMix_kpp_params_user=CS%KPP_params ) ! safety check, Kviscosity and Kdiffusivity must be >= 0 @@ -962,7 +981,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m] - ! Variables for passing to CVMix routines, often in MKS units real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1] real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3] @@ -997,6 +1015,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: Uk, Vk ! Layer velocities relative to their averages in the surface layer [L T-1 ~> m s-1] real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth [Z ~> m] real :: hTot ! Running sum of thickness used in the surface layer average [Z ~> m] + real :: I_hTot ! The inverse of hTot [Z-1 ~> m-1] real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: delH ! Thickness of a layer [Z ~> m] real :: surfTemp ! Average of temperature over the surface layer [C ~> degC] @@ -1018,6 +1037,17 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl integer :: i, j, k, km1, kk, ksfc, ktmp ! Loop indices + real, dimension(GV%ke) :: uE_H, vE_H ! Eulerian velocities h-points, centers [L T-1 ~> m s-1] + real, dimension(GV%ke) :: uS_H, vS_H ! Stokes drift components h-points, centers [L T-1 ~> m s-1] + real, dimension(GV%ke) :: uSbar_H, vSbar_H ! Cell Average Stokes drift h-points [L T-1 ~> m s-1] + real, dimension(GV%ke+1) :: uS_Hi, vS_Hi ! Stokes Drift components at interfaces [L T-1 ~> m s-1] + real :: uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD ! Stokes at/to to Surface Layer Extent + ! [L T-1 ~> m s-1] + real :: StokesXI ! Stokes similarity parameter [nondim] + real, dimension( GV%ke ) :: StokesXI_1d , StokesVt_1d ! Parameters of TKE production ratio [nondim] + real :: Llimit ! Stable boundary Layer Limit = vonk Lstar [Z ~> m] + integer :: kbl ! index of cell containing boundary layer depth + if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, & "KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.") @@ -1046,27 +1076,36 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & !$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, vt2_1d, & !$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, & - !$OMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, & + !$OMP surfHvS, hTot, I_hTot, delH, surftemp, surfsalt, surfu, surfv, & !$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, N_col, & !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, & !$OMP deltarho, deltaBuoy, N2_1d, ws_1d, LangEnhVT2,KPP_OBL_depth, z_cell, & - !$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset) & + !$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset, uE_H, vE_H, & + !$OMP uS_H, vS_H, uSbar_H, vSbar_H , uS_Hi, vS_Hi, & + !$OMP uS_SLD, vS_SLD, uS_SLC, vS_SLC, uSbar_SLD, vSbar_SLD, & + !$OMP StokesXI, StokesXI_1d, StokesVt_1d, kbl) & !$OMP shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, & !$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult) do j = G%jsc, G%jec do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then do k=1,GV%ke - U_H(k) = 0.5 * (u(i,j,k)+u(i-1,j,k)) - V_H(k) = 0.5 * (v(i,j,k)+v(i,j-1,k)) + U_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)) + V_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)) enddo + if (CS%StokesMOST) then + do k=1,GV%ke + uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) + vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k)) + enddo + endif ! things independent of position within the column Coriolis = 0.25*US%s_to_T*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) + & (G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) surfFricVel = US%Z_to_m*US%s_to_T * uStar(i,j) - ! Bullk Richardson number computed for each cell in a column, + ! Bulk Richardson number computed for each cell in a column, ! assuming OBLdepth = grid cell depth. After Rib(k) is ! known for the column, then CVMix interpolates to find ! the actual OBLdepth. This approach avoids need to iterate @@ -1075,8 +1114,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl iFaceHeight(1) = 0.0 ! BBL is all relative to the surface pRef = 0. ; if (associated(tv%p_surf)) pRef = tv%p_surf(i,j) hcorr = 0. - do k=1,GV%ke + if (CS%StokesMOST) call Compute_StokesDrift( i, j, h(i,j,1) , iFaceHeight(1), & + uS_Hi(1), vS_Hi(1), uS_H(1), vS_H(1), uSbar_H(1), vSbar_H(1), Waves) + + do k=1,GV%ke ! cell center and cell bottom in meters (negative values in the ocean) dh = dz(i,j,k) ! Nominal thickness to use for increment dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) @@ -1095,53 +1137,99 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif enddo - ! average temperature, salinity, u and v over surface layer - ! use C-grid average to get u and v on T-points. - surfHtemp = 0.0 - surfHsalt = 0.0 - surfHu = 0.0 - surfHv = 0.0 - surfHuS = 0.0 - surfHvS = 0.0 - hTot = 0.0 - do ktmp = 1,ksfc - - ! SLdepth_0d can be between cell interfaces - delH = min( max(0.0, SLdepth_0d - hTot), dz(i,j,ktmp) ) - - ! surface layer thickness - hTot = hTot + delH - - ! surface averaged fields - surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH - surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH - surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH - surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + if (CS%StokesMOST) then + surfBuoyFlux = buoy_scale * & + (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,k)+buoyFlux(i,j,k+1)) ) + surfBuoyFlux2(k) = surfBuoyFlux + call Compute_StokesDrift(i,j, iFaceHeight(k),iFaceHeight(k+1), & + uS_Hi(k+1), vS_Hi(k+1), uS_H(k), vS_H(k), uSbar_H(k), vSbar_H(k), Waves) + call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, & + uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves) + call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d,surfBuoyFlux, & + surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, vS_Hi, uSbar_H, vSbar_H, uS_SLD,& + vS_SLD, uSbar_SLD, vSbar_SLD, StokesXI, CVMix_kpp_params_user=CS%KPP_params ) + StokesXI_1d(k) = StokesXI + StokesVt_1d(k) = StokesXI + + ! average temperature, salinity, u and v over surface layer starting at ksfc + delH = SLdepth_0d + iFaceHeight(ksfc-1) + surfHtemp = Temp(i,j,ksfc) * delH + surfHsalt = Salt(i,j,ksfc) * delH + surfHu = (uE_H(ksfc) + uSbar_SLD) * delH + surfHv = (vE_H(ksfc) + vSbar_SLD) * delH + hTot = delH + do ktmp = 1,ksfc-1 ! if ksfc >=2 + delH = h(i,j,ktmp)*GV%H_to_Z + hTot = hTot + delH + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + (uE_H(ktmp) + uSbar_H(ktmp)) * delH + surfHv = surfHv + (vE_H(ktmp) + vSbar_H(ktmp)) * delH + enddo + I_hTot = 1./hTot + surfTemp = surfHtemp * I_hTot + surfSalt = surfHsalt * I_hTot + surfU = surfHu * I_hTot + surfV = surfHv * I_hTot + Uk = uE_H(k) + uS_H(k) - surfU + Vk = vE_H(k) + vS_H(k) - surfV + + else !not StokesMOST + StokesXI_1d(k) = 0.0 + + ! average temperature, salinity, u and v over surface layer + ! use C-grid average to get u and v on T-points. + surfHtemp = 0.0 + surfHsalt = 0.0 + surfHu = 0.0 + surfHv = 0.0 + surfHuS = 0.0 + surfHvS = 0.0 + hTot = 0.0 + do ktmp = 1,ksfc + + ! SLdepth_0d can be between cell interfaces + delH = min( max(0.0, SLdepth_0d - hTot), dz(i,j,ktmp) ) + + ! surface layer thickness + hTot = hTot + delH + + ! surface averaged fields + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH + surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + if (CS%Stokes_Mixing) then + surfHus = surfHus + 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) * delH + surfHvs = surfHvs + 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) * delH + endif + + enddo + surfTemp = surfHtemp / hTot + surfSalt = surfHsalt / hTot + surfU = surfHu / hTot + surfV = surfHv / hTot + surfUs = surfHus / hTot + surfVs = surfHvs / hTot + + ! vertical shear between present layer and surface layer averaged surfU and surfV. + ! C-grid average to get Uk and Vk on T-points. + Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + if (CS%Stokes_Mixing) then - surfHus = surfHus + 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) * delH - surfHvs = surfHvs + 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) * delH + ! If momentum is mixed down the Stokes drift gradient, then + ! the Stokes drift must be included in the bulk Richardson number + ! calculation. + Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs ) + Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs ) endif - enddo - surfTemp = surfHtemp / hTot - surfSalt = surfHsalt / hTot - surfU = surfHu / hTot - surfV = surfHv / hTot - surfUs = surfHus / hTot - surfVs = surfHvs / hTot - - ! vertical shear between present layer and surface layer averaged surfU and surfV. - ! C-grid average to get Uk and Vk on T-points. - Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU - Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV - - if (CS%Stokes_Mixing) then - ! If momentum is mixed down the Stokes drift gradient, then - ! the Stokes drift must be included in the bulk Richardson number - ! calculation. - Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs ) - Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs ) - endif + ! this difference accounts for penetrating SW + surfBuoyFlux = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1)) + surfBuoyFlux2(k) = surfBuoyFlux + + endif ! StokesMOST deltaU2(k) = US%L_T_to_m_s**2 * (Uk**2 + Vk**2) @@ -1165,9 +1253,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! iterate pRef for next pass through k-loop. pRef = pRef + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k) - ! this difference accounts for penetrating SW - surfBuoyFlux2(k) = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1)) - enddo ! k-loop finishes if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then @@ -1215,11 +1300,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass ! sigma=CS%surf_layer_ext for this calculation. - call CVMix_kpp_compute_turbulent_scales( & + call CVMix_kpp_compute_turbulent_scales( & ! 1d_OBL CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext OBL_depth, & ! (in) OBL depth [m] surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3] surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + xi=StokesVt_1d, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance of Vt w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params ) @@ -1255,10 +1341,17 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl N_iface=N_col, & ! 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] + bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3] uStar=surfFricVel, & ! surface friction velocity [m s-1] CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters +! ! A hack to avoid KPP reaching the bottom. It was needed during development +! ! because KPP was unable to handle vanishingly small layers near the bottom. +! if (CS%deepOBLoffset>0.) then +! zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset, -0.1*iFaceHeight(GV%ke+1)) +! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) +! endif + zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(GV%ke+1)) call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number @@ -1267,11 +1360,39 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent zt_cntr=z_cell, & ! (in) Height of cell centers [m] surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] - surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] + surf_buoy=surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3] Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1] + Xi = StokesXI_1d, & ! (in) Stokes similarity parameter Lmob limit (1-Xi) + zBottom = zBottomMinusOffset, & ! (in) Numerical limit on OBLdepth CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters CS%OBLdepth(i,j) = US%m_to_Z * KPP_OBL_depth + if (CS%StokesMOST) then + kbl = nint(CS%kOBL(i,j)) + SLdepth_0d = CS%surf_layer_ext*CS%OBLdepth(i,j) + surfBuoyFlux = surfBuoyFlux2(kbl) + ! find ksfc for cell where "surface layer" sits + ksfc = kbl + do ktmp = 1, kbl + if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then + ksfc = ktmp + exit + endif + enddo + + call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, & + uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves) + call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d, & + surfBuoyFlux, surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, & + vS_Hi, uSbar_H, vSbar_H, uS_SLD, vS_SLD, uSbar_SLD, vSbar_SLD, & + StokesXI, CVMix_kpp_params_user=CS%KPP_params ) + CS%StokesParXI(i,j) = StokesXI + CS%Lam2(i,j) = sqrt(US_Hi(1)**2+VS_Hi(1)**2) / MAX(surfFricVel,0.0002) + + else !.not Stokes_MOST + CS%StokesParXI(i,j) = 10.0 + CS%Lam2(i,j) = sqrt(US_Hi(1)**2+VS_Hi(1)**2) / MAX(surfFricVel,0.0002) + ! A hack to avoid KPP reaching the bottom. It was needed during development ! because KPP was unable to handle vanishingly small layers near the bottom. if (CS%deepOBLoffset>0.) then @@ -1285,6 +1406,8 @@ 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) ) + endif !Stokes_MOST + ! compute unresolved squared velocity for diagnostics if (CS%id_Vt2 > 0) then Vt2_1d(:) = CVmix_kpp_compute_unresolved_shear( & @@ -1293,7 +1416,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl N_iface=N_col, & ! 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] + bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3] uStar=surfFricVel, & ! surface friction velocity [m s-1] CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters CS%Vt2(i,j,:) = US%m_to_Z*US%T_to_s * Vt2_1d(:) @@ -1307,6 +1430,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl US%Z_to_m*CS%OBLdepth(i,j), & ! (in) OBL depth [m] surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + xi=StokesXI, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params) ! KPP parameters CS%Ws(i,j,:) = US%m_to_Z*US%T_to_s*Ws_1d(:) @@ -1342,6 +1466,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl 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) + if (CS%StokesMOST) then + if (CS%id_StokesXI > 0) call post_data(CS%id_StokesXI, CS%StokesParXI, CS%diag) + if (CS%id_Lam2 > 0) call post_data(CS%id_Lam2 , CS%Lam2 , CS%diag) + endif + ! BLD smoothing: if (CS%n_smooth > 0) call KPP_smooth_BLD(CS, G, GV, US, dz) @@ -1594,6 +1723,49 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, end subroutine KPP_NonLocalTransport_saln +!> Compute Stokes Drift components at zbot < ztop <= 0 and at k=0.5*(ztop+zbot) and +!! average components from ztop to zbot <= 0 +subroutine Compute_StokesDrift(i ,j, ztop, zbot, uS_i, vS_i, uS_k, vS_k, uSbar, vSbar, waves) + + type(wave_parameters_CS), pointer :: waves !< Wave CS for Langmuir turbulence + real, intent(in) :: ztop !< cell top + real, intent(in) :: zbot !< cell bottom + real, intent(inout) :: uS_i !< Stokes u velocity at zbot interface + real, intent(inout) :: vS_i !< Stokes v velocity at zbot interface + real, intent(inout) :: uS_k !< Stokes u velocity at zk center + real, intent(inout) :: vS_k !< Stokes v at zk =0.5(ztop+zbot) + real, intent(inout) :: uSbar !< mean Stokes u (ztop to zbot) + real, intent(inout) :: vSbar !< mean Stokes v (ztop to zbot) + integer, intent(in) :: i !< Meridional index of H-point + integer, intent(in) :: j !< Zonal index of H-point + + ! local variables + integer :: b !< wavenumber band index + real :: fexp !< an exponential function + real :: WaveNum !< Wavenumber + + uS_i = 0.0 + vS_i = 0.0 + uS_k = 0.0 + vS_k = 0.0 + uSbar = 0.0 + vSbar = 0.0 + do b = 1, waves%NumBands + WaveNum = waves%WaveNum_Cen(b) + fexp = exp(2. * WaveNum * zbot) + uS_i = uS_i + waves%Ustk_Hb(i,j,b) * fexp + vS_i = vS_i + waves%Vstk_Hb(i,j,b) * fexp + fexp = exp( WaveNum * (ztop + zbot) ) + uS_k = uS_k+ waves%Ustk_Hb(i,j,b) * fexp + vS_k = vS_k+ waves%Vstk_Hb(i,j,b) * fexp + fexp = exp(2. * WaveNum * ztop) - exp(2. * WaveNum * zbot) + uSbar = uSbar + 0.5 * waves%Ustk_Hb(i,j,b) * fexp / WaveNum + vSbar = vSbar + 0.5 * waves%Vstk_Hb(i,j,b) * fexp / WaveNum + enddo + uSbar = uSbar / (ztop-zbot) + vSbar = vSbar / (ztop-zbot) + +end subroutine Compute_StokesDrift !> Clear pointers, deallocate memory subroutine KPP_end(CS) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 561ace60a7..e3560dc03e 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -528,13 +528,15 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C RmixConst = -0.5*CS%rivermix_depth * GV%g_Earth do i=is,ie TKE_river(i) = max(0.0, RmixConst * dSpV0_dS(i) * & - (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * S(i,1)) + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) + & + fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j)) * S(i,1)) enddo else RmixConst = 0.5*CS%rivermix_depth * GV%g_Earth * Irho0**2 do i=is,ie TKE_river(i) = max(0.0, RmixConst*dR0_dS(i)* & - (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * S(i,1)) + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) + & + fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j)) * S(i,1)) enddo endif else diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index aa31024b24..c54240aae2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1431,7 +1431,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t RivermixConst = -0.5*(CS%rivermix_depth*dt) * GV%Rho0 * ( US%L_to_Z**2*GV%g_Earth ) endif cTKE(i,j,k) = cTKE(i,j,k) + max(0.0, RivermixConst*dSV_dS(i,j,1) * & - (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1)) + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) + & + fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j)) * tv%S(i,j,1)) endif ! Update state diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c26ee4ac75..3f968b2101 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -31,6 +31,8 @@ module MOM_vert_friction use MOM_set_visc, only : set_v_at_u, set_u_at_v use MOM_lateral_mixing_coeffs, only : VarMix_CS +use CVMix_kpp, only : cvmix_kpp_composite_Gshape + implicit none ; private #include @@ -170,9 +172,11 @@ module MOM_vert_friction integer :: id_au_vv = -1, id_av_vv = -1, id_au_gl90_vv = -1, id_av_gl90_vv = -1 integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 - integer :: id_FPw2x = -1 !W id_FPhbl_u = -1, id_FPhbl_v = -1 - integer :: id_tauFP_u = -1, id_tauFP_v = -1 !W, id_FPtau2x_u = -1, id_FPtau2x_v = -1 - integer :: id_FPtau2s_u = -1, id_FPtau2s_v = -1, id_FPtau2w_u = -1, id_FPtau2w_v = -1 + integer :: id_Omega_w2x = -1, id_FPtau2s = -1 , id_FPtau2w = -1 + integer :: id_uE_h = -1, id_vE_h = -1 + integer :: id_uStk = -1, id_vStk = -1 + integer :: id_uStk0 = -1, id_vStk0 = -1 + integer :: id_uInc_h= -1, id_vInc_h= -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 integer :: id_Kv_gl90_u = -1, id_Kv_gl90_v = -1 @@ -191,392 +195,211 @@ module MOM_vert_friction contains -!> Add nonlocal stress increments to u^n (uold) and v^n (vold) using ui and vi. -subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) +!> Add nonlocal stress increments to ui^n and vi^n. +subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G, GV, US, CS, OBC, Waves) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1] + intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: uold !< Old Zonal velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(inout) :: vold !< Old Meridional velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h !< boundary layer depth [H ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces - real, intent(in) :: dt !< Time increment [T ~> s] - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces + real, intent(in) :: dt !< Time increment [T ~> s] + real, intent(in) :: Cemp_NL !< empirical coefficient of non-local momentum mixing [nondim] + logical, intent(in) :: lpost !< Compute and make available FPMix diagnostics + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure + type(wave_parameters_CS), & + optional, pointer :: Waves !< Container for wave/Stokes information ! local variables - real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth at u-pts [H ~> m] - real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth at v-pts [H ~> m] - integer, dimension(SZIB_(G),SZJ_(G)) :: kbl_u !< index of the BLD at u-pts [nondim] - integer, dimension(SZI_(G),SZJB_(G)) :: kbl_v !< index of the BLD at v-pts [nondim] - real, dimension(SZIB_(G),SZJ_(G)) :: ustar2_u !< ustar squared at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G)) :: ustar2_v !< ustar squared at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< zonal wind stress at u-pts [R L Z T-2 ~> Pa] - real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< meridional wind stress at v-pts [R L Z T-2 ~> Pa] - !real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u !< angle between wind and x-axis at u-pts [rad] - !real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v !< angle between wind and y-axis at v-pts [rad] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u !< kinematic zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tau_v !< kinematic mer. mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauxDG_u !< downgradient zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauyDG_u !< downgradient meri mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauxDG_v !< downgradient zonal mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauyDG_v !< downgradient meri mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2s_u !< angle between mtm flux and vert shear at u-pts [rad] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2s_v !< angle between mtm flux and vert shear at v-pts [rad] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u !< angle between mtm flux and wind at u-pts [rad] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v !< angle between mtm flux and wind at v-pts [rad] - - real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp !< constants and dummy variables [nondim] - real :: omega_tmp !< A dummy angle [radians] - real :: du, dv !< Velocity increments [L T-1 ~> m s-1] - real :: depth !< Cumulative layer thicknesses [H ~> m or kg m=2] - real :: sigma !< Fractional depth in the mixed layer [nondim] - real :: Wind_x, Wind_y !< intermediate wind stress componenents [L2 T-2 ~> m2 s-2] - real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables [L2 T-2 ~> m2 s-2] - real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables [L2 T-2 ~> m2 s-2] - real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles [radians] - integer :: kblmin, kbld, kp1, k, nz !< vertical indices - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices + real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth (u-pts) [H ~> m] + real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth (v-pts) [H ~> m] + real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< kinematic zonal wind stress (u-pts) [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< kinematic merid wind stress (v-pts) [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: uS0 !< surface zonal Stokes drift [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G)) :: vS0 !< surface zonal Stokes drift [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uE_u !< zonal Eulerian u-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uE_h !< zonal Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vE_v !< merid Eulerian v-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vE_h !< merid Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uInc_u !< zonal Eulerian u-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uInc_h !< zonal Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vInc_v !< merid Eulerian v-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vInc_h !< merid Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uStk !< zonal Stokes Drift (h-pts) [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vStk !< merid Stokes Drift (h-pts) [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2s !< angle stress to shear (h-pts) [rad] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2w !< angle stress to wind (h-pts) [rad] + real :: pi, tmp_u, tmp_v, omega_tmp, Irho0, fexp !< constants and dummy variables + real :: sigma,Gat1,Gsig,dGdsig !< Shape parameters + real :: du, dv, depth, Wind_x, Wind_y !< intermediate variables + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, ustar2min, tauh !< intermediate variables + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables + real :: omega_w2s, omega_s2x, omega_tau2x, omega_s2w , omega_e2x, omega_l2x !< intermediate angles + integer :: b, kbld, kp1, k, nz !< band and vertical indices + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq !< horizontal indices is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke pi = 4. * atan2(1.,1.) - Cemp_CG = 3.6 - kblmin = 1 - taux_u(:,:) = 0. - tauy_v(:,:) = 0. + Irho0 = 1.0 / GV%Rho0 - do j = js,je - do I = Isq,Ieq - taux_u(I,j) = forces%taux(I,j) / GV%H_to_RZ !W rho0=1035. - enddo - enddo - - do J = Jsq,Jeq - do i = is,ie - tauy_v(i,J) = forces%tauy(i,J) / GV%H_to_RZ - enddo - enddo + call pass_var(hbl_h , G%Domain, halo=1) - call pass_var( hbl_h ,G%Domain, halo=1 ) - call pass_vector(taux_u , tauy_v, G%Domain, To_All ) - ustar2_u(:,:) = 0. - ustar2_v(:,:) = 0. - hbl_u(:,:) = 0. - hbl_v(:,:) = 0. - kbl_u(:,:) = 0 - kbl_v(:,:) = 0 - !omega_w2x_u(:,:) = 0.0 - !omega_w2x_v(:,:) = 0.0 - tauxDG_u(:,:,:) = 0.0 - tauyDG_v(:,:,:) = 0.0 + ! u-points do j = js,je do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - tmp = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) - hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j) * hbl_h(i+1,j)) /tmp - tmp = MAX(1.0, (G%mask2dCv(i,j) + G%mask2dCv(i,j-1) + G%mask2dCv(i+1,j) + G%mask2dCv(i+1,j-1) ) ) - tauy = ( G%mask2dCv(i ,j )*tauy_v(i ,j ) + G%mask2dCv(i ,j-1)*tauy_v(i ,j-1) & - + G%mask2dCv(i+1,j )*tauy_v(i+1,j ) + G%mask2dCv(i+1,j-1)*tauy_v(i+1,j-1) ) / tmp - ustar2_u(I,j) = sqrt( taux_u(I,j)*taux_u(I,j) + tauy*tauy ) - !omega_w2x_u(I,j) = atan2( tauy , taux_u(I,j) ) - tauxDG_u(I,j,1) = taux_u(I,j) - depth = 0.0 - do k = 1, nz - depth = depth + CS%h_u(I,j,k) - if( (depth >= hbl_u(I,j)) .and. (kbl_u(I,j) == 0 ) .and. (k > (kblmin-1)) ) then - kbl_u(I,j) = k - hbl_u(I,j) = depth - endif - enddo + taux_u(I,j) = forces%taux(I,j) * Irho0 + if ( (G%mask2dCu(I,j) > 0.5) ) then + ! h to u-pts + tmp_u = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) + hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j) * hbl_h(i+1,j)) / tmp_u endif + depth = 0. + Gat1 = 0. + do k=1, nz + ! cell center + depth = depth + 0.5*CS%h_u(I,j,k) + uE_u(I,j,k) = ui(I,j,k) - waves%Us_x(I,j,k) + if ( depth < hbl_u(I,j) ) then + sigma = depth / hbl_u(i,j) + ! cell bottom + depth = depth + 0.5*CS%h_u(I,j,k) + call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig) + ! nonlocal boundary-layer increment + uInc_u(I,j,k) = dt * Cemp_NL * taux_u(I,j) * dGdsig / hbl_u(I,j) + ui(I,j,k) = ui(I,j,k) + uInc_u(I,j,k) + else + uInc_u(I,j,k) = 0.0 + endif + enddo enddo enddo + + ! v-points do J = Jsq,Jeq do i = is,ie - if( (G%mask2dCv(i,J) > 0.5) ) then - tmp = max( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1))) - hbl_v(i,J) = (G%mask2dT(i,j) * hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) /tmp - tmp = max(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i,j+1) + G%mask2dCu(i-1,j) + G%mask2dCu(i-1,j+1))) - taux = ( G%mask2dCu(i ,j) * taux_u(i ,j) + G%mask2dCu(i ,j+1) * taux_u(i ,j+1) & - + G%mask2dCu(i-1,j) * taux_u(i-1,j) + G%mask2dCu(i-1,j+1) * taux_u(i-1,j+1)) / tmp - ustar2_v(i,J) = sqrt(tauy_v(i,J)*tauy_v(i,J) + taux*taux) - !omega_w2x_v(i,J) = atan2( tauy_v(i,J), taux ) - tauyDG_v(i,J,1) = tauy_v(i,J) - depth = 0.0 - do k = 1, nz - depth = depth + CS%h_v(i,J,k) - if( (depth >= hbl_v(i,J)) .and. (kbl_v(i,J) == 0) .and. (k > (kblmin-1))) then - kbl_v(i,J) = k - hbl_v(i,J) = depth - endif - enddo + tauy_v(i,J) = forces%tauy(i,J) * Irho0 + if ( (G%mask2dCv(i,J) > 0.5) ) then + ! h to v-pts + tmp_v = max( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1))) + hbl_v(i,J) = (G%mask2dT(i,j) * hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) / tmp_v endif - enddo - enddo - - if (CS%debug) then - !### These checksum calls are missing necessary dimensional scaling factors. - call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=1, scalar_pair=.true.) - call uvchksum("ustar2", ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) - call uvchksum(" hbl", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) - endif - - ! Compute downgradient stresses - do k = 1, nz - kp1 = min( k+1 , nz) - do j = js ,je - do I = Isq , Ieq - tauxDG_u(I,j,k+1) = CS%a_u(I,j,kp1) * (ui(I,j,k) - ui(I,j,kp1)) - enddo - enddo - do J = Jsq , Jeq - do i = is , ie - tauyDG_v(i,J,k+1) = CS%a_v(i,J,kp1) * (vi(i,J,k) - vi(i,J,kp1)) - enddo - enddo - enddo - - call pass_vector(tauxDG_u, tauyDG_v , G%Domain, To_All) - call pass_vector(ui,vi, G%Domain, To_All) - tauxDG_v(:,:,:) = 0. - tauyDG_u(:,:,:) = 0. - - ! Thickness weighted interpolations - do k = 1, nz - ! v to u points - do j = js , je - do I = Isq, Ieq - tauyDG_u(I,j,k) = set_v_at_u(tauyDG_v, h, G, GV, I, j, k, G%mask2dCv, OBC) - enddo - enddo - ! u to v points - do J = Jsq, Jeq - do i = is, ie - tauxDG_v(i,J,k) = set_u_at_v(tauxDG_u, h, G, GV, i, J, k, G%mask2dCu, OBC) + depth = 0. + Gat1 = 0. + do k=1, nz + ! cell center + depth = depth + 0.5* CS%h_v(i,J,k) + vE_v(i,J,k) = vi(i,J,k) - waves%Us_y(i,J,k) + if ( depth < hbl_v(i,J) ) then + sigma = depth / hbl_v(i,J) + ! cell bottom + depth = depth + 0.5* CS%h_v(i,J,k) + call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig) + ! nonlocal boundary-layer increment + vInc_v(i,J,k) = dt * Cemp_NL * tauy_v(i,J) * dGdsig / hbl_v(i,J) + vi(i,J,k) = vi(i,J,k) + vInc_v(i,J,k) + else + vInc_v(i,J,k) = 0.0 + endif enddo enddo enddo - if (CS%debug) then - call uvchksum(" tauyDG_u tauxDG_v",tauyDG_u,tauxDG_v, G%HI, haloshift=0, scalar_pair=.true.) - endif - ! compute angles, tau2x_[u,v], tau2w_[u,v], tau2s_[u,v], s2w_[u,v] and stress mag tau_[u,v] - omega_tau2w_u(:,:,:) = 0.0 - omega_tau2w_v(:,:,:) = 0.0 - omega_tau2s_u(:,:,:) = 0.0 - omega_tau2s_v(:,:,:) = 0.0 - tau_u(:,:,:) = 0.0 - tau_v(:,:,:) = 0.0 - - ! stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] - do j = js,je - do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - ! SURFACE - tauyDG_u(I,j,1) = ustar2_u(I,j) !* cos(omega_w2x_u(I,j)) - tau_u(I,j,1) = ustar2_u(I,j) - Omega_tau2w_u(I,j,1) = 0.0 - Omega_tau2s_u(I,j,1) = 0.0 + ! Compute and store diagnostics, only during the corrector step. + if (lpost) then + call pass_vector(uE_u , vE_v , G%Domain, To_All) + call pass_vector(uInc_u, vInc_v , G%Domain, To_All) + uStk = 0.0 + vStk = 0.0 + uS0 = 0.0 + vS0 = 0.0 + + do j = js,je + do i = is,ie + if (G%mask2dT(i,j) > 0.5) then + ! u to h-pts + tmp_u = max( 1.0 ,(G%mask2dCu(i,j) + G%mask2dCu(i-1,j))) + ! v to h-pts + tmp_v = max( 1.0 ,(G%mask2dCv(i,j) + G%mask2dCv(i,j-1))) + do k = 1,nz + uE_h(i,j,k) = (G%mask2dCu(i,j) * uE_u(i,j,k) + G%mask2dCu(i-1,j) * uE_u(i-1,j,k)) / tmp_u + uInc_h(i,j,k) = (G%mask2dCu(i,j) * uInc_u(i,j,k) + G%mask2dCu(i-1,j) * uInc_u(i-1,j,k)) / tmp_u + vE_h(i,j,k) = (G%mask2dCv(i,j) * vE_v(i,j,k) + G%mask2dCv(i,j-1) * vE_v(i,j-1,k)) / tmp_v + vInc_h(i,j,k) = (G%mask2dCv(i,j) * vInc_v(i,j,k) + G%mask2dCv(i,j-1) * vInc_v(i,j-1,k)) / tmp_v + enddo + ! Wind, Stress and Shear align at surface + Omega_tau2w(i,j,1) = 0.0 + Omega_tau2s(i,j,1) = 0.0 + do k = 1,nz + kp1 = min( nz , k+1) + du = uE_h(i,j,k) - uE_h(i,j,kp1) + dv = vE_h(i,j,k) - vE_h(i,j,kp1) + omega_s2x = atan2( dv , du ) + + du = du + uInc_h(i,j,k) - uInc_h(i,j,kp1) + dv = dv + vInc_h(i,j,k) - vInc_h(i,j,kp1) + omega_tau2x = atan2( dv , du ) + + omega_tmp = omega_tau2x - forces%omega_w2x(i,j) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2w(i,j,kp1) = omega_tmp + + omega_tmp = omega_tau2x - omega_s2x + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2s(i,j,kp1) = omega_tmp - do k=1,nz - kp1 = MIN(k+1 , nz) - tau_u(I,j,k+1) = sqrt( tauxDG_u(I,j,k+1)*tauxDG_u(I,j,k+1) + tauyDG_u(I,j,k+1)*tauyDG_u(I,j,k+1)) - Omega_tau2x = atan2( tauyDG_u(I,j,k+1) , tauxDG_u(I,j,k+1) ) - omega_tmp = Omega_tau2x !- omega_w2x_u(I,j) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2w_u(I,j,k+1) = omega_tmp - Omega_tau2s_u(I,j,k+1) = 0.0 - enddo - endif - enddo - enddo - do J = Jsq, Jeq - do i = is, ie - if( (G%mask2dCv(i,J) > 0.5) ) then - ! SURFACE - tauxDG_v(i,J,1) = ustar2_v(i,J) !* sin(omega_w2x_v(i,J)) - tau_v(i,J,1) = ustar2_v(i,J) - Omega_tau2w_v(i,J,1) = 0.0 - Omega_tau2s_v(i,J,1) = 0.0 - - do k=1,nz-1 - kp1 = MIN(k+1 , nz) - tau_v(i,J,k+1) = sqrt ( tauxDG_v(i,J,k+1)*tauxDG_v(i,J,k+1) + tauyDG_v(i,J,k+1)*tauyDG_v(i,J,k+1) ) - omega_tau2x = atan2( tauyDG_v(i,J,k+1) , tauxDG_v(i,J,k+1) ) - omega_tmp = omega_tau2x !- omega_w2x_v(i,J) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2w_v(i,J,k+1) = omega_tmp - Omega_tau2s_v(i,J,k+1) = 0.0 - enddo - endif - enddo - enddo + enddo + endif - ! Parameterized stress orientation from the wind at interfaces (tau2x) - ! and centers (tau2x) OVERWRITE to kbl-interface above hbl - do j = js,je - do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - kbld = min( (kbl_u(I,j)) , (nz-2) ) - if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 - - !### This expression is dimensionally inconsistent. - tauh = tau_u(I,j,kbld+1) + GV%H_subroundoff - ! surface boundary conditions - depth = 0. - tauNLup = 0.0 - do k=1, kbld - depth = depth + CS%h_u(I,j,k) - sigma = min( 1.0 , depth / hbl_u(i,j) ) - - ! linear stress mag - tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) - !### The following expressions are dimensionally inconsistent. - cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - - ! rotate to wind coordinates - Wind_x = ustar2_u(I,j) !* cos(omega_w2x_u(I,j)) - Wind_y = ustar2_u(I,j) !* sin(omega_w2x_u(I,j)) - tauNL_DG = (Wind_x * cos_tmp + Wind_y * sin_tmp) - tauNL_CG = (Wind_y * cos_tmp - Wind_x * sin_tmp) - omega_w2s = atan2(tauNL_CG, tauNL_DG) - omega_s2w = 0.0-omega_w2s - tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG - tau_MAG = max(tau_MAG, tauNL_CG) - tauNL_DG = sqrt(tau_MAG*tau_MAG - tauNL_CG*tauNL_CG) - tau_u(I,j,k+1) - - ! back to x,y coordinates - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp) - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp) - tauNLdn = tauNL_X - - ! nonlocal increment and update to uold - !### The following expression is dimensionally inconsistent and missing parentheses. - du = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) - ui(I,j,k) = uold(I,j,k) + du - uold(I,j,k) = du - tauNLup = tauNLdn - - ! diagnostics - Omega_tau2s_u(I,j,k+1) = atan2(tauNL_CG , (tau_u(I,j,k+1)+tauNL_DG)) - tau_u(I,j,k+1) = sqrt((tauxDG_u(I,j,k+1) + tauNL_X)**2 + (tauyDG_u(I,j,k+1) + tauNL_Y)**2) - omega_tau2x = atan2((tauyDG_u(I,j,k+1) + tauNL_Y), (tauxDG_u(I,j,k+1) + tauNL_X)) - omega_tau2w = omega_tau2x !- omega_w2x_u(I,j) - if (omega_tau2w >= pi ) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w <= (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - Omega_tau2w_u(I,j,k+1) = omega_tau2w + ! Stokes drift + do b=1,waves%NumBands + uS0(i,j) = uS0(i,j) + waves%UStk_Hb(i,j,b) ! or forces%UStkb(i,j,b) + vS0(i,j) = vS0(i,j) + waves%VStk_Hb(i,j,b) ! or forces%VStkb(i,j,b) enddo - do k= kbld+1, nz - ui(I,j,k) = uold(I,j,k) - uold(I,j,k) = 0.0 + depth = 0.0 + do k = 1,nz + do b = 1, waves%NumBands + ! cell center + fexp = exp(-2. * waves%WaveNum_Cen(b) * (depth+0.5*h(i,j,k)) ) + uStk(i,j,k) = uStk(i,j,k) + waves%UStk_Hb(i,j,b) * fexp + vStk(i,j,k) = vStk(i,j,k) + waves%VStk_Hb(i,j,b) * fexp + enddo + ! cell bottom + depth = depth + h(i,j,k) enddo - endif + enddo enddo - enddo - ! v-point dv increment - do J = Jsq,Jeq - do i = is,ie - if( (G%mask2dCv(i,J) > 0.5) ) then - kbld = min((kbl_v(i,J)), (nz-2)) - if (tau_v(i,J,kbld+2) > tau_v(i,J,kbld+1)) kbld = kbld + 1 - tauh = tau_v(i,J,kbld+1) - - !surface boundary conditions - depth = 0. - tauNLup = 0.0 - do k=1, kbld - depth = depth + CS%h_v(i,J,k) - sigma = min(1.0, depth/ hbl_v(I,J)) - - ! linear stress - tau_MAG = (ustar2_v(i,J) * (1.-sigma)) + (tauh * sigma) - !### The following expressions are dimensionally inconsistent. - cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - - ! rotate into wind coordinate - Wind_x = ustar2_v(i,J) !* cos(omega_w2x_v(i,J)) - Wind_y = ustar2_v(i,J) !* sin(omega_w2x_v(i,J)) - tauNL_DG = (Wind_x * cos_tmp + Wind_y * sin_tmp) - tauNL_CG = (Wind_y * cos_tmp - Wind_x * sin_tmp) - omega_w2s = atan2(tauNL_CG , tauNL_DG) - omega_s2w = 0.0 - omega_w2s - tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG - tau_MAG = max( tau_MAG , tauNL_CG ) - tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt(tau_MAG*tau_MAG - tauNL_CG*tauNL_CG) - - ! back to x,y coordinate - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp) - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp) - tauNLdn = tauNL_Y - !### The following expression is dimensionally inconsistent, [L T-1] vs. [L2 H-1 T-1] on the right, - ! and it is inconsistent with the counterpart expression for du. - dv = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) - vi(i,J,k) = vold(i,J,k) + dv - vold(i,J,k) = dv - tauNLup = tauNLdn - - ! diagnostics - Omega_tau2s_v(i,J,k+1) = atan2(tauNL_CG, tau_v(i,J,k+1) + tauNL_DG) - tau_v(i,J,k+1) = sqrt((tauxDG_v(i,J,k+1) + tauNL_X)**2 + (tauyDG_v(i,J,k+1) + tauNL_Y)**2) - !omega_tau2x = atan2((tauyDG_v(i,J,k+1) + tauNL_Y) , (tauxDG_v(i,J,k+1) + tauNL_X)) - !omega_tau2w = omega_tau2x - omega_w2x_v(i,J) - if (omega_tau2w > pi) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - Omega_tau2w_v(i,J,k+1) = omega_tau2w - enddo + ! post FPmix diagnostics + if (CS%id_uE_h > 0) call post_data(CS%id_uE_h , uE_h , CS%diag) + if (CS%id_vE_h > 0) call post_data(CS%id_vE_h , vE_h , CS%diag) + if (CS%id_uInc_h > 0) call post_data(CS%id_uInc_h , uInc_h , CS%diag) + if (CS%id_vInc_h > 0) call post_data(CS%id_vInc_h , vInc_h , CS%diag) + if (CS%id_FPtau2s > 0) call post_data(CS%id_FPtau2s, Omega_tau2s, CS%diag) + if (CS%id_FPtau2w > 0) call post_data(CS%id_FPtau2w, Omega_tau2w, CS%diag) + if (CS%id_uStk0 > 0) call post_data(CS%id_uStk0 , uS0 , CS%diag) + if (CS%id_vStk0 > 0) call post_data(CS%id_vStk0 , vS0 , CS%diag) + if (CS%id_uStk > 0) call post_data(CS%id_uStk , uStk , CS%diag) + if (CS%id_vStk > 0) call post_data(CS%id_vStk , vStk , CS%diag) + if (CS%id_Omega_w2x > 0) call post_data(CS%id_Omega_w2x, forces%omega_w2x, CS%diag) - do k= kbld+1, nz - vi(i,J,k) = vold(i,J,k) - vold(i,J,k) = 0.0 - enddo - endif - enddo - enddo - - if (CS%debug) then - call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) endif - if (CS%id_tauFP_u > 0) call post_data(CS%id_tauFP_u, tau_u, CS%diag) - if (CS%id_tauFP_v > 0) call post_data(CS%id_tauFP_v, tau_v, CS%diag) - if (CS%id_FPtau2s_u > 0) call post_data(CS%id_FPtau2s_u, omega_tau2s_u, CS%diag) - if (CS%id_FPtau2s_v > 0) call post_data(CS%id_FPtau2s_v, omega_tau2s_v, CS%diag) - if (CS%id_FPtau2w_u > 0) call post_data(CS%id_FPtau2w_u, omega_tau2w_u, CS%diag) - if (CS%id_FPtau2w_v > 0) call post_data(CS%id_FPtau2w_v, omega_tau2w_v, CS%diag) - !if (CS%id_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) - end subroutine vertFPmix -!> Returns the empirical shape-function given sigma [nondim] -real function G_sig(sigma) - real , intent(in) :: sigma !< Normalized boundary layer depth [nondim] - - ! local variables - real :: p1, c2, c3 !< Parameters used to fit and match empirical shape-functions [nondim] - - ! parabola - p1 = 0.287 - ! cubic function - c2 = 1.74392 - c3 = 2.58538 - G_sig = min( p1 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (c2*sigma - c3) ) ) -end function G_sig - !> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb !! (1990), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme !! redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, @@ -701,7 +524,7 @@ end subroutine find_coupling_coef_gl90 !! if DIRECT_STRESS is true, applied to the surface layer. subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & - taux_bot, tauy_bot, Waves) + taux_bot, tauy_bot, fpmix, Waves) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -725,6 +548,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, dimension(SZI_(G),SZJB_(G)), & optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to !! rock [R L Z T-2 ~> Pa] + logical, optional, intent(in) :: fpmix !< fpmix along Eulerian shear type(wave_parameters_CS), & optional, pointer :: Waves !< Container for wave/Stokes information @@ -765,6 +589,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & logical :: do_i(SZIB_(G)) logical :: DoStokesMixing + logical :: lfpmix integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec @@ -802,6 +627,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call MOM_error(FATAL,"Stokes Mixing called without allocated"//& "Waves Control Structure") endif + lfpmix = .false. + if ( present(fpmix) ) lfpmix = fpmix do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo @@ -814,11 +641,17 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo + ! WGL: Brandon Reichl says the following is obsolete. u(I,j,k) already + ! includes Stokes. ! When mixing down Eulerian current + Stokes drift add before calling solver if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) + enddo ; enddo ; endif + if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif @@ -976,6 +809,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) + enddo ; enddo ; endif + enddo ! end u-component j loop ! Now work on the meridional velocity component. @@ -991,6 +828,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(i)) v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,j,k) = v(i,j,k) - Waves%Us_y(i,j,k) + enddo ; enddo ; endif + if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif @@ -1119,6 +960,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(i)) v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,J,k) = v(i,J,k) + Waves%Us_y(i,J,k) + enddo ; enddo ; endif + enddo ! end of v-component J loop ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3]. @@ -2618,7 +2463,7 @@ end subroutine vertvisc_limit_vel !> Initialize the vertical friction module subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & - ntrunc, CS) + ntrunc, CS, fpmix) type(ocean_internal_state), & target, intent(in) :: MIS !< The "MOM Internal State", a set of pointers !! to the fields and accelerations that make @@ -2633,6 +2478,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & type(directories), intent(in) :: dirs !< Relevant directory paths integer, target, intent(inout) :: ntrunc !< Number of velocity truncations type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure + logical, optional, intent(in) :: fpmix !< Nonlocal momentum mixing ! Local variables @@ -2640,6 +2486,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & real :: Kv_back_z ! A background kinematic viscosity [Z2 T-1 ~> m2 s-1] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz + logical :: lfpmix character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname ! This include declares and sets the variable "version". # include "version_variable.h" @@ -2664,6 +2511,9 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & CS%diag => diag ; CS%ntrunc => ntrunc ; ntrunc = 0 + lfpmix = .false. + if (present(fpmix)) lfpmix = fpmix + ! Default, read and log parameters call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & @@ -2966,20 +2816,29 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=US%Z_to_m) - CS%id_FPw2x = register_diag_field('ocean_model', 'FPw2x', diag%axesT1, Time, & - 'Wind direction from x-axis','radians') - CS%id_tauFP_u = register_diag_field('ocean_model', 'tauFP_u', diag%axesCui, Time, & - 'Stress Mag Profile (u-points)', 'm2 s-2') - CS%id_tauFP_v = register_diag_field('ocean_model', 'tauFP_v', diag%axesCvi, Time, & - 'Stress Mag Profile (v-points)', 'm2 s-2') - CS%id_FPtau2s_u = register_diag_field('ocean_model', 'FPtau2s_u', diag%axesCui, Time, & - 'stress from shear direction (u-points)', 'radians ') - CS%id_FPtau2s_v = register_diag_field('ocean_model', 'FPtau2s_v', diag%axesCvi, Time, & - 'stress from shear direction (v-points)', 'radians') - CS%id_FPtau2w_u = register_diag_field('ocean_model', 'FPtau2w_u', diag%axesCui, Time, & - 'stress from wind direction (u-points)', 'radians') - CS%id_FPtau2w_v = register_diag_field('ocean_model', 'FPtau2w_v', diag%axesCvi, Time, & - 'stress from wind direction (v-points)', 'radians') + if (lfpmix) then + CS%id_uE_h = register_diag_field('ocean_model', 'uE_h' , CS%diag%axesTL, & + Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vE_h = register_diag_field('ocean_model', 'vE_h' , CS%diag%axesTL, & + Time, 'y-merid Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_uInc_h = register_diag_field('ocean_model','uInc_h',CS%diag%axesTL, & + Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vInc_h = register_diag_field('ocean_model','vInc_h',CS%diag%axesTL, & + Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_uStk = register_diag_field('ocean_model', 'uStk' , CS%diag%axesTL, & + Time, 'x-FP du increment' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vStk = register_diag_field('ocean_model', 'vStk' , CS%diag%axesTL, & + Time, 'y-FP dv increment' , 'm s-1', conversion=US%L_T_to_m_s) + + CS%id_FPtau2s = register_diag_field('ocean_model','Omega_tau2s',CS%diag%axesTi, & + Time, 'Stress direction from shear','radians') + CS%id_FPtau2w = register_diag_field('ocean_model','Omega_tau2w',CS%diag%axesTi, & + Time, 'Stress direction from wind','radians') + CS%id_uStk0 = register_diag_field('ocean_model', 'uStk0' , diag%axesT1, & + Time, 'Zonal Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vStk0 = register_diag_field('ocean_model', 'vStk0' , diag%axesT1, & + Time, 'Merid Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s) + endif CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) diff --git a/src/tracer/MARBL_tracers.F90 b/src/tracer/MARBL_tracers.F90 index baf7931e51..0896917f2c 100644 --- a/src/tracer/MARBL_tracers.F90 +++ b/src/tracer/MARBL_tracers.F90 @@ -82,6 +82,8 @@ module MARBL_tracers integer :: alk_alt_co2_ind !< ALK_ALT_CO2 index integer :: dic_ind !< DIC index integer :: dic_alt_co2_ind !< DIC_ALT_CO2 index + integer :: abio_dic_ind !< ABIO_DIC index + integer :: abio_di14c_ind !< ABIO_DI14C index end type tracer_ind_type !> MOM needs to store some information about saved_state; besides providing these @@ -183,6 +185,7 @@ module MARBL_tracers !! because we already copy data into CS%STF; latter requires copying data and indices !! so currently using temp_MARBL_diag for that. integer, allocatable :: id_surface_flux_out(:) !< register_diag indices for surface_flux output + integer, allocatable :: id_surface_flux_from_salt_flux(:) !< register_diag indices for surface_flux from salt_flux type(temp_MARBL_diag), allocatable :: interior_tendency_out(:) !< collect interior tendencies for diagnostic output type(temp_MARBL_diag), allocatable :: interior_tendency_out_zint(:) !< vertical integral of interior tendencies !! (full column) @@ -192,6 +195,9 @@ module MARBL_tracers integer, allocatable :: fracr_cat_id(:) !< register_diag index for per-category ice fraction integer, allocatable :: qsw_cat_id(:) !< register_diag index for per-category shortwave + real :: DIC_salt_ratio !< ratio to convert salt surface flux to DIC surface flux [conc ppt-1] + real :: ALK_salt_ratio !< ratio to convert salt surface flux to ALK surface flux [conc ppt-1] + real, allocatable :: STF(:,:,:) !< surface fluxes returned from MARBL to use in tracer_vertdiff !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1] real, allocatable :: SFO(:,:,:) !< surface flux output returned from MARBL for use in GCM @@ -703,7 +709,14 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS, call log_param(param_file, mdl, "INPUTDIR/D14C_FILE", CS%d14c_dataset(m)%file_name) endif enddo -endif + endif + + call get_param(param_file, mdl, "DIC_SALT_RATIO", CS%DIC_salt_ratio, & + "Ratio to convert salt surface flux to DIC surface flux", units="conc ppt-1", & + default=64.0) + call get_param(param_file, mdl, "ALK_SALT_RATIO", CS%ALK_salt_ratio, & + "Ratio to convert salt surface flux to ALK surface flux", units="conc ppt-1", & + default=70.0) ! ** Tracer Restoring call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_SOURCE", CS%restoring_source, & @@ -858,6 +871,7 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag ! Register per-tracer diagnostics computed from MARBL surface flux / interior tendency values allocate(CS%id_surface_flux_out(CS%ntr)) + allocate(CS%id_surface_flux_from_salt_flux(CS%ntr)) allocate(CS%interior_tendency_out(CS%ntr)) allocate(CS%interior_tendency_out_zint(CS%ntr)) allocate(CS%interior_tendency_out_zint_100m(CS%ntr)) @@ -869,6 +883,12 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag diag%axesT1, & ! T => tracer grid? 1 => no vertical grid day, trim(longname), trim(units), conversion=US%Z_to_m*US%s_to_T) + write(name, "(2A)") "STF_SALT_", trim(MARBL_instances%tracer_metadata(m)%short_name) + write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), " Surface Flux from Salt Flux" + CS%id_surface_flux_from_salt_flux(m) = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T => tracer grid? 1 => no vertical grid + day, trim(longname), trim(units), conversion=US%Z_to_m*US%s_to_T) + write(name, "(2A)") "J_", trim(MARBL_instances%tracer_metadata(m)%short_name) write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), " Source Sink Term" write(units, "(2A)") trim(MARBL_instances%tracer_metadata(m)%units), "/s" @@ -1237,8 +1257,12 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which !! fluxes can be applied [m] -! Local variables + ! Local variables character(len=256) :: log_message + real, dimension(SZI_(G),SZJ_(G)) :: net_salt_rate ! Surface salt flux into the ocean + ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]. + real, dimension(SZI_(G),SZJ_(G)) :: flux_from_salt_flux ! Surface tracer flux from salt flux + ! [conc Z T-1 ~> conc m s-1]. real, dimension(SZI_(G),SZJ_(G)) :: ref_mask ! Mask for 2D MARBL diags using ref_depth real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_loc ! Local copy of CS%RIV_FLUXES*dt real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified @@ -1368,6 +1392,69 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, enddo enddo + ! convert salt flux to tracer fluxes and add to STF + do j=js,je ; do i=is,ie + net_salt_rate(i,j) = (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j)) * GV%RZ_to_H + enddo ; enddo + + ! DIC related tracers + do j=js,je ; do i=is,ie + flux_from_salt_flux(i,j) = (CS%DIC_salt_ratio * GV%H_to_Z) * net_salt_rate(i,j) + enddo ; enddo + m = CS%tracer_inds%dic_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + m = CS%tracer_inds%dic_alt_co2_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + m = CS%tracer_inds%abio_dic_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + m = CS%tracer_inds%abio_di14c_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + + ! ALK related tracers + do j=js,je ; do i=is,ie + flux_from_salt_flux(i,j) = (CS%ALK_salt_ratio * GV%H_to_Z) * net_salt_rate(i,j) + enddo ; enddo + m = CS%tracer_inds%alk_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + m = CS%tracer_inds%alk_alt_co2_ind + if (m > 0) then + do j=js,je ; do i=is,ie + CS%STF(i,j,m) = CS%STF(i,j,m) + flux_from_salt_flux(i,j) + enddo ; enddo + if (CS%id_surface_flux_from_salt_flux(m) > 0) & + call post_data(CS%id_surface_flux_from_salt_flux(m), flux_from_salt_flux, CS%diag) + endif + if (CS%debug) then do m=1,CS%ntr call hchksum(CS%STF(:,:,m), & @@ -1923,7 +2010,7 @@ function MARBL_tracers_stock(h, stocks, G, GV, CS, names, units, stock_index) integer :: MARBL_tracers_stock !< Return value: the number of stocks !! calculated here. -! Local variables + ! Local variables integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -2061,6 +2148,8 @@ subroutine set_riv_flux_tracer_inds(CS) CS%tracer_inds%alk_alt_co2_ind = 0 CS%tracer_inds%dic_ind = 0 CS%tracer_inds%dic_alt_co2_ind = 0 + CS%tracer_inds%abio_dic_ind = 0 + CS%tracer_inds%abio_di14c_ind = 0 do m=1,CS%ntr name = MARBL_instances%tracer_metadata(m)%short_name if (trim(name) == "NO3") then @@ -2091,6 +2180,10 @@ subroutine set_riv_flux_tracer_inds(CS) CS%tracer_inds%dic_ind = m elseif (trim(name) == "DIC_ALT_CO2") then CS%tracer_inds%dic_alt_co2_ind = m + elseif (trim(name) == "ABIO_DIC") then + CS%tracer_inds%abio_dic_ind = m + elseif (trim(name) == "ABIO_DI14C") then + CS%tracer_inds%abio_di14c_ind = m endif enddo diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 656ff5b569..3744469891 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -715,7 +715,7 @@ subroutine Update_Surface_Waves(G, GV, US, Time_present, dt, CS, forces) enddo do j=G%jsc,G%jec do i=G%isc,G%iec - !CS%Omega_w2x(i,j) = forces%omega_w2x(i,j) + CS%Omega_w2x(i,j) = forces%omega_w2x(i,j) do b=1,CS%NumBands CS%UStk_Hb(i,j,b) = forces%UStkb(i,j,b) CS%VStk_Hb(i,j,b) = forces%VStkb(i,j,b)