diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index fb98a7b2bf..6583c21ff0 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -80,16 +80,8 @@ module MOM_ocean_model_mct public ocean_model_restart public ice_ocn_bnd_type_chksum public ocean_public_type_chksum -public ocean_model_data_get public get_ocean_grid -!> This interface extracts a named scalar field or array from the ocean surface or public type -interface ocean_model_data_get - module procedure ocean_model_data1D_get - module procedure ocean_model_data2D_get -end interface - - !> This type is used for communication with other components via the FMS coupler. !! The element names and types can be changed only with great deliberation, hence !! the persistnce of things like the cutsy element name "avg_kount". @@ -1052,98 +1044,6 @@ subroutine Ocean_stock_pe(OS, index, value, time_index) end subroutine Ocean_stock_pe -!> This subroutine extracts a named 2-D field from the ocean surface or public type -subroutine ocean_model_data2D_get(OS,Ocean, name, array2D,isc,jsc) - use MOM_constants, only : CELSIUS_KELVIN_OFFSET - type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the - !! internal ocean state (intent in). - type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly - !! visible ocean surface fields. - character(len=*) , intent(in) :: name !< The name of the field to extract - real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must - !! cover only the computational domain - integer , intent(in) :: isc !< The starting i-index of array2D - integer , intent(in) :: jsc !< The starting j-index of array2D - - integer :: g_isc, g_iec, g_jsc, g_jec,g_isd, g_ied, g_jsd, g_jed, i, j - - if (.not.associated(OS)) return - if (.not.OS%is_ocean_pe) return - -! The problem is %areaT is on MOM domain but Ice_Ocean_Boundary%... is on mpp domain. -! We want to return the MOM data on the mpp (compute) domain -! Get MOM domain extents - call mpp_get_compute_domain(OS%grid%Domain%mpp_domain, g_isc, g_iec, g_jsc, g_jec) - call mpp_get_data_domain (OS%grid%Domain%mpp_domain, g_isd, g_ied, g_jsd, g_jed) - - g_isc = g_isc-g_isd+1 ; g_iec = g_iec-g_isd+1 ; g_jsc = g_jsc-g_jsd+1 ; g_jec = g_jec-g_jsd+1 - - - select case(name) - case('area') - array2D(isc:,jsc:) = OS%US%L_to_m**2*OS%grid%areaT(g_isc:g_iec,g_jsc:g_jec) - case('mask') - array2D(isc:,jsc:) = OS%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec) -!OR same result -! do j=g_jsc,g_jec ; do i=g_isc,g_iec -! array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j) -! enddo ; enddo - case('t_surf') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_pme') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_runoff') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_calving') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('btfHeat') - array2D(isc:,jsc:) = 0 - case('tlat') - array2D(isc:,jsc:) = OS%grid%geoLatT(g_isc:g_iec,g_jsc:g_jec) - case('tlon') - array2D(isc:,jsc:) = OS%grid%geoLonT(g_isc:g_iec,g_jsc:g_jec) - case('ulat') - array2D(isc:,jsc:) = OS%grid%geoLatCu(g_isc:g_iec,g_jsc:g_jec) - case('ulon') - array2D(isc:,jsc:) = OS%grid%geoLonCu(g_isc:g_iec,g_jsc:g_jec) - case('vlat') - array2D(isc:,jsc:) = OS%grid%geoLatCv(g_isc:g_iec,g_jsc:g_jec) - case('vlon') - array2D(isc:,jsc:) = OS%grid%geoLonCv(g_isc:g_iec,g_jsc:g_jec) - case('geoLatBu') - array2D(isc:,jsc:) = OS%grid%geoLatBu(g_isc:g_iec,g_jsc:g_jec) - case('geoLonBu') - array2D(isc:,jsc:) = OS%grid%geoLonBu(g_isc:g_iec,g_jsc:g_jec) - case('cos_rot') - array2D(isc:,jsc:) = OS%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1 - case('sin_rot') - array2D(isc:,jsc:) = OS%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0 - case default - call MOM_error(FATAL,'get_ocean_grid_data2D: unknown argument name='//name) - end select -end subroutine ocean_model_data2D_get - -!> This subroutine extracts a named scalar field from the ocean surface or public type -subroutine ocean_model_data1D_get(OS, Ocean, name, value) - type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the - !! internal ocean state (intent in). - type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly - !! visible ocean surface fields. - character(len=*) , intent(in) :: name !< The name of the field to extract - real , intent(out):: value !< The value of the named field - - if (.not.associated(OS)) return - if (.not.OS%is_ocean_pe) return - - select case(name) - case('c_p') - value = OS%C_p - case default - call MOM_error(FATAL,'get_ocean_grid_data1D: unknown argument name='//name) - end select - -end subroutine ocean_model_data1D_get - !> Write out FMS-format checsums on fields from the ocean surface state subroutine ocean_public_type_chksum(id, timestep, ocn) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index a6aebade08..b1ce9a60c0 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -36,7 +36,7 @@ module ocn_comp_mct use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) use MOM_time_manager, only: operator(==), operator(/=), operator(>), get_time -use MOM_file_parser, only: get_param, log_version, param_file_type +use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file use MOM_get_input, only: Get_MOM_Input, directories use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct use MOM_constants, only: CELSIUS_KELVIN_OFFSET @@ -281,6 +281,9 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) glb%c1 = 0.0; glb%c2 = 0.0; glb%c3 = 0.0; glb%c4 = 0.0 endif + ! Close param file before it gets opened by ocean_model_init again. + call close_param_file(param_file) + ! Initialize the MOM6 model runtype = get_runtype() if (runtype == "initial") then diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index cbbefd6afd..219245e473 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -123,6 +123,7 @@ module MOM_cap_mod integer :: export_slice = 1 character(len=256) :: tmpstr logical :: write_diagnostics = .false. +logical :: overwrite_timeslice = .false. character(len=32) :: runtype !< run type integer :: logunit !< stdout logging unit number logical :: profile_memory = .true. @@ -278,6 +279,21 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return + overwrite_timeslice = .false. + call NUOPC_CompAttributeGet(gcomp, name="OverwriteSlice", value=value, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return + if (isPresent .and. isSet) overwrite_timeslice=(trim(value)=="true") + write(logmsg,*) overwrite_timeslice + call ESMF_LogWrite('MOM_cap:OverwriteSlice = '//trim(logmsg), ESMF_LOGMSG_INFO, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return + profile_memory = .false. call NUOPC_CompAttributeGet(gcomp, name="ProfileMemory", value=value, & isPresent=isPresent, isSet=isSet, rc=rc) @@ -708,12 +724,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), & Ice_ocean_boundary% mi (isc:iec,jsc:jec), & Ice_ocean_boundary% p (isc:iec,jsc:jec), & - Ice_ocean_boundary% runoff (isc:iec,jsc:jec), & - Ice_ocean_boundary% calving (isc:iec,jsc:jec), & - Ice_ocean_boundary% runoff_hflx (isc:iec,jsc:jec), & - Ice_ocean_boundary% calving_hflx (isc:iec,jsc:jec), & - Ice_ocean_boundary% rofl_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% rofi_flux (isc:iec,jsc:jec)) + Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec), & + Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec), & + Ice_ocean_boundary% lrunoff (isc:iec,jsc:jec), & + Ice_ocean_boundary% frunoff (isc:iec,jsc:jec)) Ice_ocean_boundary%u_flux = 0.0 Ice_ocean_boundary%v_flux = 0.0 @@ -731,12 +745,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%seaice_melt_heat= 0.0 Ice_ocean_boundary%mi = 0.0 Ice_ocean_boundary%p = 0.0 - Ice_ocean_boundary%runoff = 0.0 - Ice_ocean_boundary%calving = 0.0 - Ice_ocean_boundary%runoff_hflx = 0.0 - Ice_ocean_boundary%calving_hflx = 0.0 - Ice_ocean_boundary%rofl_flux = 0.0 - Ice_ocean_boundary%rofi_flux = 0.0 + Ice_ocean_boundary%lrunoff_hflx = 0.0 + Ice_ocean_boundary%frunoff_hflx = 0.0 + Ice_ocean_boundary%lrunoff = 0.0 + Ice_ocean_boundary%frunoff = 0.0 ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state call ESMF_GridCompSetInternalState(gcomp, ocean_internalstate, rc) @@ -745,11 +757,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return ! bail out + if (len_trim(scalar_field_name) > 0) then + call fld_list_add(fldsToOcn_num, fldsToOcn, trim(scalar_field_name), "will_provide") + call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide") + end if + if (cesm_coupled) then - if (len_trim(scalar_field_name) > 0) then - call fld_list_add(fldsToOcn_num, fldsToOcn, trim(scalar_field_name), "will_provide") - call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide") - endif !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes" , "will provide") @@ -780,11 +793,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi" , "will provide") !-> ice runoff call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn" , "will provide") - - !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_rate" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_rate" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") + !These are not currently used and changing requires a nuopc dictionary change + !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") + !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") !--------- export fields ------------- call fld_list_add(fldsFrOcn_num, fldsFrOcn, "ocean_mask" , "will provide") @@ -1494,13 +1505,11 @@ subroutine DataInitialize(gcomp, rc) ocean_state => ocean_internalstate%ptr%ocean_state_type_ptr call get_ocean_grid(ocean_state, ocean_grid) - if (cesm_coupled) then - call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - endif + call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out call ESMF_StateGet(exportState, itemCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -1543,7 +1552,7 @@ subroutine DataInitialize(gcomp, rc) if(write_diagnostics) then call NUOPC_Write(exportState, fileNamePrefix='field_init_ocn_export_', & - timeslice=import_slice, relaxedFlag=.true., rc=rc) + overwrite=overwrite_timeslice,timeslice=import_slice, relaxedFlag=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -1697,7 +1706,7 @@ subroutine ModelAdvance(gcomp, rc) if (write_diagnostics) then call NUOPC_Write(importState, fileNamePrefix='field_ocn_import_', & - timeslice=import_slice, relaxedFlag=.true., rc=rc) + overwrite=overwrite_timeslice,timeslice=import_slice, relaxedFlag=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -1853,7 +1862,7 @@ subroutine ModelAdvance(gcomp, rc) if (write_diagnostics) then call NUOPC_Write(exportState, fileNamePrefix='field_ocn_export_', & - timeslice=export_slice, relaxedFlag=.true., rc=rc) + overwrite=overwrite_timeslice,timeslice=export_slice, relaxedFlag=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index 2f872c7da5..70915d0e95 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -214,64 +214,42 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, return ! bail out !---- - ! runoff and heat content of runoff + ! mass and heat content of liquid and frozen runoff !---- ! Note - preset values to 0, if field does not exist in importState, then will simply return ! and preset value will be used ! liquid runoff - ice_ocean_boundary%rofl_flux (:,:) = 0._ESMF_KIND_R8 + ice_ocean_boundary%lrunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofl', & - isc, iec, jsc, jec, ice_ocean_boundary%rofl_flux,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out ! ice runoff - ice_ocean_boundary%rofi_flux (:,:) = 0._ESMF_KIND_R8 + ice_ocean_boundary%frunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofi', & - isc, iec, jsc, jec, ice_ocean_boundary%rofi_flux,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%frunoff,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! total runoff - ice_ocean_boundary%runoff (:,:) = 0._ESMF_KIND_R8 - call state_getimport(importState, 'mean_runoff_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%runoff, rc=rc) + ! heat content of lrunoff + ice_ocean_boundary%lrunoff_hflx(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'mean_runoff_heat_flx', & + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! heat content of runoff - ice_ocean_boundary%runoff_hflx(:,:) = 0._ESMF_KIND_R8 - call state_getimport(importState, 'mean_runoff_heat_flux', & - isc, iec, jsc, jec, ice_ocean_boundary%runoff_hflx, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - - !---- - ! calving rate and heat flux - !---- - ! Note - preset values to 0, if field does not exist in importState, then will simply return - ! and preset value will be used - - ice_ocean_boundary%calving(:,:) = 0._ESMF_KIND_R8 - call state_getimport(importState, 'mean_calving_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%calving, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - - ice_ocean_boundary%calving_hflx(:,:) = 0._ESMF_KIND_R8 - call state_getimport(importState, 'mean_calving_heat_flux', & - isc, iec, jsc, jec, ice_ocean_boundary%calving_hflx, rc=rc) + ! heat content of frunoff + ice_ocean_boundary%frunoff_hflx(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'mean_calving_heat_flx', & + isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 240b576669..202fa5791d 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -77,17 +77,9 @@ module MOM_ocean_model_nuopc public ocean_model_restart public ice_ocn_bnd_type_chksum public ocean_public_type_chksum -public ocean_model_data_get public get_ocean_grid public get_eps_omesh -!> This interface extracts a named scalar field or array from the ocean surface or public type -interface ocean_model_data_get - module procedure ocean_model_data1D_get - module procedure ocean_model_data2D_get -end interface - - !> This type is used for communication with other components via the FMS coupler. !! The element names and types can be changed only with great deliberation, hence !! the persistnce of things like the cutsy element name "avg_kount". @@ -1047,99 +1039,6 @@ subroutine Ocean_stock_pe(OS, index, value, time_index) end subroutine Ocean_stock_pe -!> This subroutine extracts a named 2-D field from the ocean surface or public type -subroutine ocean_model_data2D_get(OS, Ocean, name, array2D, isc, jsc) - use MOM_constants, only : CELSIUS_KELVIN_OFFSET - type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the - !! internal ocean state (intent in). - type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly - !! visible ocean surface fields. - character(len=*) , intent(in) :: name !< The name of the field to extract - real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must - !! cover only the computational domain - integer , intent(in) :: isc !< The starting i-index of array2D - integer , intent(in) :: jsc !< The starting j-index of array2D - - integer :: g_isc, g_iec, g_jsc, g_jec,g_isd, g_ied, g_jsd, g_jed, i, j - - if (.not.associated(OS)) return - if (.not.OS%is_ocean_pe) return - -! The problem is %areaT is on MOM domain but Ice_Ocean_Boundary%... is on mpp domain. -! We want to return the MOM data on the mpp (compute) domain -! Get MOM domain extents - call mpp_get_compute_domain(OS%grid%Domain%mpp_domain, g_isc, g_iec, g_jsc, g_jec) - call mpp_get_data_domain (OS%grid%Domain%mpp_domain, g_isd, g_ied, g_jsd, g_jed) - - g_isc = g_isc-g_isd+1 ; g_iec = g_iec-g_isd+1 ; g_jsc = g_jsc-g_jsd+1 ; g_jec = g_jec-g_jsd+1 - - - select case(name) - case('area') - array2D(isc:,jsc:) = OS%US%L_to_m**2*OS%grid%areaT(g_isc:g_iec,g_jsc:g_jec) - case('mask') - array2D(isc:,jsc:) = OS%grid%mask2dT(g_isc:g_iec,g_jsc:g_jec) -!OR same result -! do j=g_jsc,g_jec ; do i=g_isc,g_iec -! array2D(isc+i-g_isc,jsc+j-g_jsc) = OS%grid%mask2dT(i,j) -! enddo ; enddo - case('t_surf') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_pme') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_runoff') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('t_calving') - array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET - case('btfHeat') - array2D(isc:,jsc:) = 0 - case('tlat') - array2D(isc:,jsc:) = OS%grid%geoLatT(g_isc:g_iec,g_jsc:g_jec) - case('tlon') - array2D(isc:,jsc:) = OS%grid%geoLonT(g_isc:g_iec,g_jsc:g_jec) - case('ulat') - array2D(isc:,jsc:) = OS%grid%geoLatCu(g_isc:g_iec,g_jsc:g_jec) - case('ulon') - array2D(isc:,jsc:) = OS%grid%geoLonCu(g_isc:g_iec,g_jsc:g_jec) - case('vlat') - array2D(isc:,jsc:) = OS%grid%geoLatCv(g_isc:g_iec,g_jsc:g_jec) - case('vlon') - array2D(isc:,jsc:) = OS%grid%geoLonCv(g_isc:g_iec,g_jsc:g_jec) - case('geoLatBu') - array2D(isc:,jsc:) = OS%grid%geoLatBu(g_isc:g_iec,g_jsc:g_jec) - case('geoLonBu') - array2D(isc:,jsc:) = OS%grid%geoLonBu(g_isc:g_iec,g_jsc:g_jec) - case('cos_rot') - array2D(isc:,jsc:) = OS%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1 - case('sin_rot') - array2D(isc:,jsc:) = OS%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0 - case default - call MOM_error(FATAL,'get_ocean_grid_data2D: unknown argument name='//name) - end select - -end subroutine ocean_model_data2D_get - -!> This subroutine extracts a named scalar field from the ocean surface or public type -subroutine ocean_model_data1D_get(OS, Ocean, name, value) - type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the - !! internal ocean state (intent in). - type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly - !! visible ocean surface fields. - character(len=*) , intent(in) :: name !< The name of the field to extract - real , intent(out):: value !< The value of the named field - - if (.not.associated(OS)) return - if (.not.OS%is_ocean_pe) return - - select case(name) - case('c_p') - value = OS%C_p - case default - call MOM_error(FATAL,'get_ocean_grid_data1D: unknown argument name='//name) - end select - -end subroutine ocean_model_data1D_get - !> Write out FMS-format checsums on fields from the ocean surface state subroutine ocean_public_type_chksum(id, timestep, ocn) diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index 270d4e9f4c..b75ff96532 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -107,6 +107,7 @@ module MOM_surface_forcing_nuopc !! sea-ice viscosity becomes effective, in kg m-2, !! typically of order 1000 [kg m-2]. logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments + logical :: liquid_runoff_from_data !< If true, use data_override to obtain liquid runoff real :: Flux_const !< piston velocity for surface restoring [m/s] logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux @@ -152,8 +153,8 @@ module MOM_surface_forcing_nuopc !> Structure corresponding to forcing, but with the elements, units, and conventions !! that exactly conform to the use for MOM-based coupled models. type, public :: ice_ocean_boundary_type - real, pointer, dimension(:,:) :: rofl_flux =>NULL() !< liquid runoff [kg/m2/s] - real, pointer, dimension(:,:) :: rofi_flux =>NULL() !< ice runoff [kg/m2/s] + real, pointer, dimension(:,:) :: lrunoff =>NULL() !< liquid runoff [kg/m2/s] + real, pointer, dimension(:,:) :: frunoff =>NULL() !< ice runoff [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] @@ -168,13 +169,11 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() !< diffuse Near InfraRed sw radiation [W/m2] real, pointer, dimension(:,:) :: lprec =>NULL() !< mass flux of liquid precip [kg/m2/s] real, pointer, dimension(:,:) :: fprec =>NULL() !< mass flux of frozen precip [kg/m2/s] - real, pointer, dimension(:,:) :: runoff =>NULL() !< mass flux of liquid runoff [kg/m2/s] - real, pointer, dimension(:,:) :: calving =>NULL() !< mass flux of frozen runoff [kg/m2/s] real, pointer, dimension(:,:) :: ustar_berg =>NULL() !< frictional velocity beneath icebergs [m/s] real, pointer, dimension(:,:) :: area_berg =>NULL() !< area covered by icebergs[m2/m2] real, pointer, dimension(:,:) :: mass_berg =>NULL() !< mass of icebergs(kg/m2) - real, pointer, dimension(:,:) :: runoff_hflx =>NULL() !< heat content of liquid runoff [W/m2] - real, pointer, dimension(:,:) :: calving_hflx =>NULL() !< heat content of frozen runoff [W/m2] + real, pointer, dimension(:,:) :: lrunoff_hflx =>NULL() !< heat content of liquid runoff [W/m2] + real, pointer, dimension(:,:) :: frunoff_hflx =>NULL() !< heat content of frozen runoff [W/m2] real, pointer, dimension(:,:) :: p =>NULL() !< pressure of overlying ice and atmosphere !< on ocean surface [Pa] real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] @@ -422,6 +421,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo endif + ! Check that liquid runoff has a place to go + if (CS%liquid_runoff_from_data .and. .not. associated(IOB%lrunoff)) then + call MOM_error(FATAL, "liquid runoff is being added via data_override but "// & + "there is no associated runoff in the IOB") + return + end if + ! obtain fluxes from IOB; note the staggering of indices i0 = is - isc_bnd ; j0 = js - jsc_bnd do j=js,je ; do i=is,ie @@ -436,17 +442,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%evap(i,j) = kg_m2_s_conversion * IOB%q_flux(i-i0,j-j0) * G%mask2dT(i,j) ! liquid runoff flux - if (associated(IOB%rofl_flux)) then - fluxes%lrunoff(i,j) = kg_m2_s_conversion * IOB%rofl_flux(i-i0,j-j0) * G%mask2dT(i,j) - else if (associated(IOB%runoff)) then - fluxes%lrunoff(i,j) = kg_m2_s_conversion * IOB%runoff(i-i0,j-j0) * G%mask2dT(i,j) + if (associated(IOB%lrunoff)) then + if(CS%liquid_runoff_from_data)call data_override('OCN', 'runoff', IOB%lrunoff, Time) + fluxes%lrunoff(i,j) = kg_m2_s_conversion * IOB%lrunoff(i-i0,j-j0) * G%mask2dT(i,j) endif ! ice runoff flux - if (associated(IOB%rofi_flux)) then - fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%rofi_flux(i-i0,j-j0) * G%mask2dT(i,j) - elseif (associated(IOB%calving)) then - fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%calving(i-i0,j-j0) * G%mask2dT(i,j) + if (associated(IOB%frunoff)) then + fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff(i-i0,j-j0) * G%mask2dT(i,j) endif if (associated(IOB%ustar_berg)) & @@ -458,11 +461,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%mass_berg)) & fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) - if (associated(IOB%runoff_hflx)) & - fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) + if (associated(IOB%lrunoff_hflx)) & + fluxes%heat_content_lrunoff(i,j) = IOB%lrunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) - if (associated(IOB%calving_hflx)) & - fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion*IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j) + if (associated(IOB%frunoff_hflx)) & + fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * IOB%frunoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) if (associated(IOB%lw_flux)) & fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j) @@ -483,9 +486,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion endif - if (associated(IOB%calving)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + if (associated(IOB%frunoff)) then + fluxes%latent(i,j) = fluxes%latent(i,j) + IOB%frunoff(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%frunoff(i-i0,j-j0)*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)*CS%latent_heat_vapor @@ -1284,7 +1287,12 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & "If true, allows flux adjustments to specified via the "//& "data_table using the component name 'OCN'.", default=.false.) - if (CS%allow_flux_adjustments) then + + call get_param(param_file, mdl, "LIQUID_RUNOFF_FROM_DATA", CS%liquid_runoff_from_data, & + "If true, allows liquid river runoff to be specified via the "//& + "data_table using the component name 'OCN'.", default=.false.) + + if (CS%allow_flux_adjustments .or. CS%liquid_runoff_from_data) then call data_override_init(Ocean_domain_in=G%Domain%mpp_domain) endif @@ -1374,8 +1382,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) write(outunit,100) 'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif) write(outunit,100) 'iobt%lprec ' , mpp_chksum( iobt%lprec ) write(outunit,100) 'iobt%fprec ' , mpp_chksum( iobt%fprec ) - write(outunit,100) 'iobt%runoff ' , mpp_chksum( iobt%runoff ) - write(outunit,100) 'iobt%calving ' , mpp_chksum( iobt%calving ) + write(outunit,100) 'iobt%lrunoff ' , mpp_chksum( iobt%lrunoff ) + write(outunit,100) 'iobt%frunoff ' , mpp_chksum( iobt%frunoff ) write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p ) if (associated(iobt%ustar_berg)) & write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 690e5250db..0926867cce 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -77,7 +77,7 @@ module MOM use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_interface_heights, only : find_eta use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init -use MOM_lateral_mixing_coeffs, only : calc_resoln_function, VarMix_CS +use MOM_lateral_mixing_coeffs, only : calc_resoln_function, calc_depth_function, VarMix_CS use MOM_MEKE, only : MEKE_init, MEKE_alloc_register_restart, step_forward_MEKE, MEKE_CS use MOM_MEKE_types, only : MEKE_type use MOM_mixed_layer_restrat, only : mixedlayer_restrat, mixedlayer_restrat_init, mixedlayer_restrat_CS @@ -566,6 +566,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & if (associated(CS%VarMix)) then call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag) call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix) + call calc_depth_function(G, CS%VarMix) call disable_averaging(CS%diag) endif endif @@ -1406,6 +1407,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (associated(CS%VarMix)) then call pass_var(CS%h, G%Domain) call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix) + call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix) endif call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, & @@ -1431,6 +1433,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (associated(CS%VarMix)) then call pass_var(CS%h, G%Domain) call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix) + call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix) endif call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, & @@ -2373,7 +2376,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp) - call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, & + call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & CS%tracer_diff_CSp) call lock_tracer_registry(CS%tracer_Reg) @@ -2894,7 +2897,7 @@ subroutine extract_surface_state(CS, sfc_state) if (allocated(sfc_state%melt_potential)) then - !$OMP parallel do default(shared) + !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, delT) do j=js,je do i=is,ie depth(i) = 0.0 diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8c016b11b0..ca94af2225 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -681,7 +681,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, & - OBC=CS%OBC, BT=CS%barotropic_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") @@ -1149,7 +1149,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param .not. query_initialized(CS%diffv,"diffv",restart_CS)) then call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, & G, GV, US, CS%hor_visc_CSp, & - OBC=CS%OBC, BT=CS%barotropic_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, & + TD=thickness_diffuse_CSp) else if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & (US%m_to_L * US%s_to_T_restart**2 /= US%m_to_L_restart * US%s_to_T**2) ) then diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index ff5a93a62c..4197cfea3f 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -5,10 +5,11 @@ module MOM_unit_tests use MOM_error_handler, only : MOM_error, FATAL, is_root_pe -use MOM_string_functions, only : string_functions_unit_tests -use MOM_remapping, only : remapping_unit_tests -use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests -use MOM_diag_vkernels, only : diag_vkernels_unit_tests +use MOM_string_functions, only : string_functions_unit_tests +use MOM_remapping, only : remapping_unit_tests +use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests +use MOM_diag_vkernels, only : diag_vkernels_unit_tests +use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests implicit none ; private @@ -35,6 +36,9 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: neutralDiffusionUnitTests FAILED") if (diag_vkernels_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: diag_vkernels_unit_tests FAILED") + if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: near_boundary_unit_tests FAILED") + endif end subroutine unit_tests diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 95c3ad6916..77b36f85db 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -600,7 +600,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_rhopot2 > 0) call post_data(CS%id_rhopot2, Rcv, CS%diag) endif if (CS%id_rhoinsitu > 0) then -!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d,h,GV) +!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,h,GV) private(pressure_1d) do j=js,je pressure_1d(:) = 0. ! Start at p=0 Pa at surface do k=1,nz diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index a2257369a8..1e785fa930 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -30,6 +30,8 @@ module MOM_MEKE !> Control structure that contains MEKE parameters and diagnostics handles type, public :: MEKE_CS ; private ! Parameters + real, dimension(:,:), pointer :: equilibrium_value => NULL() !< The equilbrium value + !! of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2] real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim] real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim] real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim] @@ -43,8 +45,12 @@ module MOM_MEKE logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag. logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC !! framework (Marshall et al., 2012) + real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of the + !! GEOMETRIC thickness diffusion. logical :: MEKE_equilibrium_alt !< If true, use an alternative calculation for the !! equilibrium value of MEKE. + logical :: MEKE_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value, + !! which is calculated at each time step. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather !! than the streamfunction for the MEKE GM source term. logical :: Rd_as_max_scale !< If true the length scale can not exceed the @@ -75,6 +81,8 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. logical :: debug !< If true, write out checksums of data for debugging @@ -87,6 +95,7 @@ module MOM_MEKE integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1 integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1 integer :: id_Lrhines = -1, id_Leady = -1 + integer :: id_MEKE_equilibrium = -1 !!@} ! Infrastructure @@ -277,7 +286,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Calculates bottomFac2, barotrFac2 and LmixScale call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, bottomFac2, barotrFac2, LmixScale) if (CS%debug) then - call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, scale=US%Z_to_m*US%s_to_T) + if (CS%visc_drag) & + call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, scale=US%Z_to_m*US%s_to_T) call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, scale=US%R_to_kg_m3*US%Z_to_m) call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, scale=US%L_T_to_m_s) call hchksum(bottomFac2, 'MEKE bottomFac2', G%HI) @@ -320,6 +330,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif endif + if (CS%MEKE_equilibrium_restoring) then + call MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) + do j=js,je ; do i=is,ie + src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - CS%equilibrium_value(i,j)) + enddo ; enddo + endif + ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -651,6 +668,8 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2]. logical :: useSecant, debugIteration + real :: Lgrid, Ldeform, Lfrict + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec debugIteration = .false. @@ -665,112 +684,112 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & - (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points - - ! Since zero-bathymetry cells are masked, this avoids calculations on land - if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.) then - beta_topo_x = 0. ; beta_topo_y = 0. + if (CS%MEKE_equilibrium_alt) then + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. - beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & - + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) - beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & - (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) - endif - beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & - (G%dF_dy(i,j) + beta_topo_y)**2 ) - - I_H = US%L_to_Z*GV%Rho0 * I_mass(i,j) - - if (KhCoeff*SN*I_H>0.) then - ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E - EKEmin = 0. ! Use the trivial root as the left bracket - ResMin = 0. ! Need to detect direction of left residual - EKEmax = 0.01*US%m_s_to_L_T**2 ! First guess at right bracket - useSecant = .false. ! Start using a bisection method - - ! First find right bracket for which resid<0 - resid = 1.0*US%m_to_L**2*US%T_to_s**3 ; n1 = 0 - do while (resid>0.) - n1 = n1 + 1 - EKE = EKEmax - call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & - MEKE%Rd_dx_h(i,j), SN, EKE, & - bottomFac2, barotrFac2, LmixScale, LRhines, LEady) - ! TODO: Should include resolution function in Kh - Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) - src = Kh * (SN * SN) - drag_rate = I_H * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) - ldamping = CS%MEKE_damping + drag_rate * bottomFac2 - resid = src - ldamping * EKE - ! if (debugIteration) then - ! write(0,*) n1, 'EKE=',EKE,'resid=',resid - ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin - ! write(0,*) 'src=',src,'ldamping=',ldamping - ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2 - ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2 - ! endif - if (resid>0.) then ! EKE is to the left of the root - EKEmin = EKE ! so we move the left bracket here - EKEmax = 10. * EKE ! and guess again for the right bracket - if (resid 2.e17*US%m_s_to_L_T**2) then - if (debugIteration) stop 'Something has gone very wrong' - debugIteration = .true. - resid = 1. ; n1 = 0 - EKEmin = 0. ; ResMin = 0. - EKEmax = 0.01*US%m_s_to_L_T**2 - useSecant = .false. + FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points + + ! Since zero-bathymetry cells are masked, this avoids calculations on land + if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.) then + beta_topo_x = 0. ; beta_topo_y = 0. + else + !### Consider different combinations of these estimates of topographic beta, and the use + ! of the water column thickness instead of the bathymetric depth. + beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & + (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & + / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & + (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & + / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & + / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + endif + beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & + (G%dF_dy(i,j) + beta_topo_y)**2 ) + + I_H = US%L_to_Z*GV%Rho0 * I_mass(i,j) + + if (KhCoeff*SN*I_H>0.) then + ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E + EKEmin = 0. ! Use the trivial root as the left bracket + ResMin = 0. ! Need to detect direction of left residual + EKEmax = 0.01*US%m_s_to_L_T**2 ! First guess at right bracket + useSecant = .false. ! Start using a bisection method + + ! First find right bracket for which resid<0 + resid = 1.0*US%m_to_L**2*US%T_to_s**3 ; n1 = 0 + do while (resid>0.) + n1 = n1 + 1 + EKE = EKEmax + call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & + MEKE%Rd_dx_h(i,j), SN, EKE, & + bottomFac2, barotrFac2, LmixScale, LRhines, LEady) + ! TODO: Should include resolution function in Kh + Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) + src = Kh * (SN * SN) + drag_rate = I_H * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) + ldamping = CS%MEKE_damping + drag_rate * bottomFac2 + resid = src - ldamping * EKE + ! if (debugIteration) then + ! write(0,*) n1, 'EKE=',EKE,'resid=',resid + ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin + ! write(0,*) 'src=',src,'ldamping=',ldamping + ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2 + ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2 + ! endif + if (resid>0.) then ! EKE is to the left of the root + EKEmin = EKE ! so we move the left bracket here + EKEmax = 10. * EKE ! and guess again for the right bracket + if (resid 2.e17*US%m_s_to_L_T**2) then + if (debugIteration) stop 'Something has gone very wrong' + debugIteration = .true. + resid = 1. ; n1 = 0 + EKEmin = 0. ; ResMin = 0. + EKEmax = 0.01*US%m_s_to_L_T**2 + useSecant = .false. + endif endif - endif - enddo ! while(resid>0.) searching for right bracket - ResMax = resid - - ! Bisect the bracket - n2 = 0 ; EKEerr = EKEmax - EKEmin - do while (EKEerr > tolerance) - n2 = n2 + 1 - if (useSecant) then - EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax)) - else - EKE = 0.5 * (EKEmin + EKEmax) - endif - EKEerr = min( EKE-EKEmin, EKEmax-EKE ) - ! TODO: Should include resolution function in Kh - Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) - src = Kh * (SN * SN) - drag_rate = I_H * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) - ldamping = CS%MEKE_damping + drag_rate * bottomFac2 - resid = src - ldamping * EKE - if (useSecant .and. resid>ResMin) useSecant = .false. - if (resid>0.) then ! EKE is to the left of the root - EKEmin = EKE ! so we move the left bracket here - if (resid EKE is exactly at the root - endif - if (n2>200) stop 'Failing to converge?' - enddo ! while(EKEmax-EKEmin>tolerance) + enddo ! while(resid>0.) searching for right bracket + ResMax = resid + + ! Bisect the bracket + n2 = 0 ; EKEerr = EKEmax - EKEmin + do while (EKEerr > tolerance) + n2 = n2 + 1 + if (useSecant) then + EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax)) + else + EKE = 0.5 * (EKEmin + EKEmax) + endif + EKEerr = min( EKE-EKEmin, EKEmax-EKE ) + ! TODO: Should include resolution function in Kh + Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) + src = Kh * (SN * SN) + drag_rate = I_H * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) + ldamping = CS%MEKE_damping + drag_rate * bottomFac2 + resid = src - ldamping * EKE + if (useSecant .and. resid>ResMin) useSecant = .false. + if (resid>0.) then ! EKE is to the left of the root + EKEmin = EKE ! so we move the left bracket here + if (resid EKE is exactly at the root + endif + if (n2>200) stop 'Failing to converge?' + enddo ! while(EKEmax-EKEmin>tolerance) - else - EKE = 0. - endif - if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (US%Z_to_L*G%bathyT(i,j) * SN / (8*CS%cdrag))**2 - else + else + EKE = 0. + endif MEKE%MEKE(i,j) = EKE endif enddo ; enddo @@ -778,6 +797,37 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m end subroutine MEKE_equilibrium +!< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into +!! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value +subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type. + type(MEKE_CS), pointer :: CS !< MEKE control structure. + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1]. + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1]. + ! Local variables + real :: SN ! The local Eady growth rate [T-1 ~> s-1] + integer :: i, j, is, ie, js, je ! local indices + real :: cd2 ! bottom drag + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + cd2 = CS%cdrag**2 + + if (.not. associated(CS%equilibrium_value)) allocate(CS%equilibrium_value(SZI_(G),SZJ_(G))) + CS%equilibrium_value(:,:) = 0.0 + +!$OMP do + do j=js,je ; do i=is,ie + ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) + ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v + SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + enddo ; enddo + + if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) + +end subroutine MEKE_equilibrium_restoring + !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$ !! functions that are ratios of either bottom or barotropic eddy energy to the !! column eddy energy, respectively. See \ref section_MEKE_equations. @@ -942,6 +992,8 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) ! run to the representation in a restart file. real :: L_rescale ! A rescaling factor for length from the internal representation in this ! run to the representation in a restart file. + real :: MEKE_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value. + real :: cdrag ! The default bottom drag coefficient [nondim]. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed logical :: laplacian, biharmonic, useVarMix, coldStart ! This include declares and sets the variable "version". @@ -1001,9 +1053,22 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//& "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) + call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & + "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//& + "thickness diffusion.", units="nondim", default=0.05) call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", CS%MEKE_equilibrium_alt, & "If true, use an alternative formula for computing the (equilibrium)"//& "initial value of MEKE.", default=.false.) + call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", CS%MEKE_equilibrium_restoring, & + "If true, restore MEKE back to its equilibrium value, which is calculated at"//& + "each time step.", default=.false.) + if (CS%MEKE_equilibrium_restoring) then + call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & + "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & + default=1e6, scale=US%T_to_s) + CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale + endif + call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, & "The efficiency of the conversion of mean energy into "//& "MEKE. If MEKE_FRCOEFF is negative, this conversion "//& @@ -1121,10 +1186,14 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) units="nondim", default=0.0) ! Nonlocal module parameters - call get_param(param_file, mdl, "CDRAG", CS%cdrag, & + call get_param(param_file, mdl, "CDRAG", cdrag, & "CDRAG is the drag coefficient relating the magnitude of "//& "the velocity field to the bottom stress.", units="nondim", & default=0.003) + call get_param(param_file, mdl, "MEKE_CDRAG", CS%cdrag, & + "Drag coefficient relating the magnitude of the velocity "//& + "field to the bottom stress in MEKE.", units="nondim", & + default=cdrag) call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.) @@ -1198,6 +1267,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) 'Meridional diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) endif + if (CS%MEKE_equilibrium_restoring) then + CS%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, Time, & + 'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2) + endif + CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE) ! Detect whether this instance of MEKE_init() is at the beginning of a run diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 63811e14d7..c3ec878bc1 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -75,6 +75,8 @@ module MOM_hor_visc !! Default is False to maintain answers with legacy experiments !! but should be changed to True for new experiments. logical :: anisotropic !< If true, allow anisotropic component to the viscosity. + logical :: add_LES_viscosity!< If true, adds the viscosity from Smagorinsky and Leith to + !! the background viscosity instead of taking the maximum. real :: Kh_aniso !< The anisotropic viscosity [L2 T-1 ~> m2 s-1]. logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function !! of state. This is set depending on ANISOTROPIC_MODE. @@ -201,7 +203,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT) + CS, OBC, BT, TD) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & @@ -226,7 +228,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type type(barotropic_CS), optional, pointer :: BT !< Pointer to a structure containing !! barotropic velocities. - + type(thickness_diffuse_CS), optional, pointer :: TD !< Pointer to a structure containing + !! thickness diffusivities. ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & Del2u, & ! The u-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] @@ -409,15 +412,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call barotropic_get_tav(BT, ubtav, vbtav, G, US) call pass_vector(ubtav, vbtav, G%Domain) - do j=js-1,je+1 ; do i=is-1,ie+1 + do j=js-1,je+2 ; do i=is-1,ie+2 dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - & G%IdyCu(I-1,j) * ubtav(I-1,j)) dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - & G%IdxCv(i,J-1) * vbtav(i,J-1)) enddo; enddo - call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j) enddo ; enddo @@ -430,6 +431,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, - ubtav(I,j)*G%IdxCu(I,j)) enddo ; enddo + call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE) call pass_vector(dvdx_bt, dudy_bt, G%Domain, stagger=AGRID) if (CS%no_slip) then @@ -466,7 +468,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI5, H0_GME, & !$OMP diffu, diffv, max_diss_rate_h, max_diss_rate_q, & !$OMP Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & - !$OMP div_xx_h, vort_xy_q, GME_coeff_h, GME_coeff_q & + !$OMP div_xx_h, vort_xy_q, GME_coeff_h, GME_coeff_q, & + !$OMP TD, KH_u_GME, KH_v_GME & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & @@ -689,8 +692,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif - call pass_var(vort_xy, G%Domain, position=CORNER, complete=.true.) - ! Vorticity gradient do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) @@ -702,23 +703,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo - call pass_vector(vort_xy_dy, vort_xy_dx, G%Domain) - if (CS%modified_Leith) then ! Divergence do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - div_xx(i,j) = 0.5*((G%dyCu(I,j) * u(I,j,k) * (h(i+1,j,k)+h(i,j,k)) - & - G%dyCu(I-1,j) * u(I-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + & - (G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - & - G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k))))*G%IareaT(i,j) / & - (h(i,j,k) + GV%H_subroundoff) + div_xx(i,j) = dudx(i,j) + dvdy(i,j) enddo ; enddo ! Divergence gradient - do j=Jsq,Jeq+1 ; do I=Isq-1,Ieq+1 + do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j)) enddo ; enddo - do J=Jsq-1,Jeq+1 ; do i=Isq,Ieq+1 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j)) enddo ; enddo @@ -727,7 +722,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + !do J=js-1,Jeq ; do I=is-1,Ieq + do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) enddo ; enddo @@ -737,13 +733,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1 div_xx_dx(I,j) = 0.0 enddo ; enddo - do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 div_xx_dy(i,J) = 0.0 enddo ; enddo do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grad_div_mag_h(i,j) = 0.0 enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 grad_div_mag_q(I,J) = 0.0 enddo ; enddo @@ -812,8 +808,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Determine the Laplacian viscosity at h points, using the ! largest value from several parameterizations. Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + if (CS%add_LES_viscosity) then + if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag + if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + else + if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) + if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + endif ! All viscosity contributions above are subject to resolution scaling if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) @@ -837,7 +838,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) -! if (CS%debug) sh_xx_3d(i,j,k) = sh_xx(i,j) str_xx(i,j) = -Kh * sh_xx(i,j) else ! not Laplacian @@ -975,8 +975,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Determine the Laplacian viscosity at q points, using the ! largest value from several parameterizations. Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + if (CS%add_LES_viscosity) then + if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag + if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + else + if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) + if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + endif ! All viscosity contributions above are subject to resolution scaling if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j) @@ -1003,7 +1008,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J) -! if (CS%debug) sh_xy_3d(I,J,k) = sh_xy(I,J) str_xy(I,J) = -Kh * sh_xy(I,J) else ! not Laplacian @@ -1060,48 +1064,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo if (CS%use_GME) then - if (CS%answers_2018) then - do j=js,je ; do i=is,ie - grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + & - (0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2) - max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) - enddo ; enddo - else ! This form is invariant to 90-degree rotations. - do j=js,je ; do i=is,ie - grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + & - ((0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2)) - max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) - enddo ; enddo - endif - - if (CS%answers_2018) then - do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB - grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * (dudx(i,j)**2 + dvdy(i,j)**2 + & - (0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2) - - max_diss_rate_q(I,J,k) = 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)+ & - MEKE%MEKE(i,j+1)+MEKE%MEKE(i+1,j+1)) * sqrt(grad_vel_mag_q(I,J)) - enddo ; enddo - else ! This form is rotationally invariant - do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB - grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * ((dudx(i,j)**2 + dvdy(i,j)**2) + & - ((0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2)) - - max_diss_rate_q(I,J,k) = 0.5*((MEKE%MEKE(i,j) + MEKE%MEKE(i+1,j+1)) + & - (MEKE%MEKE(i+1,j) + MEKE%MEKE(i,j+1))) * sqrt(grad_vel_mag_q(I,J)) - enddo ; enddo - endif + call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G) + call pass_vector(KH_u_GME, KH_v_GME, G%Domain) - do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((grad_vel_mag_bt_h(i,j)>0) .and. (max_diss_rate_h(i,j,k)>0)) then - GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_h(i,j,k) / & - grad_vel_mag_bt_h(i,j) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (grad_vel_mag_bt_h(i,j)>0) then + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) else - GME_coeff = 1.0 + GME_coeff = 0.0 endif ! apply mask @@ -1114,10 +1085,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - - if ((grad_vel_mag_bt_q(I,J)>0) .and. (max_diss_rate_q(I,J,k)>0)) then - GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_q(I,J,k) / & - grad_vel_mag_bt_q(I,J) + if (grad_vel_mag_bt_q(i,j)>0) then + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) else GME_coeff = 0.0 endif @@ -1150,7 +1120,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (associated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie - FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_RZ * grad_vel_mag_bt_h(i,j) + FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * grad_vel_mag_bt_h(i,j) enddo ; enddo endif @@ -1322,8 +1292,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%Laplacian) then call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) -! call Bchksum(sh_xy_3d, "shear_xy", G%HI, haloshift=0, scale=US%s_to_T) -! call hchksum(sh_xx_3d, "shear_xx", G%HI, haloshift=0, scale=US%s_to_T) endif if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) @@ -1394,8 +1362,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) ! valid parameters. logical :: split ! If true, use the split time stepping scheme. ! If false and USE_GME = True, issue a FATAL error. - logical :: use_MEKE ! If true, use the MEKE module for calculating eddy kinetic energy. - ! If false and USE_GME = True, issue a FATAL error. logical :: default_2018_answers character(len=64) :: inputdir, filename @@ -1532,6 +1498,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& "horizontal viscosity.", default=.false.) + call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, & + "If true, adds the viscosity from Smagorinsky and Leith to the "//& + "background viscosity instead of taking the maximum.", default=.false.) endif if (CS%anisotropic .or. get_all) then call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, & @@ -1663,14 +1632,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") - call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & - "If true, turns on the MEKE scheme which calculates\n"// & - "a sub-grid mesoscale eddy kinetic energy budget.", & - default=.false.) - - if (.not. use_MEKE) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & - "cannot be used with USE_MEKE=False.") - call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & "The strength of GME tapers quadratically to zero when the bathymetric "//& "depth is shallower than GME_H0.", units="m", scale=US%m_to_Z, & diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 710012c837..d9543322c9 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -29,6 +29,8 @@ module MOM_lateral_mixing_coeffs !! when the deformation radius is well resolved. logical :: Resoln_scaled_KhTh !< If true, scale away the thickness diffusivity !! when the deformation radius is well resolved. + logical :: Depth_scaled_KhTh !< If true, KHTH is scaled away when the depth is + !! shallower than a reference depth. logical :: Resoln_scaled_KhTr !< If true, scale away the tracer diffusivity !! when the deformation radius is well resolved. logical :: interpolate_Res_fn !< If true, interpolate the resolution function @@ -48,6 +50,8 @@ module MOM_lateral_mixing_coeffs !! This parameter is set depending on other parameters. logical :: calculate_res_fns !< If true, calculate all the resolution factors. !! This parameter is set depending on other parameters. + logical :: calculate_depth_fns !< If true, calculate all the depth factors. + !! This parameter is set depending on other parameters. logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rate. !! This parameter is set depending on other parameters. real, dimension(:,:), pointer :: & @@ -64,6 +68,10 @@ module MOM_lateral_mixing_coeffs !! deformation radius to the grid spacing at u points [nondim]. Res_fn_v => NULL(), & !< Non-dimensional function of the ratio the first baroclinic !! deformation radius to the grid spacing at v points [nondim]. + Depth_fn_u => NULL(), & !< Non-dimensional function of the ratio of the depth to + !! a reference depth (maximum 1) at u points [nondim] + Depth_fn_v => NULL(), & !< Non-dimensional function of the ratio of the depth to + !! a reference depth (maximum 1) at v points [nondim] beta_dx2_h => NULL(), & !< The magnitude of the gradient of the Coriolis parameter !! times the grid spacing squared at h points [L T-1 ~> m s-1]. beta_dx2_q => NULL(), & !< The magnitude of the gradient of the Coriolis parameter @@ -111,6 +119,8 @@ module MOM_lateral_mixing_coeffs real :: Res_coef_visc !< A non-dimensional number that determines the function !! of resolution, used for lateral viscosity, as: !! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power) + real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m] + real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim] real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] integer :: Res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any !! positive integer power may be used, but even powers @@ -140,10 +150,44 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, calc_slope_functions, calc_resoln_function -public calc_QG_Leith_viscosity +public calc_QG_Leith_viscosity, calc_depth_function contains +!> Calculates the non-dimensional depth functions. +subroutine calc_depth_function(G, CS) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(VarMix_CS), pointer :: CS !< Variable mixing coefficients + + ! Local variables + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq + integer :: i, j + real :: H0 ! local variable for reference depth + real :: expo ! exponent used in the depth dependent scaling + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + if (.not. associated(CS)) call MOM_error(FATAL, "calc_depth_function:"// & + "Module must be initialized before it is used.") + if (.not. CS%calculate_depth_fns) return + if (.not. associated(CS%Depth_fn_u)) call MOM_error(FATAL, & + "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.") + if (.not. associated(CS%Depth_fn_v)) call MOM_error(FATAL, & + "calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.") + + H0 = CS%depth_scaled_khth_h0 + expo = CS%depth_scaled_khth_exp +!$OMP do + do j=js,je ; do I=is-1,Ieq + CS%Depth_fn_u(I,j) = (MIN(1.0, 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j))/H0))**expo + enddo ; enddo +!$OMP do + do J=js-1,Jeq ; do i=is,ie + CS%Depth_fn_v(i,J) = (MIN(1.0, 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1))/H0))**expo + enddo ; enddo + +end subroutine calc_depth_function + !> Calculates and stores the non-dimensional resolution functions subroutine calc_resoln_function(h, tv, G, GV, US, CS) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure @@ -904,7 +948,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) CS%calculate_Rd_dx = .false. CS%calculate_res_fns = .false. CS%calculate_Eady_growth_rate = .false. - + CS%calculate_depth_fns = .false. ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "USE_VARIABLE_MIXING", CS%use_variable_mixing,& @@ -920,6 +964,12 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, the Laplacian lateral viscosity is scaled away "//& "when the first baroclinic deformation radius is well "//& "resolved.", default=.false.) + call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", CS%Depth_scaled_KhTh, & + "If true, KHTH is scaled away when the depth is shallower"//& + "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//& + "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//& + "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",& + default=.false.) call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", CS%Resoln_scaled_KhTh, & "If true, the interface depth diffusivity is scaled away "//& "when the first baroclinic deformation radius is well "//& @@ -969,6 +1019,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct) then in_use = .true. call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, & @@ -1152,6 +1203,18 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) endif + if (CS%Depth_scaled_KhTh) then + CS%calculate_depth_fns = .true. + allocate(CS%Depth_fn_u(IsdB:IedB,jsd:jed)) ; CS%Depth_fn_u(:,:) = 0.0 + allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB)) ; CS%Depth_fn_v(:,:) = 0.0 + call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, & + "The depth above which KHTH is scaled away.",& + units="m", default=1000.) + call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, & + "The exponent used in the depth dependent scaling function for KHTH.",& + units="nondim", default=3.0) + endif + ! Resolution %Rd_dx_h CS%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, Time, & 'Ratio between deformation radius and grid spacing', 'm m-1') diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index a567edb4be..53250b7023 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -137,11 +137,12 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real, dimension(SZI_(G), SZJB_(G)) :: & KH_v_CFL ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1] real :: Khth_Loc_u(SZIB_(G), SZJ_(G)) + real :: Khth_Loc_v(SZI_(G), SZJB_(G)) real :: Khth_Loc(SZIB_(G), SZJB_(G)) ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real, dimension(:,:), pointer :: cg1 => null() !< Wave speed [L T-1 ~> m s-1] - logical :: use_VarMix, Resoln_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck + logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck logical :: use_QG_Leith integer :: i, j, k, is, ie, js, je, nz real :: hu(SZI_(G), SZJ_(G)) ! u-thickness [H ~> m or kg m-2] @@ -166,10 +167,12 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp use_VarMix = .false. ; Resoln_scaled = .false. ; use_stored_slopes = .false. khth_use_ebt_struct = .false. ; use_Visbeck = .false. ; use_QG_Leith = .false. + Depth_scaled = .false. if (associated(VarMix)) then use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0.) Resoln_scaled = VarMix%Resoln_scaled_KhTh + Depth_scaled = VarMix%Depth_scaled_KhTh use_stored_slopes = VarMix%use_stored_slopes khth_use_ebt_struct = VarMix%khth_use_ebt_struct use_Visbeck = VarMix%use_Visbeck @@ -198,7 +201,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp !$OMP parallel default(none) shared(is,ie,js,je,Khth_Loc_u,CS,use_VarMix,VarMix, & !$OMP MEKE,Resoln_scaled,KH_u,G,use_QG_Leith,use_Visbeck,& !$OMP KH_u_CFL,nz,Khth_Loc,KH_v,KH_v_CFL,int_slope_u, & -!$OMP int_slope_v,khth_use_ebt_struct) +!$OMP int_slope_v,khth_use_ebt_struct, Depth_scaled, & +!$OMP Khth_loc_v) !$OMP do do j=js,je; do I=is-1,ie Khth_loc_u(I,j) = CS%Khth @@ -236,6 +240,13 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo ; enddo endif + if (Depth_scaled) then +!$OMP do + do j=js,je; do I=is-1,ie + Khth_loc_u(I,j) = Khth_loc_u(I,j) * VarMix%Depth_fn_u(I,j) + enddo ; enddo + endif + if (CS%Khth_Max > 0) then !$OMP do do j=js,je; do I=is-1,ie @@ -282,55 +293,62 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = CS%Khth + Khth_loc_v(i,J) = CS%Khth enddo ; enddo if (use_VarMix) then if (use_Visbeck) then !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = Khth_loc(i,j) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) + Khth_loc_v(i,J) = Khth_loc_v(i,J) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) enddo ; enddo endif endif if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then if (CS%MEKE_GEOMETRIC) then !$OMP do - do j=js-1,je ; do I=is,ie - Khth_loc(I,j) = Khth_loc(I,j) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * & + do J=js-1,je ; do i=is,ie + Khth_loc_v(i,J) = Khth_loc_v(i,J) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * & 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / & (VarMix%SN_v(i,J) + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo else do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = Khth_loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) + Khth_loc_v(i,J) = Khth_loc_v(i,J) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) enddo ; enddo endif endif ; endif if (Resoln_scaled) then +!$OMP do + do J=js-1,je; do i=is,ie + Khth_loc_v(i,J) = Khth_loc_v(i,J) * VarMix%Res_fn_v(i,J) + enddo ; enddo + endif + + if (Depth_scaled) then !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = Khth_loc(i,j) * VarMix%Res_fn_v(i,J) + Khth_loc_v(i,J) = Khth_loc_v(i,J) * VarMix%Depth_fn_v(i,J) enddo ; enddo endif if (CS%Khth_Max > 0) then !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = max(CS%Khth_Min, min(Khth_loc(i,j), CS%Khth_Max)) + Khth_loc_v(i,J) = max(CS%Khth_Min, min(Khth_loc_v(i,J), CS%Khth_Max)) enddo ; enddo else !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc(i,j) = max(CS%Khth_Min, Khth_loc(i,j)) + Khth_loc_v(i,J) = max(CS%Khth_Min, Khth_loc_v(i,J)) enddo ; enddo endif if (CS%max_Khth_CFL > 0.0) then !$OMP do do J=js-1,je ; do i=is,ie - KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_loc(i,j)) + KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_loc_v(i,J)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index f5ee25c743..5ed9e2a7a4 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -196,7 +196,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) !! False => compute G'(1) as in LMD94 if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, KPP_init: '// & 'Control structure has already been initialized') - allocate(CS) ! Read parameters call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // & @@ -207,6 +206,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) default=.false.) ! Forego remainder of initialization if not using this scheme if (.not. KPP_init) return + allocate(CS) call openParameterBlock(paramFile,'KPP') call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, & @@ -640,7 +640,12 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & buoy_scale = US%L_to_m**2*US%s_to_T**3 - !$OMP parallel do default(shared) firstprivate(nonLocalTrans) + !$OMP parallel do default(none) firstprivate(nonLocalTrans) & + !$OMP private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & + !$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, & + !$OMP sigmaRatio) & + !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, & + !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves) ! loop over horizontal points on processor do j = G%jsc, G%jec do i = G%isc, G%iec @@ -957,7 +962,16 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF buoy_scale = US%L_to_m**2*US%s_to_T**3 ! loop over horizontal points on processor - !$OMP parallel do default(shared) + !GOMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & + !GOMP surfBuoyFlux, U_H, V_H, u, v, Coriolis, pRef, SLdepth_0d, & + !GOMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, & + !GOMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, & + !GOMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, & + !GOMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D, & + !GOMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, & + !GOMP BulkRi_1d, zBottomMinusOffset) & + !GOMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & + !GOMP Temp, Salt, waves, EOS, GoRho) do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1463,7 +1477,7 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & dtracer(:,:,:) = 0.0 - !$OMP parallel do default(shared) + !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1522,7 +1536,7 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, dtracer(:,:,:) = 0.0 - !$OMP parallel do default(shared) + !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index 57400e31bf..6abd126ea2 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -145,8 +145,11 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) CS%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,Time, & 'Double-diffusion density ratio', 'nondim') - if (CS%id_R_rho > 0) & - allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)); CS%R_rho(:,:,:) = 0.0 + + if (CS%id_R_rho > 0) then + allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)) + CS%R_rho(:,:,:) = 0.0 + endif call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & kappa_ddiff_s=CS%kappa_ddiff_s, & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index f65a0e8eae..38cecf0425 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1332,6 +1332,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, nkmb = GV%nk_rho_varies h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 + ea_s(:,:,:) = 0.0; eb_s(:,:,:) = 0.0; ea_t(:,:,:) = 0.0; eb_t(:,:,:) = 0.0 showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic_ALE(), MOM_diabatic_driver.F90") @@ -2860,12 +2861,14 @@ end subroutine layered_diabatic !> Returns pointers or values of members within the diabatic_CS type. For extensibility, !! each returned argument is an optional argument -subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, & - evap_CFL_limit, minimum_forcing_depth, diabatic_aux_CSp) - type(diabatic_CS), intent(in ) :: CS !< module control structure +subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, & + minimum_forcing_depth, KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp) + type(diabatic_CS), intent(in ) :: CS !< module control structure ! All output arguments are optional type(opacity_CS), optional, pointer :: opacity_CSp !< A pointer to be set to the opacity control structure type(optics_type), optional, pointer :: optics_CSp !< A pointer to be set to the optics control structure + type(KPP_CS), optional, pointer :: KPP_CSp !< A pointer to be set to the KPP CS + type(energetic_PBL_CS), optional, pointer :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS real, optional, intent( out) :: evap_CFL_limit ! CS%opacity_CSp - if (present(optics_CSp)) optics_CSp => CS%optics - if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp + if (present(opacity_CSp)) opacity_CSp => CS%opacity_CSp + if (present(optics_CSp)) optics_CSp => CS%optics + if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp + if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL_CSp ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 0aaba9d3cf..921769091b 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -365,10 +365,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) if (.not.use_BBL_EOS) Rml_vel(:,:) = 0.0 - !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,is,ie,js,je,nz,nkmb, & + !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,nz,nkmb, & !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, & !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, & - !$OMP OBC,maxitt,Vol_quit,D_u,D_v,mask_u,mask_v) + !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v) & + !$OMP firstprivate(Vol_quit) do j=Jsq,Jeq ; do m=1,2 if (m==1) then diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 new file mode 100644 index 0000000000..82e0d6a559 --- /dev/null +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -0,0 +1,1126 @@ +!> Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by +!! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean. + +module MOM_lateral_boundary_diffusion + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_domains, only : pass_var +use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : post_data, register_diag_field +use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_grid, only : ocean_grid_type +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d +use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme +use MOM_tracer_registry, only : tracer_registry_type, tracer_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member + +implicit none ; private + +public near_boundary_unit_tests, lateral_boundary_diffusion, lateral_boundary_diffusion_init +public boundary_k_range + +! Private parameters to avoid doing string comparisons for bottom or top boundary layer +integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary +integer, public, parameter :: BOTTOM = 1 !< Set a value that corresponds to the bottom boundary +#include + +!> Sets parameters for lateral boundary mixing module. +type, public :: lateral_boundary_diffusion_CS ; private + integer :: method !< Determine which of the three methods calculate + !! and apply near boundary layer fluxes + !! 1. Bulk-layer approach + !! 2. Along layer + integer :: deg !< Degree of polynomial reconstruction + integer :: surface_boundary_scheme !< Which boundary layer scheme to use + !! 1. ePBL; 2. KPP + type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. +end type lateral_boundary_diffusion_CS + +! This include declares and sets the variable "version". +#include "version_variable.h" +character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module + +contains + +!> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be +!! needed for lateral boundary diffusion. +logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS) + type(time_type), target, intent(in) :: Time !< Time structure + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure + type(diabatic_CS), pointer :: diabatic_CSp !< KPP control structure needed to get BLD + type(lateral_boundary_diffusion_CS), pointer :: CS !< Lateral boundary mixing control structure + + ! local variables + character(len=80) :: string ! Temporary strings + logical :: boundary_extrap + + if (ASSOCIATED(CS)) then + call MOM_error(FATAL, "lateral_boundary_diffusion_init called with associated control structure.") + return + endif + + ! Log this module and master switch for turning it on/off + call log_version(param_file, mdl, version, & + "This module implements lateral diffusion of tracers near boundaries") + call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, & + "If true, enables the lateral boundary tracer's diffusion module.", & + default=.false.) + + if (.not. lateral_boundary_diffusion_init) then + return + endif + + allocate(CS) + CS%diag => diag + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + + CS%surface_boundary_scheme = -1 + if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then + call MOM_error(FATAL,"Lateral boundary diffusion is true, but no valid boundary layer scheme was found") + endif + + ! Read all relevant parameters and write them to the model log. + call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & + "Determine how to apply boundary lateral diffusion of tracers: \n"//& + "1. Bulk layer approach \n"//& + "2. Along layer approach", default=1) + call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & + "Use boundary extrapolation in LBD code", & + default=.false.) + call get_param(param_file, mdl, "LBD_REMAPPING_SCHEME", string, & + "This sets the reconstruction scheme used "//& + "for vertical remapping for all variables. "//& + "It can be one of the following schemes: "//& + trim(remappingSchemesDoc), default=remappingDefaultScheme) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) + +end function lateral_boundary_diffusion_init + +!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. +!! Two different methods are available: +!! Method 1: lower order representation, calculate fluxes from bulk layer integrated quantities. +!! Method 2: more straight forward, diffusion is applied layer by layer using only information +!! from neighboring cells. +subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) + type(ocean_grid_type), intent(inout) :: G !< Grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [m2] + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [m2] + real, intent(in) :: dt !< Tracer time step * I_numitts + !! (I_numitts in tracer_hordiff) + type(tracer_registry_type), pointer :: Reg !< Tracer registry + type(lateral_boundary_diffusion_CS), intent(in) :: CS !< Control structure for this module + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] + real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial + real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_E !< Edge values from reconstructions + real, dimension(SZK_(G),CS%deg+1) :: ppoly_S !< Slopes from reconstruction (placeholder) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uFlx !< Zonal flux of tracer [conc m^3] + real, dimension(SZIB_(G),SZJ_(G)) :: uFLx_bulk !< Total calculated bulk-layer u-flux for the tracer + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vFlx !< Meridional flux of tracer [conc m^3] + real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer + real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport + real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport + real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency !< tendency array for diagn + real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d !< depth integrated content tendency for diagn + type(tracer_type), pointer :: Tracer => NULL() !< Pointer to the current tracer + integer :: remap_method !< Reconstruction method + integer :: i,j,k,m !< indices to loop over + real :: Idt !< inverse of the time step [s-1] + + Idt = 1./dt + hbl(:,:) = 0. + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + + call pass_var(hbl,G%Domain) + do m = 1,Reg%ntr + tracer => Reg%tr(m) + + ! for diagnostics + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then + tendency(:,:,:) = 0.0 + endif + + do j = G%jsc-1, G%jec+1 + ! Interpolate state to interface + do i = G%isc-1, G%iec+1 + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & + ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) + enddo + enddo + ! Diffusive fluxes in the i-direction + uFlx(:,:,:) = 0. + vFlx(:,:,:) = 0. + uFlx_bulk(:,:) = 0. + vFlx_bulk(:,:) = 0. + + ! Method #1 + if ( CS%method == 1 ) then + do j=G%jsc,G%jec + do i=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), & + ppoly0_coefs(I,j,:,:), ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), & + ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:)) + endif + enddo + enddo + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), & + ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & + ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:)) + endif + enddo + enddo + ! Post tracer bulk diags + if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) + if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) + + ! Method #2 + elseif (CS%method == 2) then + do j=G%jsc,G%jec + do i=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), ppoly0_coefs(I,j,:,:), & + ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx(I,j,:)) + endif + enddo + enddo + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), & + ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx(i,J,:)) + endif + enddo + enddo + endif + + ! Update the tracer fluxes + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%mask2dT(i,j)>0.) then + tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & + (G%IareaT(i,j)/( h(i,j,k) + GV%H_subroundoff)) + + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then + tendency(i,j,k) = (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) )) * G%IareaT(i,j) * Idt + endif + + endif + enddo ; enddo ; enddo + + ! Post the tracer diagnostics + if (tracer%id_lbd_dfx>0) call post_data(tracer%id_lbd_dfx, uFlx*Idt, CS%diag) + if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) + if (tracer%id_lbd_dfx_2d>0) then + uwork_2d(:,:) = 0. + do k=1,GV%ke; do j=G%jsc,G%jec; do I=G%isc-1,G%iec + uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) + enddo; enddo; enddo + call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) + endif + + if (tracer%id_lbd_dfy_2d>0) then + vwork_2d(:,:) = 0. + do k=1,GV%ke; do J=G%jsc-1,G%jec; do i=G%isc,G%iec + vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) + enddo; enddo; enddo + call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) + endif + + ! post tendency of tracer content + if (tracer%id_lbdxy_cont > 0) then + call post_data(tracer%id_lbdxy_cont, tendency(:,:,:), CS%diag) + endif + + ! post depth summed tendency for tracer content + if (tracer%id_lbdxy_cont_2d > 0) then + tendency_2d(:,:) = 0. + do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k = 1, GV%ke + tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) + enddo + enddo ; enddo + call post_data(tracer%id_lbdxy_cont_2d, tendency_2d(:,:), CS%diag) + endif + + ! post tendency of tracer concentration; this step must be + ! done after posting tracer content tendency, since we alter + ! the tendency array. + if (tracer%id_lbdxy_conc > 0) then + do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) + enddo ; enddo ; enddo + call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) + endif + + enddo + +end subroutine lateral_boundary_diffusion + +!< Calculate bulk layer value of a scalar quantity as the thickness weighted average +real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, & + zeta_bot) + integer :: boundary !< SURFACE or BOTTOM [nondim] + integer :: nk !< Number of layers [nondim] + integer :: deg !< Degree of polynomial [nondim] + real, dimension(nk) :: h !< Layer thicknesses [m] + real :: hBLT !< Depth of the boundary layer [m] + real, dimension(nk) :: phi !< Scalar quantity + real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial + integer :: method !< Remapping scheme to use + + integer :: k_top !< Index of the first layer within the boundary + real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer + !! (0 if none, 1. if all). For the surface, this is always 0. because + !! integration starts at the surface [nondim] + integer :: k_bot !< Index of the last layer within the boundary + real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer + !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. + !! because integration starts at the bottom [nondim] + ! Local variables + real :: htot !< Running sum of the thicknesses (top to bottom) + integer :: k !< k indice + + + htot = 0. + bulk_average = 0. + if (hblt == 0.) return + if (boundary == SURFACE) then + htot = (h(k_bot) * zeta_bot) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot + do k = k_bot-1,1,-1 + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + elseif (boundary == BOTTOM) then + htot = (h(k_top) * zeta_top) + ! (note 1-zeta_top because zeta_top is the fraction of the layer) + bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot + do k = k_top+1,nk + bulk_average = bulk_average + phi(k)*h(k) + htot = htot + h(k) + enddo + else + call MOM_error(FATAL, "bulk_average: a valid boundary type must be provided.") + endif + + bulk_average = bulk_average / hBLT + +end function bulk_average + +!> Calculate the harmonic mean of two quantities +!! See \ref section_harmonic_mean. +real function harmonic_mean(h1,h2) + real :: h1 !< Scalar quantity + real :: h2 !< Scalar quantity + if (h1 + h2 == 0.) then + harmonic_mean = 0. + else + harmonic_mean = 2.*(h1*h2)/(h1+h2) + endif +end function harmonic_mean + +!> Find the k-index range corresponding to the layers that are within the boundary-layer region +subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot) + integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [m] + real, intent(in ) :: hbl !< Thickness of the boundary layer [m] + !! If surface, with respect to zbl_ref = 0. + !! If bottom, with respect to zbl_ref = SUM(h) + integer, intent( out) :: k_top !< Index of the first layer within the boundary + real, intent( out) :: zeta_top !< Distance from the top of a layer to the intersection of the + !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim] + integer, intent( out) :: k_bot !< Index of the last layer within the boundary + real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth + !! (0 at top, 1 at bottom) [nondim] + ! Local variables + real :: htot + integer :: k + ! Surface boundary layer + if ( boundary == SURFACE ) then + k_top = 1 + zeta_top = 0. + htot = 0. + k_bot = 1 + zeta_bot = 0. + if (hbl == 0.) return + if (hbl >= SUM(h(:))) then + k_bot = nk + zeta_bot = 0. + return + endif + do k=1,nk + htot = htot + h(k) + if ( htot >= hbl) then + k_bot = k + zeta_bot = 1 - (htot - hbl)/h(k) + return + endif + enddo + ! Bottom boundary layer + elseif ( boundary == BOTTOM ) then + k_top = nk + zeta_top = 1. + k_bot = nk + zeta_bot = 1. + htot = 0. + if (hbl == 0.) return + if (hbl >= SUM(h(:))) then + k_top = 1 + zeta_top = 0. + return + endif + do k=nk,1,-1 + htot = htot + h(k) + if (htot >= hbl) then + k_top = k + zeta_top = 1 - (htot - hbl)/h(k) + return + endif + enddo + else + call MOM_error(FATAL,"Houston, we've had a problem in boundary_k_range") + endif + +end subroutine boundary_k_range + + +!> Calculate the lateral boundary diffusive fluxes using the layer by layer method. +!! See \ref section_method2 +subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, & + ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, intent(in ) :: hbl_L !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: hbl_R !< Thickness of the boundary boundary + !! layer (right) [m] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] + integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point [m^3 conc] + + ! Local variables + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! This is just to remind developers that khtr_avg should be + !! computed once khtr is 3D. + real :: heff !< Harmonic mean of layer thicknesses [m] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) + !! [conc m^-3 ] + real :: htot !< Total column thickness [m] + integer :: k, k_bot_min, k_top_max !< k-indices, min and max for top and bottom, respectively + integer :: k_top_L, k_bot_L !< k-indices left + integer :: k_top_R, k_bot_R !< k-indices right + real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary + !! layer depth [nondim] + real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary + !!layer depth [nondim] + real :: h_work_L, h_work_R !< dummy variables + real :: hbl_min !< minimum BLD (left and right) [m] + + F_layer(:) = 0.0 + if (hbl_L == 0. .or. hbl_R == 0.) then + return + endif + + ! Calculate vertical indices containing the boundary layer + call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + + if (boundary == SURFACE) then + k_bot_min = MIN(k_bot_L, k_bot_R) + ! make sure left and right k indices span same range + if (k_bot_min .ne. k_bot_L) then + k_bot_L = k_bot_min + zeta_bot_L = 1.0 + endif + if (k_bot_min .ne. k_bot_R) then + k_bot_R= k_bot_min + zeta_bot_R = 1.0 + endif + + h_work_L = (h_L(k_bot_L) * zeta_bot_L) + h_work_R = (h_R(k_bot_R) * zeta_bot_R) + + phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_bot_L, 0., zeta_bot_L) + phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_bot_R, 0., zeta_bot_R) + heff = harmonic_mean(h_work_L, h_work_R) + ! tracer flux where the minimum BLD intersets layer + ! GMM, khtr_avg should be computed once khtr is 3D + F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) + + do k = k_bot_min-1,1,-1 + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + enddo + endif + + if (boundary == BOTTOM) then + k_top_max = MAX(k_top_L, k_top_R) + ! make sure left and right k indices span same range + if (k_top_max .ne. k_top_L) then + k_top_L = k_top_max + zeta_top_L = 1.0 + endif + if (k_top_max .ne. k_top_R) then + k_top_R= k_top_max + zeta_top_R = 1.0 + endif + + h_work_L = (h_L(k_top_L) * zeta_top_L) + h_work_R = (h_R(k_top_R) * zeta_top_R) + + phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, 1.0-zeta_top_L, 1.0) + phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, 1.0-zeta_top_R, 1.0) + heff = harmonic_mean(h_work_L, h_work_R) + + ! tracer flux where the minimum BLD intersets layer + F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) + + do k = k_top_max+1,nk + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + enddo + endif + +end subroutine fluxes_layer_method + +!> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' +!! See \ref section_method1 +subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, & + ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit) + + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: nk !< Number of layers [nondim] + integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, intent(in ) :: hbl_L !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: hbl_R !< Thickness of the boundary boundary + !! layer (left) [m] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] + integer, intent(in ) :: method !< Method of polynomial integration [nondim] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] + real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux [m^3 conc] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^3 conc] + real, optional, dimension(nk), intent( out) :: F_limit !< The amount of flux not applied due to limiter + !! F_layer(k) - F_max [m^3 conc] + ! Local variables + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! This is just to remind developers that khtr_avg should be + !! computed once khtr is 3D. + real :: heff !< Harmonic mean of layer thicknesses [m] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) + !! [conc m^-3 ] + real :: htot ! Total column thickness [m] + integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively + integer :: k_top_L, k_bot_L !< k-indices left + integer :: k_top_R, k_bot_R !< k-indices right + real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the + !! boundary layer [nondim] + real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the + !! boundary layer [nondim] + real :: h_work_L, h_work_R !< dummy variables + real :: F_max !< The maximum amount of flux that can leave a + !! cell [m^3 conc] + logical :: limited !< True if the flux limiter was applied + real :: hfrac, F_bulk_remain + + if (hbl_L == 0. .or. hbl_R == 0.) then + F_bulk = 0. + F_layer(:) = 0. + return + endif + + ! Calculate vertical indices containing the boundary layer + call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + + ! Calculate bulk averages of various quantities + phi_L_avg = bulk_average(boundary, nk, deg, h_L, hbl_L, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, & + zeta_top_L, k_bot_L, zeta_bot_L) + phi_R_avg = bulk_average(boundary, nk, deg, h_R, hbl_R, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, & + zeta_top_R, k_bot_R, zeta_bot_R) + + ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities + ! GMM, khtr_avg should be computed once khtr is 3D + heff = harmonic_mean(hbl_L, hbl_R) + F_bulk = -(khtr_u * heff) * (phi_R_avg - phi_L_avg) + F_bulk_remain = F_bulk + ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated + ! above, but is used as a way to decompose the fluxes onto the individual layers + h_means(:) = 0. + + if (boundary == SURFACE) then + k_min = MIN(k_bot_L, k_bot_R) + + ! left hand side + if (k_bot_L == k_min) then + h_work_L = h_L(k_min) * zeta_bot_L + else + h_work_L = h_L(k_min) + endif + + ! right hand side + if (k_bot_R == k_min) then + h_work_R = h_R(k_min) * zeta_bot_R + else + h_work_R = h_R(k_min) + endif + + h_means(k_min) = harmonic_mean(h_work_L,h_work_R) + + do k=1,k_min-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + + elseif (boundary == BOTTOM) then + k_max = MAX(k_top_L, k_top_R) + ! left hand side + if (k_top_L == k_max) then + h_work_L = h_L(k_max) * zeta_top_L + else + h_work_L = h_L(k_max) + endif + + ! right hand side + if (k_top_R == k_max) then + h_work_R = h_R(k_max) * zeta_top_R + else + h_work_R = h_R(k_max) + endif + + h_means(k_max) = harmonic_mean(h_work_L,h_work_R) + + do k=nk,k_max+1,-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + endif + + if ( SUM(h_means) == 0. ) then + return + ! Decompose the bulk flux onto the individual layers + else + ! Initialize remaining thickness + inv_heff = 1./SUM(h_means) + do k=1,nk + if (h_means(k) > 0.) then + hfrac = h_means(k)*inv_heff + F_layer(k) = F_bulk * hfrac + ! limit the flux to 0.2 of the tracer *gradient* + ! Why 0.2? + ! t=0 t=inf + ! 0 .2 + ! 0 1 0 .2.2.2 + ! 0 .2 + ! + F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) + + ! check if bulk flux (or F_layer) and F_max have same direction + if ( SIGN(1.,F_bulk) == SIGN(1., F_max)) then + ! Distribute bulk flux onto layers + if ( ((boundary == SURFACE) .and. (k == k_min)) .or. ((boundary == BOTTOM) .and. (k == nk)) ) then + F_layer(k) = F_bulk_remain ! GMM, are not using F_bulk_remain for now. Should we keep it? + endif + F_bulk_remain = F_bulk_remain - F_layer(k) + + ! Apply flux limiter calculated above + if (F_max >= 0.) then + limited = F_layer(k) > F_max + F_layer(k) = MIN(F_layer(k),F_max) + else + limited = F_layer(k) < F_max + F_layer(k) = MAX(F_layer(k),F_max) + endif + + ! GMM, again we are not using F_limit. Should we delete it? + if (PRESENT(F_limit)) then + if (limited) then + F_limit(k) = F_layer(k) - F_max + else + F_limit(k) = 0. + endif + endif + else + ! do not apply a flux on this layer + F_bulk_remain = F_bulk_remain - F_layer(k) + F_layer(k) = 0. + endif + else + F_layer(k) = 0. + endif + enddo + endif + +end subroutine fluxes_bulk_method + +!> Unit tests for near-boundary horizontal mixing +logical function near_boundary_unit_tests( verbose ) + logical, intent(in) :: verbose !< If true, output additional information for debugging unit tests + + ! Local variables + integer, parameter :: nk = 2 ! Number of layers + integer, parameter :: deg = 1 ! Degree of reconstruction (linear here) + integer, parameter :: method = 1 ! Method used for integrating polynomials + real, dimension(nk) :: phi_L, phi_R ! Tracer values (left and right column) [ nondim m^-3 ] + real, dimension(nk) :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) + real, dimension(nk,deg+1) :: phi_pp_L, phi_pp_R ! Coefficients for the linear pseudo-reconstructions + ! [ nondim m^-3 ] + + real, dimension(nk,2) :: ppoly0_E_L, ppoly0_E_R! Polynomial edge values (left and right) [concentration] + real, dimension(nk) :: h_L, h_R ! Layer thickness (left and right) [m] + real :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1] + real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] + real :: F_bulk ! Total diffusive flux across the U point [nondim s^-1] + real, dimension(nk) :: F_layer ! Diffusive flux within each layer at U-point [nondim s^-1] + real :: h_u, hblt_u ! Thickness at the u-point [m] + real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real :: heff ! Harmonic mean of layer thicknesses [m] + real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] + character(len=120) :: test_name ! Title of the unit test + integer :: k_top ! Index of cell containing top of boundary + real :: zeta_top ! Nondimension position + integer :: k_bot ! Index of cell containing bottom of boundary + real :: zeta_bot ! Nondimension position + real :: area_L,area_R ! Area of grid cell [m^2] + area_L = 1.; area_R = 1. ! Set to unity for all unit tests + + near_boundary_unit_tests = .false. + + ! Unit tests for boundary_k_range + test_name = 'Surface boundary spans the entire top cell' + h_L = (/5.,5./) + call boundary_k_range(SURFACE, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose) + + test_name = 'Surface boundary spans the entire column' + h_L = (/5.,5./) + call boundary_k_range(SURFACE, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., test_name, verbose) + + test_name = 'Bottom boundary spans the entire bottom cell' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0., 2, 1., test_name, verbose) + + test_name = 'Bottom boundary spans the entire column' + h_L = (/5.,5./) + call boundary_k_range(BOTTOM, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose) + + test_name = 'Surface boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose) + + test_name = 'Surface boundary intersects first layer' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose) + + test_name = 'Surface boundary is deeper than column thickness' + h_L = (/10.,10./) + call boundary_k_range(SURFACE, nk, h_L, 21.0, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., test_name, verbose) + + test_name = 'Bottom boundary intersects first layer' + h_L = (/10.,10./) + call boundary_k_range(BOTTOM, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 1., test_name, verbose) + + test_name = 'Bottom boundary intersects second layer' + h_L = (/10.,10./) + call boundary_k_range(BOTTOM, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot) + near_boundary_unit_tests = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 1., test_name, verbose) + + ! All cases in this section have hbl which are equal to the column thicknesses + test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R, & + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-5.0,-5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/0.,0./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 0.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 1. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. + ppoly0_E_R(2,1) = 0.; ppoly0_E_R(2,2) = 0. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/5.0,5.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (no gradient)' + hbl_L = 10; hbl_R = 10 + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,1./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Equal hbl and different layer thicknesses (gradient right to left)' + hbl_L = 16.; hbl_R = 16. + h_L = (/10.,6./) ; h_R = (/6.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-8.0,-8.0/) ) + + test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)' + hbl_L = 10.; hbl_R = 10. + h_L = (/5.,5./) ; h_R = (/5.,5./) + phi_L = (/1.,0./) ; phi_R = (/0.,1./) + phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + + test_name = 'Different hbl and different layer thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) + + ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction) + + test_name = 'hbl < column thickness, hbl same, constant concentration each column' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + khtr_u = 1. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + + test_name = 'hbl < column thickness, hbl same, linear profile right' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/0.5,2./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. + khtr_u = 1. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. + call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& + ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + + test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2' + hbl_L = 2; hbl_R = 2 + h_L = (/1.,2./) ; h_R = (/1.,2./) + phi_L = (/0.,0./) ; phi_R = (/0.5,2./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. + khtr_u = 2. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. + call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & + phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.,-2./) ) + + ! unit tests for layer by layer method + test_name = 'Different hbl and different column thicknesses (gradient from right to left)' + hbl_L = 12; hbl_R = 20 + h_L = (/6.,6./) ; h_R = (/10.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,1./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. + phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. + ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + khtr_u = 1. + call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & + phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,0.0/) ) + + test_name = 'Different hbl and different column thicknesses (linear profile right)' + + hbl_L = 15; hbl_R = 6 + h_L = (/10.,10./) ; h_R = (/12.,10./) + phi_L = (/0.,0./) ; phi_R = (/1.,3./) + phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. + phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. + phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 2. + phi_pp_R(2,1) = 2.; phi_pp_R(2,2) = 2. + ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. + ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. + ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 2. + ppoly0_E_R(2,1) = 2.; ppoly0_E_R(2,2) = 4. + khtr_u = 1. + call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & + phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + near_boundary_unit_tests = test_layer_fluxes( verbose, nk, test_name, F_layer, (/-3.75,0.0/) ) +end function near_boundary_unit_tests + +!> Returns true if output of near-boundary unit tests does not match correct computed values +!! and conditionally writes results to stream +logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + integer, intent(in) :: nk !< Number of layers + real, dimension(nk), intent(in) :: F_calc !< Fluxes of the unitless tracer from the algorithm [s^-1] + real, dimension(nk), intent(in) :: F_ans !< Fluxes of the unitless tracer calculated by hand [s^-1] + ! Local variables + integer :: k + integer, parameter :: stdunit = 6 + + test_layer_fluxes = .false. + do k=1,nk + if ( F_calc(k) /= F_ans(k) ) then + test_layer_fluxes = .true. + write(stdunit,*) "UNIT TEST FAILED: ", test_name + write(stdunit,10) k, F_calc(k), F_ans(k) + elseif (verbose) then + write(stdunit,10) k, F_calc(k), F_ans(k) + endif + enddo + +10 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16) +end function test_layer_fluxes + +!> Return true if output of unit tests for boundary_k_range does not match answers +logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,& + k_bot_ans, zeta_bot_ans, test_name, verbose) + integer :: k_top !< Index of cell containing top of boundary + real :: zeta_top !< Nondimension position + integer :: k_bot !< Index of cell containing bottom of boundary + real :: zeta_bot !< Nondimension position + integer :: k_top_ans !< Index of cell containing top of boundary + real :: zeta_top_ans !< Nondimension position + integer :: k_bot_ans !< Index of cell containing bottom of boundary + real :: zeta_bot_ans !< Nondimension position + character(len=80) :: test_name !< Name of the unit test + logical :: verbose !< If true always print output + + integer, parameter :: stdunit = 6 + + test_boundary_k_range = k_top .ne. k_top_ans + test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans) + test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans) + test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans) + + if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name + if (test_boundary_k_range .or. verbose) then + write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans + write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans + write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans + write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans + endif + + 20 format(A,"=",i3,X,A,"=",i3) + 30 format(A,"=",f20.16,X,A,"=",f20.16) + + +end function test_boundary_k_range + +!> \namespace mom_lateral_boundary_diffusion +!! +!! \section section_LBD The Lateral Boundary Diffusion (LBD) framework +!! +!! The LBD framework accounts for the effects of diabatic mesoscale fluxes +!! within surface and bottom boundary layers. Unlike the equivalent adiabatic +!! fluxes, which is applied along neutral density surfaces, LBD is purely +!! horizontal. +!! +!! The bottom boundary layer fluxes remain to be implemented, although most +!! of the steps needed to do so have already been added and tested. +!! +!! Boundary lateral diffusion can be applied using one of the three methods: +!! +!! * [Method #1: Bulk layer](@ref section_method1) (default); +!! * [Method #2: Along layer](@ref section_method2); +!! +!! A brief summary of these methods is provided below. +!! +!! \subsection section_method1 Bulk layer approach (Method #1) +!! +!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This +!! is a lower order representation (Kraus-Turner like approach) which assumes that +!! eddies are acting along well mixed layers (i.e., eddies do not know care about +!! vertical tracer gradients within the boundary layer). +!! +!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! For the TOP boundary layer, these are: +!! +!! k_top, k_bot, zeta_top, zeta_bot +!! +!! Step #2: compute bulk averages (thickness weighted) tracer averages (phi_L and phi_R), +!! then calculate the bulk diffusive flux (F_{bulk}): +!! +!! \f[ F_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L), \f] +!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the boundary layer depth +!! in the left and right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively). +!! +!! Step #3: decompose F_bulk onto individual layers: +!! +!! \f[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) , \f] +!! +!! where h_{frac} is +!! +!! \f[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}. \f] +!! +!! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer. +!! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R). +!! +!! Step #4: limit the tracer flux so that 1) only down-gradient fluxes are applied, +!! and 2) the flux cannot be larger than F_max, which is defined using the tracer +!! gradient: +!! +!! \f[ F_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))], \f] +!! where V is the cell volume. Why 0.2? +!! t=0 t=inf +!! 0 .2 +!! 0 1 0 .2.2.2 +!! 0 .2 +!! +!! \subsection section_method2 Along layer approach (Method #2) +!! +!! This is a more straight forward method where diffusion is applied layer by layer using +!! only information from neighboring cells. +!! +!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! For the TOP boundary layer, these are: +!! +!! k_top, k_bot, zeta_top, zeta_bot +!! +!! Step #2: calculate the diffusive flux at each layer: +!! +!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f] +!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness +!! in the left and right columns. Special care (layer reconstruction) must be taken at +!! k_min = min(k_botL, k_bot_R). This method does not require a limiter since KHTR +!! is already limted based on a diffusive CFL condition prior to the call of this +!! module. +!! +!! \subsection section_harmonic_mean Harmonic Mean +!! +!! The harmonic mean (HM) betwen h1 and h2 is defined as: +!! +!! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f] +!! +end module MOM_lateral_boundary_diffusion diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 0ffff7409d..f569c81bbc 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -5,6 +5,7 @@ module MOM_neutral_diffusion use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_domains, only : pass_var use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_diag_mediator, only : post_data, register_diag_field use MOM_EOS, only : EOS_type, EOS_manual_init, calculate_compress, calculate_density_derivs @@ -23,7 +24,10 @@ module MOM_neutral_diffusion use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation use regrid_edge_values, only : edge_values_implicit_h4 - +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member +use MOM_lateral_boundary_diffusion, only : boundary_k_range, SURFACE, BOTTOM implicit none ; private #include @@ -39,11 +43,13 @@ module MOM_neutral_diffusion integer :: deg = 2 !< Degree of polynomial used for reconstructions logical :: continuous_reconstruction = .true. !< True if using continuous PPM reconstruction at interfaces logical :: debug = .false. !< If true, write verbose debugging messages + logical :: hard_fail_heff !< Bring down the model if a problem with heff is detected integer :: max_iter !< Maximum number of iterations if refine_position is defined real :: drho_tol !< Convergence criterion representing difference from true neutrality real :: x_tol !< Convergence criterion for how small an update of the position can be real :: ref_pres !< Reference pressure, negative if using locally referenced neutral density - + logical :: interior_only !< If true, only applies neutral diffusion in the ocean interior. + !! That is, the algorithm will exclude the surface and bottom boundary layers. ! Positions of neutral surfaces in both the u, v directions real, allocatable, dimension(:,:,:) :: uPoL !< Non-dimensional position with left layer uKoL-1, u-point real, allocatable, dimension(:,:,:) :: uPoR !< Non-dimensional position with right layer uKoR-1, u-point @@ -88,6 +94,8 @@ module MOM_neutral_diffusion real :: C_p !< heat capacity of seawater (J kg-1 K-1) type(EOS_type), pointer :: EOS !< Equation of state parameters type(remapping_CS) :: remap_CS !< Remapping control structure used to create sublayers + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get MLD end type neutral_diffusion_CS ! This include declares and sets the variable "version". @@ -97,12 +105,13 @@ module MOM_neutral_diffusion contains !> Read parameters and allocate control structure for neutral_diffusion module. -logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) +logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< Time structure type(ocean_grid_type), intent(in) :: G !< Grid structure type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), target, intent(in) :: EOS !< Equation of state + type(diabatic_CS), pointer :: diabatic_CSp!< KPP control structure needed to get BLD type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure ! Local variables @@ -115,6 +124,7 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) return endif + ! Log this module and master switch for turning it on/off call log_version(param_file, mdl, version, & "This module implements neutral diffusion of tracers") @@ -143,6 +153,11 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) "the equation of state. If negative (default), local "//& "pressure is used.", & default = -1.) + call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", CS%interior_only, & + "If true, only applies neutral diffusion in the ocean interior."//& + "That is, the algorithm will exclude the surface and bottom"//& + "boundary layers.",default = .false.) + ! Initialize and configure remapping if (CS%continuous_reconstruction .eqv. .false.) then call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, & @@ -191,6 +206,17 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS) "Turns on verbose output for discontinuous neutral "//& "diffusion routines.", & default = .false.) + call get_param(param_file, mdl, "HARD_FAIL_HEFF", CS%hard_fail_heff, & + "Bring down the model if a problem with heff is detected",& + default = .true.) + endif + + if (CS%interior_only) then + call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) + call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then + call MOM_error(FATAL,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found") + endif endif ! call get_param(param_file, mdl, "KHTR", CS%KhTr, & @@ -234,9 +260,10 @@ end function neutral_diffusion_init !> Calculate remapping factors for u/v columns used to map adjoining columns to !! a shared coordinate space. -subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) +subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS) 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 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T !< Potential temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S !< Salinity [ppt] @@ -247,14 +274,36 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) ! Variables used for reconstructions real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] integer :: iMethod real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta - real, dimension(SZI_(G)) :: rho_tmp ! Routiine to calculate drho_dp, returns density which is not used + real, dimension(SZI_(G)) :: rho_tmp ! Routine to calculate drho_dp, returns density which is not used real :: h_neglect, h_neglect_edge + integer, dimension(SZI_(G), SZJ_(G)) :: k_top !< Index of the first layer within the boundary + real, dimension(SZI_(G), SZJ_(G)) :: zeta_top !< Distance from the top of a layer to the intersection of the + !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim] + integer, dimension(SZI_(G), SZJ_(G)) :: k_bot !< Index of the last layer within the boundary + real, dimension(SZI_(G), SZJ_(G)) :: zeta_bot !< Distance of the lower layer to the boundary layer depth real :: pa_to_H pa_to_H = 1. / GV%H_to_pa + k_top(:,:) = 1 ; k_bot(:,:) = 1 + zeta_top(:,:) = 0. ; zeta_bot(:,:) = 1. + + ! check if hbl needs to be extracted + if (CS%interior_only) then + hbl(:,:) = 0. + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + call pass_var(hbl,G%Domain) + ! get k-indices and zeta + do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 + call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) + enddo; enddo + ! TODO: add similar code for BOTTOM boundary layer + endif + !### Try replacing both of these with GV%H_subroundoff if (GV%Boussinesq) then h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 @@ -346,7 +395,14 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) if (.not. CS%continuous_reconstruction) then do j = G%jsc-1, G%jec+1 ; do i = G%isc-1, G%iec+1 call mark_unstable_cells( CS, G%ke, CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%P_i(i,j,:,:), CS%stable_cell(i,j,:) ) - enddo ; enddo + if (CS%interior_only) then + if (.not. CS%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1. + ! set values in the surface and bottom boundary layer to false. + do k = 1, k_bot(i,j)-1 + CS%stable_cell(i,j,k) = .false. + enddo + endif + enddo ; enddo endif CS%uhEff(:,:,:) = 0. @@ -367,14 +423,16 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) call find_neutral_surface_positions_continuous(G%ke, & CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i+1,j,:), CS%Tint(i+1,j,:), CS%Sint(i+1,j,:), CS%dRdT(i+1,j,:), CS%dRdS(i+1,j,:), & - CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:) ) + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & + k_bot(I,j), k_bot(I+1,j), 1.-zeta_bot(I,j), 1.-zeta_bot(I+1,j)) else call find_neutral_surface_positions_discontinuous(CS, G%ke, & CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), & CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), & CS%P_i(i+1,j,:,:), h(i+1,j,:), CS%T_i(i+1,j,:,:), CS%S_i(i+1,j,:,:), CS%ppoly_coeffs_T(i+1,j,:,:), & - CS%ppoly_coeffs_S(i+1,j,:,:), CS%stable_cell(i+1,j,:), & - CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:)) + CS%ppoly_coeffs_S(i+1,j,:,:), CS%stable_cell(i+1,j,:), & + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & + hard_fail_heff = CS%hard_fail_heff) endif endif enddo ; enddo @@ -383,17 +441,19 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS) do J = G%jsc-1, G%jec ; do i = G%isc, G%iec if (G%mask2dCv(i,J) > 0.) then if (CS%continuous_reconstruction) then - call find_neutral_surface_positions_continuous(G%ke, & + call find_neutral_surface_positions_continuous(G%ke, & CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), & - CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:) ) + CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), & + k_bot(i,J), k_bot(i,J+1), 1.-zeta_bot(i,J), 1.-zeta_bot(i,J+1)) else call find_neutral_surface_positions_discontinuous(CS, G%ke, & CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), & CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), & CS%P_i(i,j+1,:,:), h(i,j+1,:), CS%T_i(i,j+1,:,:), CS%S_i(i,j+1,:,:), CS%ppoly_coeffs_T(i,j+1,:,:), & CS%ppoly_coeffs_S(i,j+1,:,:), CS%stable_cell(i,j+1,:), & - CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:)) + CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:), & + hard_fail_heff = CS%hard_fail_heff) endif endif enddo ; enddo @@ -830,7 +890,7 @@ end function fvlsq_slope !> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, & - dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff) + dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff, bl_kl, bl_kr, bl_zl, bl_zr) integer, intent(in) :: nk !< Number of levels real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa] real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature [degC] @@ -849,6 +909,10 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa] + integer, optional, intent(in) :: bl_kl !< Layer index of the boundary layer (left) + integer, optional, intent(in) :: bl_kr !< Layer index of the boundary layer (right) + real, optional, intent(in) :: bl_zl !< Nondimensional position of the boundary layer (left) + real, optional, intent(in) :: bl_zr !< Nondimensional position of the boundary layer (right) ! Local variables integer :: ns ! Number of neutral surfaces @@ -863,13 +927,22 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS real :: dRho, dRhoTop, dRhoBot, hL, hR integer :: lastK_left, lastK_right real :: lastP_left, lastP_right + logical :: interior_limit ns = 2*nk+2 + ! Initialize variables for the search - kr = 1 ; lastK_right = 1 ; lastP_right = 0. - kl = 1 ; lastK_left = 1 ; lastP_left = 0. + kr = 1 ; + kl = 1 ; + lastP_right = 0. + lastP_left = 0. + lastK_right = 1 + lastK_left = 1 reached_bottom = .false. + ! Check to see if we should limit the diffusion to the interior + interior_limit = PRESENT(bl_kl) .and. PRESENT(bl_kr) .and. PRESENT(bl_zr) .and. PRESENT(bl_zl) + ! Loop over each neutral surface, working from top to bottom neutral_surfaces: do k_surface = 1, ns klm1 = max(kl-1, 1) @@ -988,10 +1061,23 @@ subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdS else stop 'Else what?' endif + if (interior_limit) then + if (KoL(k_surface)<=bl_kl) then + KoL(k_surface) = bl_kl + if (PoL(k_surface) 0) then - if (x_first) then +!$OMP parallel default(private) shared(nz,domore_k,x_first,Tr,hprev,uhr,uh_neglect, & +!$OMP OBC,domore_u,ntr,Idt,isv,iev,jsv,jev,stencil, & +!$OMP G,GV,CS,vhr,vh_neglect,domore_v,US) + if (x_first) then + + !$OMP do ordered + do k=1,nz ; if (domore_k(k) > 0) then ! First, advect zonally. call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & isv, iev, jsv-stencil, jev+stencil, k, G, GV, US, CS%usePPM, CS%useHuynh) + endif ; enddo + !$OMP do ordered + do k=1,nz ; if (domore_k(k) > 0) then ! Next, advect meridionally. call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + ! Update domore_k(k) for the next iteration domore_k(k) = 0 do j=jsv-stencil,jev+stencil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo do J=jsv-1,jev ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo - else + endif ; enddo + else + + !$OMP do ordered + do k=1,nz ; if (domore_k(k) > 0) then ! First, advect meridionally. call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & isv-stencil, iev+stencil, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + endif ; enddo + !$OMP do ordered + do k=1,nz ; if (domore_k(k) > 0) then ! Next, advect zonally. call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + ! Update domore_k(k) for the next iteration domore_k(k) = 0 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo do J=jsv-1,jev ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo + endif ; enddo - endif - + endif ! x_first - endif ; enddo ! End of k-loop +!$OMP end parallel ! If the advection just isn't finishing after max_iter, move on. if (itt >= max_iter) then @@ -334,7 +347,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & !! tracer change [H L2 ~> m3 or kg] real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhr !< accumulated volume/mass flux through !! the zonal face [H L2 ~> m3 or kg] - real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: uh_neglect !< A tiny zonal mass flux that can + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_neglect !< A tiny zonal mass flux that can !! be neglected [H L2 ~> m3 or kg] type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used logical, dimension(SZJ_(G),SZK_(G)), intent(inout) :: domore_u !< If true, there is more advection to be @@ -353,7 +366,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & real, dimension(SZI_(G),ntr) :: & slope_x ! The concentration slope per grid point [conc]. - real, dimension(SZIB_(G),ntr) :: & + real, dimension(SZIB_(G),SZJ_(G),ntr) :: & flux_x ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc]. real, dimension(SZI_(G),ntr) :: & T_tmp ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc]. @@ -374,13 +387,18 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! any of the passes [H ~> m or kg m-2]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. - logical :: do_i(SZIB_(G)) ! If true, work on given points. + logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points. logical :: do_any_i integer :: i, j, m, n, i_up, stencil real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6 real :: fac1,u_L_in,u_L_out ! terms used for time-stepping OBC reservoirs type(OBC_segment_type), pointer :: segment=>NULL() logical :: usePLMslope + logical, dimension(SZJ_(G),SZK_(G)) :: domore_u_initial + + ! keep a local copy of the initial values of domore_u, which is to be used when computing ad2d_x + ! diagnostic at the end of this subroutine. + domore_u_initial = domore_u usePLMslope = .not. (usePPM .and. useHuynh) ! stencil for calculating slope values @@ -537,10 +555,10 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & a6 = 6.*Tc - 3. * (aR + aL) ! Curvature if (uhh(I) >= 0.0) then - flux_x(I,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( & + flux_x(I,j,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( & ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) else - flux_x(I,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( & + flux_x(I,j,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( & ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) endif enddo ; enddo @@ -550,28 +568,28 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! Indirect implementation of PLM !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m) !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) - !flux_x(I,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) + !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) - !flux_x(I,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) ) + !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) ) ! Alternative implementation of PLM Tc = T_tmp(i,m) - flux_x(I,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) + flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) ! Original implementation of PLM - !flux_x(I,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I)) + !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I)) else ! Indirect implementation of PLM !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) - !flux_x(I,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) + !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) - !flux_x(I,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) ) + !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) ) ! Alternative implementation of PLM Tc = T_tmp(i+1,m) - flux_x(I,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) + flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) ! Original implementation of PLM - !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) + !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) endif !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j))) enddo ; enddo @@ -593,8 +611,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! should the reservoir evolve for this case Kate ?? - Nope do m=1,ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then - flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k) - else ; flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif + flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k) + else ; flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif enddo endif endif @@ -616,8 +634,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & uhh(I) = uhr(I,j,k) do m=1,ntr if (associated(segment%tr_Reg%Tr(m)%tres)) then - flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k) - else; flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif + flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k) + else; flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif enddo endif endif @@ -633,16 +651,16 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo do i=is,ie if ((uhh(I) /= 0.0) .or. (uhh(I-1) /= 0.0)) then - do_i(i) = .true. + do_i(i,j) = .true. hlst(i) = hprev(i,j,k) hprev(i,j,k) = hprev(i,j,k) - (uhh(I) - uhh(I-1)) - if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false. + if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false. elseif (hprev(i,j,k) < h_neglect*G%areaT(i,j)) then hlst(i) = hlst(i) + (h_neglect*G%areaT(i,j) - hprev(i,j,k)) Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) else ; Ihnew(i) = 1.0 / hprev(i,j,k) ; endif else - do_i(i) = .false. + do_i(i,j) = .false. endif enddo @@ -651,34 +669,46 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! update tracer do i=is,ie - if (do_i(i)) then + if (do_i(i,j)) then if (Ihnew(i) > 0.0) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & - (flux_x(I,m) - flux_x(I-1,m))) * Ihnew(i) + (flux_x(I,j,m) - flux_x(I-1,j,m))) * Ihnew(i) endif endif enddo ! diagnostics - if (associated(Tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i)) then - Tr(m)%ad_x(I,j,k) = Tr(m)%ad_x(I,j,k) + flux_x(I,m)*Idt - endif ; enddo ; endif - if (associated(Tr(m)%ad2d_x)) then ; do i=is,ie ; if (do_i(i)) then - Tr(m)%ad2d_x(I,j) = Tr(m)%ad2d_x(I,j) + flux_x(I,m)*Idt + if (associated(Tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i,j)) then + Tr(m)%ad_x(I,j,k) = Tr(m)%ad_x(I,j,k) + flux_x(I,j,m)*Idt endif ; enddo ; endif ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic). ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. if (associated(Tr(m)%advection_xy)) then - do i=is,ie ; if (do_i(i)) then - Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - (flux_x(I,m) - flux_x(I-1,m)) * & + do i=is,ie ; if (do_i(i,j)) then + Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - (flux_x(I,j,m) - flux_x(I-1,j,m)) * & Idt * G%IareaT(i,j) endif ; enddo endif enddo + endif + + + enddo ! End of j-loop. + + ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active. + + !$OMP ordered + do j=js,je ; if (domore_u_initial(j,k)) then + do m=1,ntr + if (associated(Tr(m)%ad2d_x)) then ; do i=is,ie ; if (do_i(i,j)) then + Tr(m)%ad2d_x(I,j) = Tr(m)%ad2d_x(I,j) + flux_x(I,j,m)*Idt + endif ; enddo ; endif + enddo endif ; enddo ! End of j-loop. + !$OMP end ordered end subroutine advect_x @@ -733,7 +763,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. - logical :: do_i(SZIB_(G)) ! If true, work on given points. + logical :: do_i(SZIB_(G), SZJ_(G)) ! If true, work on given points. logical :: do_any_i integer :: i, j, j2, m, n, j_up, stencil real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6 @@ -1010,36 +1040,33 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & do j=js,je ; if (do_j_tr(j)) then do i=is,ie if ((vhh(i,J) /= 0.0) .or. (vhh(i,J-1) /= 0.0)) then - do_i(i) = .true. + do_i(i,j) = .true. hlst(i) = hprev(i,j,k) hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,J) - vhh(i,J-1)), 0.0) - if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false. + if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false. elseif (hprev(i,j,k) < h_neglect*G%areaT(i,j)) then hlst(i) = hlst(i) + (h_neglect*G%areaT(i,j) - hprev(i,j,k)) Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) else ; Ihnew(i) = 1.0 / hprev(i,j,k) ; endif - else ; do_i(i) = .false. ; endif + else ; do_i(i,j) = .false. ; endif enddo ! update tracer and save some diagnostics do m=1,ntr - do i=is,ie ; if (do_i(i)) then + do i=is,ie ; if (do_i(i,j)) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & (flux_y(i,m,J) - flux_y(i,m,J-1))) * Ihnew(i) endif ; enddo ! diagnostics - if (associated(Tr(m)%ad_y)) then ; do i=is,ie ; if (do_i(i)) then + if (associated(Tr(m)%ad_y)) then ; do i=is,ie ; if (do_i(i,j)) then Tr(m)%ad_y(i,J,k) = Tr(m)%ad_y(i,J,k) + flux_y(i,m,J)*Idt endif ; enddo ; endif - if (associated(Tr(m)%ad2d_y)) then ; do i=is,ie ; if (do_i(i)) then - Tr(m)%ad2d_y(i,J) = Tr(m)%ad2d_y(i,J) + flux_y(i,m,J)*Idt - endif ; enddo ; endif ! diagnose convergence of flux_y and add to convergence of flux_x. ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. if (associated(Tr(m)%advection_xy)) then - do i=is,ie ; if (do_i(i)) then + do i=is,ie ; if (do_i(i,j)) then Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - (flux_y(i,m,J) - flux_y(i,m,J-1))* Idt * & G%IareaT(i,j) endif ; enddo @@ -1048,6 +1075,18 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & enddo endif ; enddo ! End of j-loop. + ! compute ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active. + + !$OMP ordered + do j=js,je ; if (do_j_tr(j)) then + do m=1,ntr + if (associated(Tr(m)%ad2d_y)) then ; do i=is,ie ; if (do_i(i,j)) then + Tr(m)%ad2d_y(i,J) = Tr(m)%ad2d_y(i,J) + flux_y(i,m,J)*Idt + endif ; enddo ; endif + enddo + endif ; enddo ! End of j-loop. + !$OMP end ordered + end subroutine advect_y !> Initialize lateral tracer advection module diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 2d42483c49..ff98431736 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -11,6 +11,7 @@ module MOM_tracer_hor_diff use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : pass_vector use MOM_debugging, only : hchksum, uvchksum +use MOM_diabatic_driver, only : diabatic_CS use MOM_EOS, only : calculate_density, EOS_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -22,6 +23,8 @@ module MOM_tracer_hor_diff use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end use MOM_neutral_diffusion, only : neutral_diffusion_CS use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion +use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion_CS, lateral_boundary_diffusion_init +use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs @@ -56,9 +59,13 @@ module MOM_tracer_hor_diff !! the CFL limit is not violated. logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. + logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within + !! tracer_hor_diff. logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been !! exceeded type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. + type(lateral_boundary_diffusion_CS), pointer :: lateral_boundary_diffusion_CSp => NULL() !< Control structure for + !! lateral boundary mixing. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -380,6 +387,31 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online endif enddo + if (CS%use_lateral_boundary_diffusion) then + + if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary mixing (tracer_hordiff)") + + call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) + + do J=js-1,je ; do i=is,ie + Coef_y(i,J) = I_numitts * khdt_y(i,J) + enddo ; enddo + do j=js,je + do I=is-1,ie + Coef_x(I,j) = I_numitts * khdt_x(I,j) + enddo + enddo + + do itt=1,num_itts + if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary diffusion (tracer_hordiff)",itt) + if (itt>1) then ! Update halos for subsequent iterations + call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) + endif + call lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, & + CS%lateral_boundary_diffusion_CSp) + enddo ! itt + endif + if (CS%use_neutral_diffusion) then if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)") @@ -389,7 +421,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs() ! would be inside the itt-loop. -AJA - call neutral_diffusion_calc_coeffs(G, GV, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) do J=js-1,je ; do i=is,ie Coef_y(i,J) = I_numitts * khdt_y(i,J) enddo ; enddo @@ -404,7 +436,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online if (itt>1) then ! Update halos for subsequent iterations call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) if (CS%recalc_neutral_surf) then - call neutral_diffusion_calc_coeffs(G, GV, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) endif endif call neutral_diffusion(G, GV, h, Coef_x, Coef_y, I_numitts*dt, Reg, US, CS%neutral_diffusion_CSp) @@ -1386,12 +1418,13 @@ end subroutine tracer_epipycnal_ML_diff !> Initialize lateral tracer diffusion module -subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, CS) +subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< current model time type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control type(EOS_type), target, intent(in) :: EOS !< Equation of state CS + type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS type(param_file_type), intent(in) :: param_file !< parameter file type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure @@ -1462,9 +1495,14 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, CS) units="nondim", default=1.0) endif - CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, CS%neutral_diffusion_CSp) + CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, param_file, diag, EOS, diabatic_CSp, & + CS%neutral_diffusion_CSp ) if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") + CS%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, & + CS%lateral_boundary_diffusion_CSp) + if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 01d15fb887..9229074099 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -57,6 +57,18 @@ module MOM_tracer_registry !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_y => NULL() !< diagnostic array for y-diffusive tracer flux !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: lbd_dfx => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: lbd_dfy => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_dfx_2d => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_dfy_2d => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_bulk_df_x => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: lbd_bulk_df_y => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_x => NULL() !< diagnostic vertical sum x-diffusive flux !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_y => NULL() !< diagnostic vertical sum y-diffusive flux @@ -110,9 +122,12 @@ module MOM_tracer_registry !>@{ Diagnostic IDs integer :: id_tr = -1 integer :: id_adx = -1, id_ady = -1, id_dfx = -1, id_dfy = -1 + integer :: id_lbd_bulk_dfx = -1, id_lbd_bulk_dfy = -1, id_lbd_dfx = -1, id_lbd_dfy = -1 + integer :: id_lbd_dfx_2d = -1 , id_lbd_dfy_2d = -1 integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1 integer :: id_adv_xy = -1, id_adv_xy_2d = -1 integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1 + integer :: id_lbdxy_cont = -1, id_lbdxy_cont_2d = -1, id_lbdxy_conc = -1 integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1 integer :: id_tendency = -1, id_trxh_tendency = -1, id_trxh_tendency_2d = -1 integer :: id_tr_vardec = -1 @@ -410,6 +425,14 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) diag%axesCvL, Time, trim(flux_longname)//" diffusive zonal flux" , & trim(flux_units), v_extensive = .true., x_cell_method = 'sum', & conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & + diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux from the lateral boundary diffusion "& + "scheme", trim(flux_units), v_extensive = .true., y_cell_method = 'sum', & + conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & + diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional flux from the lateral boundary diffusion"& + " scheme", trim(flux_units), v_extensive = .true., x_cell_method = 'sum', & + conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & @@ -425,11 +448,21 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) diag%axesCvL, Time, "Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') + Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & + diag%axesCuL, Time, "Lateral Boundary Diffusive Zonal Flux of "//trim(flux_longname), & + flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + y_cell_method='sum') + Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & + diag%axesCvL, Time, "Lateral Boundary Diffusive Meridional Flux of "//trim(flux_longname), & + flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + x_cell_method='sum') endif if (Tr%id_adx > 0) call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_dfy > 0) call safe_alloc_ptr(Tr%df_y,isd,ied,JsdB,JedB,nz) + if (Tr%id_lbd_dfx > 0) call safe_alloc_ptr(Tr%lbd_dfx,IsdB,IedB,jsd,jed,nz) + if (Tr%id_lbd_dfy > 0) call safe_alloc_ptr(Tr%lbd_dfy,isd,ied,JsdB,JedB,nz) Tr%id_adx_2d = register_diag_field("ocean_model", trim(shortnm)//"_adx_2d", & diag%axesCu1, Time, & @@ -449,11 +482,33 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') + Tr%id_lbd_bulk_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffx", & + diag%axesCu1, Time, & + "Total Bulk Diffusive Zonal Flux of "//trim(flux_longname), & + flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + y_cell_method='sum') + Tr%id_lbd_bulk_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffy", & + diag%axesCv1, Time, & + "Total Bulk Diffusive Meridional Flux of "//trim(flux_longname), & + flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + x_cell_method='sum') + Tr%id_lbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx_2d", & + diag%axesCu1, Time, "Vertically-integrated zonal diffusive flux from the lateral boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + y_cell_method = 'sum') + Tr%id_lbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy_2d", & + diag%axesCv1, Time, "Vertically-integrated meridional diffusive flux from the lateral boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + x_cell_method = 'sum') if (Tr%id_adx_2d > 0) call safe_alloc_ptr(Tr%ad2d_x,IsdB,IedB,jsd,jed) if (Tr%id_ady_2d > 0) call safe_alloc_ptr(Tr%ad2d_y,isd,ied,JsdB,JedB) if (Tr%id_dfx_2d > 0) call safe_alloc_ptr(Tr%df2d_x,IsdB,IedB,jsd,jed) if (Tr%id_dfy_2d > 0) call safe_alloc_ptr(Tr%df2d_y,isd,ied,JsdB,JedB) + if (Tr%id_lbd_bulk_dfx > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_x,IsdB,IedB,jsd,jed) + if (Tr%id_lbd_bulk_dfy > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_y,isd,ied,JsdB,JedB) + if (Tr%id_lbd_dfx_2d > 0) call safe_alloc_ptr(Tr%lbd_dfx_2d,IsdB,IedB,jsd,jed) + if (Tr%id_lbd_dfy_2d > 0) call safe_alloc_ptr(Tr%lbd_dfy_2d,isd,ied,JsdB,JedB) Tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", & diag%axesTL, Time, & @@ -478,37 +533,60 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) enddo ; enddo ; enddo endif - ! Lateral diffusion convergence tendencies + ! Neutral/Lateral diffusion convergence tendencies if (Tr%diag_form == 1) then Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer content tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer concentration "//& + diag%axesT1, Time, "Depth integrated neutral diffusion tracer content "//& + "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method='sum', y_cell_method= 'sum') + + Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + + Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated lateral diffusion tracer content "//& "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method= 'sum') else cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//& ' expressed as '//trim(lowercase(flux_longname))//& - ' content due to parameterized mesoscale diffusion' + ' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & - diag%axesTL, Time, "Lateral or neutral diffusion tracer concentration tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, cmor_field_name = trim(Tr%cmor_tendprefix)//'pmdiff', & cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), & x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& - trim(lowercase(flux_longname))//' content due to parameterized mesoscale diffusion' + trim(lowercase(flux_longname))//' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral or neutral diffusion tracer "//& - "concentration tendency for "//trim(shortnm), conv_units, & + diag%axesT1, Time, "Depth integrated neutral diffusion tracer "//& + "content tendency for "//trim(shortnm), conv_units, & conversion=Tr%conv_scale*US%s_to_T, cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff_2d', & cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & x_cell_method='sum', y_cell_method='sum') + + Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.) + + Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated lateral diffusion tracer "//& + "content tendency for "//trim(shortnm), conv_units, & + conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum') endif Tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', & - diag%axesTL, Time, "Lateral (neutral) tracer concentration tendency for "//trim(shortnm), & + diag%axesTL, Time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), & + trim(units)//' s-1', conversion=US%s_to_T) + + Tr%id_lbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_conc_tendency', & + diag%axesTL, Time, "Lateral diffusion tracer concentration tendency for "//trim(shortnm), & trim(units)//' s-1', conversion=US%s_to_T) var_lname = "Net time tendency for "//lowercase(flux_longname)