diff --git a/.testing/Makefile b/.testing/Makefile index adc1a20a1e..18687a1616 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -282,7 +282,7 @@ $(BUILD)/%/MOM6: $(BUILD)/%/Makefile FORCE # Target codebase should use its own build system $(BUILD)/target/MOM6: $(BUILD)/target FORCE | $(TARGET_CODEBASE) - $(MAKE) -C $(TARGET_CODEBASE)/.testing build/symmetric/MOM6 + $(MAKE) -C $(TARGET_CODEBASE)/.testing BUILD=build build/symmetric/MOM6 $(BUILD)/target: | $(TARGET_CODEBASE) ln -s $(abspath $(TARGET_CODEBASE))/.testing/build/symmetric $@ diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 3e3abba674..fab6fe1f55 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -2371,6 +2371,11 @@ subroutine SetScalarField(field, rc) ungriddedLBound=(/1/), ungriddedUBound=(/scalar_field_count/), gridToFieldMap=(/2/), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! initialize fldptr to zero + call ESMF_FieldGet(field, farrayPtr=fldptr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + fldptr2d(:,:) = 0.0 + end subroutine SetScalarField end subroutine MOM_RealizeFields diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 47971029cc..b696213e1b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2706,7 +2706,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%tv%T => CS%T ; CS%tv%S => CS%S if (CS%tv%T_is_conT) then vd_T = var_desc(name="contemp", units="Celsius", longname="Conservative Temperature", & - cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", & + cmor_field_name="bigthetao", cmor_longname="Sea Water Conservative Temperature", & conversion=US%Q_to_J_kg*CS%tv%C_p) else vd_T = var_desc(name="temp", units="degC", longname="Potential Temperature", & @@ -2715,7 +2715,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif if (CS%tv%S_is_absS) then vd_S = var_desc(name="abssalt", units="g kg-1", longname="Absolute Salinity", & - cmor_field_name="so", cmor_longname="Sea Water Salinity", & + cmor_field_name="absso", cmor_longname="Sea Water Absolute Salinity", & conversion=0.001*US%S_to_ppt) else vd_S = var_desc(name="salt", units="psu", longname="Salinity", & @@ -2799,10 +2799,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0 !allocate porous topography variables - allocate(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0 - allocate(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0 - allocate(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0 - allocate(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0 + allocate(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz), source=1.0) + allocate(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz), source=1.0) + allocate(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1), source=1.0) + allocate(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1), source=1.0) ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 8f872ceb15..e24d4954cb 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -122,16 +122,20 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) A_layer_prev(I,j) = A_layer endif ; enddo ; enddo ; enddo else - do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then - call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), A_layer, do_I(I,j)) - if (eta_u(I,j,K) - (eta_u(I,j,K+1)+dz_min) > 0.0) then - pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) + do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq + if (do_I(I,j)) then + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), A_layer, do_I(I,j)) + if (eta_u(I,j,K) - (eta_u(I,j,K+1)+dz_min) > 0.0) then + pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1))) + else + pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice + endif + A_layer_prev(I,j) = A_layer else - pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice + pbv%por_face_areaU(I,j,k) = 1.0 endif - A_layer_prev(I,j) = A_layer - endif ; enddo ; enddo ; enddo + enddo ; enddo ; enddo endif ! v-points @@ -154,16 +158,20 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt) A_layer_prev(i,J) = A_layer endif ; enddo ; enddo ; enddo else - do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then - call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), A_layer, do_I(i,J)) - if (eta_v(i,J,K) - (eta_v(i,J,K+1)+dz_min) > 0.0) then - pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) + do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie + if (do_I(i,J)) then + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), A_layer, do_I(i,J)) + if (eta_v(i,J,K) - (eta_v(i,J,K+1)+dz_min) > 0.0) then + pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1))) + else + pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice + endif + A_layer_prev(i,J) = A_layer else - pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice + pbv%por_face_areaV(i,J,k) = 1.0 endif - A_layer_prev(i,J) = A_layer - endif ; enddo ; enddo ; enddo + enddo ; enddo ; enddo endif if (CS%debug) then @@ -231,10 +239,14 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) endif ; enddo ; enddo ; enddo else - do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then - call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & - eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) - endif ; enddo ; enddo ; enddo + do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq + if (do_I(I,j)) then + call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), & + eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j)) + else + pbv%por_layer_widthU(I,j,K) = 1.0 + endif + enddo ; enddo ; enddo endif ! v-points @@ -249,10 +261,14 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt) eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) endif ; enddo ; enddo ; enddo else - do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then - call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & - eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) - endif ; enddo ; enddo ; enddo + do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie + if (do_I(i,J)) then + call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), & + eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J)) + else + pbv%por_layer_widthV(i,J,K) = 1.0 + endif + enddo ; enddo ; enddo endif if (CS%debug) then diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 5b7740230a..65e915705a 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -323,7 +323,6 @@ module MOM_variables end type BT_cont_type !> Container for grids modifying cell metric at porous barriers -! TODO: rename porous_barrier_type to porous_barrier_type type, public :: porous_barrier_type ! Each of the following fields has nz layers. real, allocatable :: por_face_areaU(:,:,:) !< fractional open area of U-faces [nondim] diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 52b1ebdcea..bac5b0fce9 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -450,11 +450,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) I_au = 0.0 ; if (asu1 + asu2 > 0.0) I_au = 1.0 / (asu1 + asu2) I_av = 0.0 ; if (asv1 + asv2 > 0.0) I_av = 1.0 / (asv1 + asv2) if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then - taux2 = (asu1 * sfc_state%taux_shelf(I-1,j)**2 + asu2 * sfc_state%taux_shelf(I,j)**2 ) * I_au - tauy2 = (asv1 * sfc_state%tauy_shelf(i,J-1)**2 + asv2 * sfc_state%tauy_shelf(i,J)**2 ) * I_av + taux2 = (((asu1 * (sfc_state%taux_shelf(I-1,j)**2)) + (asu2 * (sfc_state%taux_shelf(I,j)**2)) ) * I_au) + tauy2 = (((asv1 * (sfc_state%tauy_shelf(i,J-1)**2)) + (asv2 * (sfc_state%tauy_shelf(i,J)**2)) ) * I_av) endif - u2_av = (asu1 * sfc_state%u(I-1,j)**2 + asu2 * sfc_state%u(I,j)**2) * I_au - v2_av = (asv1 * sfc_state%v(i,J-1)**2 + asu2 * sfc_state%v(i,J)**2) * I_av + u2_av = (((asu1 * (sfc_state%u(I-1,j)**2)) + (asu2 * sfc_state%u(I,j)**2)) * I_au) + v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asu2 * sfc_state%v(i,J)**2)) * I_av) if ((taux2 + tauy2 > 0.0) .and. .not.CS%ustar_shelf_from_vel) then if (CS%ustar_max >= 0.0) then @@ -824,7 +824,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ISS%dhdt_shelf(i,j) = (ISS%h_shelf(i,j) - ISS%dhdt_shelf(i,j))*Itime_step enddo; enddo - call IS_dynamics_post_data(time_step, Time, CS%dCS, G) + call IS_dynamics_post_data(time_step, Time, CS%dCS, ISS, G) endif if (CS%shelf_mass_is_dynamic) & @@ -2603,7 +2603,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Ifull_time_step, dh_adott, dh_adott*0.0) call disable_averaging(CS%diag) - call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) + call IS_dynamics_post_data(full_time_step, Time, CS%dCS, ISS, G) end subroutine solo_step_ice_shelf !> Post_data calls for ice-sheet scalars @@ -2648,7 +2648,9 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh endif if (CS%id_f_adott > 0 .or. CS%id_f_adot > 0) then !floating only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - tmp(:,:) = dh_adott(:,:) - tmp(:,:) + do j=js,je ; do i=is,ie + tmp(i,j) = dh_adott(i,j) - tmp(i,j) + enddo; enddo call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) @@ -2710,7 +2712,9 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh endif if (CS%id_Ant_f_adott > 0 .or. CS%id_Ant_f_adot > 0) then !floating only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - tmp(:,:) = dh_adott(:,:) - tmp(:,:) + do j=js,je ; do i=is,ie + tmp(i,j) = dh_adott(i,j) - tmp(i,j) + enddo; enddo call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) if (CS%id_Ant_f_adott > 0) call post_scalar_data(CS%id_Ant_f_adott,val ,CS%diag) if (CS%id_Ant_f_adot > 0) call post_scalar_data(CS%id_Ant_f_adot ,val*Itime_step,CS%diag) @@ -2772,7 +2776,9 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh endif if (CS%id_Gr_f_adott > 0 .or. CS%id_Gr_f_adot > 0) then !floating only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - tmp(:,:) = dh_adott(:,:) - tmp(:,:) + do j=js,je ; do i=is,ie + tmp(i,j) = dh_adott(i,j) - tmp(i,j) + enddo; enddo call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) if (CS%id_Gr_f_adott > 0) call post_scalar_data(CS%id_Gr_f_adott,val ,CS%diag) if (CS%id_Gr_f_adot > 0) call post_scalar_data(CS%id_Gr_f_adot ,val*Itime_step,CS%diag) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 3ed262e5f3..9c7dda22de 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -89,9 +89,9 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: t_shelf => NULL() !< Vertically integrated temperature in the ice shelf/stream, !! on corner-points (B grid) [C ~> degC] real, pointer, dimension(:,:) :: tmask => NULL() !< A mask on tracer points that is 1 where there is ice. - real, pointer, dimension(:,:,:) :: ice_visc => NULL() !< Glen's law ice viscosity (Pa s), - !! in [R L2 T-1 ~> kg m-1 s-1]. - !! at either 1 (cell-centered) or 4 quadrature points per cell + real, pointer, dimension(:,:,:) :: ice_visc => NULL() !< Area and depth-integrated Glen's law ice viscosity + !! (Pa m3 s) in [R L4 Z T-1 ~> kg m2 s-1]. + !! at either 1 (cell-centered) or 4 quadrature points per cell real, pointer, dimension(:,:) :: AGlen_visc => NULL() !< Ice-stiffness parameter in Glen's law ice viscosity, !! often in [Pa-3 s-1] if n_Glen is 3. real, pointer, dimension(:,:) :: u_bdry_val => NULL() !< The zonal ice velocity at inflowing boundaries @@ -227,11 +227,17 @@ module MOM_ice_shelf_dynamics logical :: module_is_initialized = .false. !< True if this module has been initialized. !>@{ Diagnostic handles - integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, & - id_taudx_shelf = -1, id_taudy_shelf = -1, id_bed_elev = -1, & + integer :: id_u_shelf = -1, id_v_shelf = -1, id_shelf_speed, id_t_shelf = -1, & + id_taudx_shelf = -1, id_taudy_shelf = -1, id_taud_shelf = -1, id_bed_elev = -1, & id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, id_float_cond = -1, & id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1, & - id_sx_shelf = -1, id_sy_shelf = -1 + id_sx_shelf = -1, id_sy_shelf = -1, id_surf_slope_mag_shelf, & + id_duHdx = -1, id_dvHdy = -1, id_fluxdiv = -1, & + id_strainrate_xx = -1, id_strainrate_yy = -1, id_strainrate_xy = -1, & + id_pstrainrate_1 = -1, id_pstrainrate_2, & + id_devstress_xx = -1, id_devstress_yy = -1, id_devstress_xy = -1, & + id_pdevstress_1 = -1, id_pdevstress_2 = -1 + !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) !>@{ Diagnostic handles for debugging @@ -282,9 +288,9 @@ function quad_area (X, Y) ! | | ! 1 - 2 - p2 = (X(4)-X(1))**2 + (Y(4)-Y(1))**2 ; q2 = (X(3)-X(2))**2 + (Y(3)-Y(2))**2 - a2 = (X(3)-X(4))**2 + (Y(3)-Y(4))**2 ; c2 = (X(1)-X(2))**2 + (Y(1)-Y(2))**2 - b2 = (X(2)-X(4))**2 + (Y(2)-Y(4))**2 ; d2 = (X(3)-X(1))**2 + (Y(3)-Y(1))**2 + p2 = ( ((X(4)-X(1))**2) + ((Y(4)-Y(1))**2) ) ; q2 = ( ((X(3)-X(2))**2) + ((Y(3)-Y(2))**2) ) + a2 = ( ((X(3)-X(4))**2) + ((Y(3)-Y(4))**2) ) ; c2 = ( ((X(1)-X(2))**2) + ((Y(1)-Y(2))**2) ) + b2 = ( ((X(2)-X(4))**2) + ((Y(2)-Y(4))**2) ) ; d2 = ( ((X(3)-X(1))**2) + ((Y(3)-Y(1))**2) ) quad_area = .25 * sqrt(4*P2*Q2-(B2+D2-A2-C2)**2) end function quad_area @@ -552,7 +558,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ units="m", default=1.e-3, scale=US%m_to_Z) call get_param(param_file, mdl, "NONLIN_SOLVE_ERR_MODE", CS%nonlin_solve_err_mode, & "Choose whether nonlin error in vel solve is based on nonlinear "//& - "residual (1), relative change since last iteration (2), or change in norm (3)", default=1) + "residual (1), relative change since last iteration (2), or change in norm (3)", default=3) call get_param(param_file, mdl, "SHELF_MOVING_FRONT", CS%moving_shelf_front, & "Specify whether to advance shelf front (and calve).", & @@ -805,14 +811,20 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesB1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) + CS%id_shelf_speed = register_diag_field('ice_shelf_model','shelf_speed',CS%diag%axesB1, Time, & + 'speed of of ice shelf', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, & 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, & 'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_taud_shelf = register_diag_field('ice_shelf_model','taud_shelf',CS%diag%axesB1, Time, & + 'magnitude of driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) CS%id_sx_shelf = register_diag_field('ice_shelf_model','sx_shelf',CS%diag%axesB1, Time, & 'x-surface slope of ice', 'none') CS%id_sy_shelf = register_diag_field('ice_shelf_model','sy_shelf',CS%diag%axesB1, Time, & 'y-surface slope of ice', 'none') + CS%id_surf_slope_mag_shelf = register_diag_field('ice_shelf_model','surf_slope_mag_shelf', CS%diag%axesB1, Time, & + 'magnitude of surface slope of ice', 'none') CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, & 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, & @@ -830,6 +842,33 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) + CS%id_duHdx = register_diag_field('ice_shelf_model','duHdx',CS%diag%axesT1, Time, & + 'x-component of ice-sheet flux divergence', 'm yr-1', conversion=365.0*86400.0*US%Z_to_m*US%s_to_T) + CS%id_dvHdy = register_diag_field('ice_shelf_model','dvHdy',CS%diag%axesT1, Time, & + 'y-component of ice-sheet flux divergence', 'm yr-1', conversion=365.0*86400.0*US%Z_to_m*US%s_to_T) + CS%id_fluxdiv = register_diag_field('ice_shelf_model','fluxdiv',CS%diag%axesT1, Time, & + 'ice-sheet flux divergence', 'm yr-1', conversion=365.0*86400.0*US%Z_to_m*US%s_to_T) + CS%id_strainrate_xx = register_diag_field('ice_shelf_model','strainrate_xx',CS%diag%axesT1, Time, & + 'x-component of ice-shelf strain-rate', 'yr-1', conversion=365.0*86400.0*US%s_to_T) + CS%id_strainrate_yy = register_diag_field('ice_shelf_model','strainrate_yy',CS%diag%axesT1, Time, & + 'y-component of ice-shelf strain-rate', 'yr-1', conversion=365.0*86400.0*US%s_to_T) + CS%id_strainrate_xy = register_diag_field('ice_shelf_model','strainrate_xy',CS%diag%axesT1, Time, & + 'xy-component of ice-shelf strain-rate', 'yr-1', conversion=365.0*86400.0*US%s_to_T) + CS%id_pstrainrate_1 = register_diag_field('ice_shelf_model','pstrainrate_1',CS%diag%axesT1, Time, & + 'max principal horizontal ice-shelf strain-rate', 'yr-1', conversion=365.0*86400.0*US%s_to_T) + CS%id_pstrainrate_2 = register_diag_field('ice_shelf_model','pstrainrate_2',CS%diag%axesT1, Time, & + 'min principal horizontal ice-shelf strain-rate', 'yr-1', conversion=365.0*86400.0*US%s_to_T) + CS%id_devstress_xx = register_diag_field('ice_shelf_model','devstress_xx',CS%diag%axesT1, Time, & + 'x-component of ice-shelf deviatoric stress', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_devstress_yy = register_diag_field('ice_shelf_model','devstress_yy',CS%diag%axesT1, Time, & + 'y-component of ice-shelf deviatoric stress', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_devstress_xy = register_diag_field('ice_shelf_model','devstress_xy',CS%diag%axesT1, Time, & + 'xy-component of ice-shelf deviatoric stress', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_pdevstress_1 = register_diag_field('ice_shelf_model','pdevstress_1',CS%diag%axesT1, Time, & + 'max principal horizontal ice-shelf deviatoric stress', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_pdevstress_2 = register_diag_field('ice_shelf_model','pdevstress_2',CS%diag%axesT1, Time, & + 'min principal ice-shelf deviatoric stress', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + !Update these variables so that they are nonzero in case !IS_dynamics_post_data is called before update_ice_shelf if (CS%id_taudx_shelf>0 .or. CS%id_taudy_shelf>0) & @@ -898,10 +937,10 @@ function ice_time_step_CFL(CS, ISS, G) min_vel = (1.0e-12/(365.0*86400.0)) * G%US%m_s_to_L_T do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (ISS%hmask(i,j) == 1.0 .or. ISS%hmask(i,j)==3) then dt_local = 2.0*G%areaT(i,j) / & - ((G%dyCu(I,j) * max(abs(CS%u_shelf(I,J) + CS%u_shelf(I,j-1)), min_vel) + & - G%dyCu(I-1,j)* max(abs(CS%u_shelf(I-1,J)+ CS%u_shelf(I-1,j-1)), min_vel)) + & - (G%dxCv(i,J) * max(abs(CS%v_shelf(i,J) + CS%v_shelf(i-1,J)), min_vel) + & - G%dxCv(i,J-1)* max(abs(CS%v_shelf(i,J-1)+ CS%v_shelf(i-1,J-1)), min_vel))) + (((G%dyCu(I,j) * max(abs(CS%u_shelf(I,J) + CS%u_shelf(I,j-1)), min_vel)) + & + (G%dyCu(I-1,j)* max(abs(CS%u_shelf(I-1,J)+ CS%u_shelf(I-1,j-1)), min_vel))) + & + ((G%dxCv(i,J) * max(abs(CS%v_shelf(i,J) + CS%v_shelf(i-1,J)), min_vel)) + & + (G%dxCv(i,J-1)* max(abs(CS%v_shelf(i,J-1)+ CS%v_shelf(i-1,J-1)), min_vel)))) min_dt = min(min_dt, dt_local) endif ; enddo ; enddo ! i- and j- loops @@ -1035,45 +1074,70 @@ subroutine masked_var_grounded(G,CS,var,varout) end subroutine masked_var_grounded !> Ice shelf dynamics post_data calls -subroutine IS_dynamics_post_data(time_step, Time, CS, G) +subroutine IS_dynamics_post_data(time_step, Time, CS, ISS, G) real :: time_step !< Length of time for post data averaging [T ~> s]. type(time_type), intent(in) :: Time !< The current model time type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure + type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe + !! the ice-shelf state type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. - real, dimension(SZDIB_(G),SZDJB_(G)) :: taud_x, taud_y ! Pa] - real, dimension(SZDI_(G),SZDJ_(G)) :: ice_visc !< area-averaged vertically integrated ice viscosity + real, dimension(SZDIB_(G),SZDJB_(G)) :: taud_x, taud_y, taud ! area-averaged driving stress [R L2 T-2 ~> Pa] + real, dimension(SZDI_(G),SZDJ_(G)) :: ice_visc ! area-averaged vertically integrated ice viscosity !! [R L2 Z T-1 ~> Pa s m] - real, dimension(SZDI_(G),SZDJ_(G)) :: basal_tr !< area-averaged taub_beta field related to basal traction, + real, dimension(SZDI_(G),SZDJ_(G)) :: basal_tr ! area-averaged taub_beta field related to basal traction, !! [R L1 T-1 ~> Pa s m-1] + real, dimension(SZDIB_(G),SZDJB_(G)) :: surf_slope ! the surface slope of the ice shelf/sheet [nondim] + real, dimension(SZDIB_(G),SZDJB_(G)) :: ice_speed ! ice sheet flow speed [L T-1 ~> m s-1] + + integer :: i,j call enable_averages(time_step, Time, CS%diag) if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) + if (CS%id_shelf_speed > 0) then + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + ice_speed(I,J) = sqrt((CS%u_shelf(I,J)**2) + (CS%v_shelf(I,J)**2)) + enddo ; enddo + call post_data(CS%id_shelf_speed, ice_speed, CS%diag) + endif ! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf, CS%t_shelf, CS%diag) if (CS%id_taudx_shelf > 0) then - taud_x(:,:) = CS%taudx_shelf(:,:)*G%IareaBu(:,:) + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + taud_x(I,J) = CS%taudx_shelf(I,J)*G%IareaBu(I,J) + enddo ; enddo call post_data(CS%id_taudx_shelf, taud_x, CS%diag) endif if (CS%id_taudy_shelf > 0) then - taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaBu(:,:) + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + taud_y(I,J) = CS%taudy_shelf(I,J)*G%IareaBu(I,J) + enddo ; enddo call post_data(CS%id_taudy_shelf, taud_y, CS%diag) endif + if (CS%id_taud_shelf > 0) then + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + taud(I,J) = sqrt((CS%taudx_shelf(I,J)**2)+(CS%taudy_shelf(I,J)**2))*G%IareaBu(I,J) + enddo ; enddo + call post_data(CS%id_taud_shelf, taud, CS%diag) + endif if (CS%id_sx_shelf > 0) call post_data(CS%id_sx_shelf, CS%sx_shelf, CS%diag) if (CS%id_sy_shelf > 0) call post_data(CS%id_sy_shelf, CS%sy_shelf, CS%diag) + if (CS%id_surf_slope_mag_shelf > 0) then + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + surf_slope(I,J) = sqrt((CS%sx_shelf(I,J)**2)+(CS%sy_shelf(I,J)**2)) + enddo ; enddo + call post_data(CS%id_surf_slope_mag_shelf, surf_slope, CS%diag) + endif if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac, CS%diag) if (CS%id_float_cond > 0) call post_data(CS%id_float_cond, CS%float_cond, CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) if (CS%id_visc_shelf > 0) then - if (CS%visc_qps==4) then - ice_visc(:,:) = (0.25 * G%IareaT(:,:)) * & - ((CS%ice_visc(:,:,1) + CS%ice_visc(:,:,4)) + (CS%ice_visc(:,:,2) + CS%ice_visc(:,:,3))) - else - ice_visc(:,:) = CS%ice_visc(:,:,1)*G%IareaT(:,:) - endif + call ice_visc_diag(CS,G,ice_visc) call post_data(CS%id_visc_shelf, ice_visc, CS%diag) endif if (CS%id_taub > 0) then - basal_tr(:,:) = CS%basal_traction(:,:)*G%IareaT(:,:) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + basal_tr(i,j) = CS%basal_traction(i,j)*G%IareaT(i,j) + enddo ; enddo call post_data(CS%id_taub, basal_tr, CS%diag) endif if (CS%id_u_mask > 0) call post_data(CS%id_u_mask, CS%umask, CS%diag) @@ -1082,9 +1146,38 @@ subroutine IS_dynamics_post_data(time_step, Time, CS, G) if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask, CS%v_face_mask_bdry, CS%diag) ! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask, CS%tmask, CS%diag) + if (CS%id_duHdx > 0 .or. CS%id_dvHdy > 0 .or. CS%id_fluxdiv > 0 .or. & + CS%id_devstress_xx > 0 .or. CS%id_devstress_yy > 0 .or. CS%id_devstress_xy > 0 .or. & + CS%id_strainrate_xx > 0 .or. CS%id_strainrate_yy > 0 .or. CS%id_strainrate_xy > 0 .or. & + CS%id_pdevstress_1 > 0 .or. CS%id_pdevstress_2 > 0 .or. & + CS%id_pstrainrate_1 > 0 .or. CS%id_pstrainrate_2 > 0) then + call IS_dynamics_post_data_2(CS, ISS, G) + endif + call disable_averaging(CS%diag) end subroutine IS_dynamics_post_data +!> Calculate cell-centered, area-averaged, vertically integrated ice viscosity for diagnostics +subroutine ice_visc_diag(CS,G,ice_visc) + type(ice_shelf_dyn_CS), intent(in) :: CS !< The ice shelf dynamics control structure + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. + real, dimension(SZDI_(G),SZDJ_(G)), intent(out) :: ice_visc !< area-averaged vertically integrated ice viscosity + !! [R L2 Z T-1 ~> Pa s m] + integer :: i,j + + ice_visc(:,:)=0.0 + if (CS%visc_qps==4) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ice_visc(i,j) = (0.25 * G%IareaT(i,j)) * & + ((CS%ice_visc(i,j,1) + CS%ice_visc(i,j,4)) + (CS%ice_visc(i,j,2) + CS%ice_visc(i,j,3))) + enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ice_visc(i,j) = CS%ice_visc(i,j,1)*G%IareaT(i,j) + enddo ; enddo + endif +end subroutine ice_visc_diag + !> Writes the total ice shelf kinetic energy and mass to an ascii file subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure @@ -1154,8 +1247,8 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) KE_scale_factor = US%L_to_m**2 * (US%RZ_to_kg_m2 * US%L_T_to_m_s**2) do j=js,je ; do i=is,ie tmp1(i,j) = (KE_scale_factor * 0.03125) * (mass(i,j) * area(i,j)) * & - (((CS%u_shelf(I-1,J-1)+CS%u_shelf(I,J))+(CS%u_shelf(I,J-1)+CS%u_shelf(I-1,J)))**2 + & - ((CS%v_shelf(I-1,J-1)+CS%v_shelf(I,J))+(CS%v_shelf(I,J-1)+CS%v_shelf(I-1,J)))**2) + ((((CS%u_shelf(I-1,J-1)+CS%u_shelf(I,J))+(CS%u_shelf(I,J-1)+CS%u_shelf(I-1,J)))**2) + & + (((CS%v_shelf(I-1,J-1)+CS%v_shelf(I,J))+(CS%v_shelf(I,J-1)+CS%v_shelf(I-1,J)))**2)) enddo; enddo KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer) @@ -1447,7 +1540,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) err_init = 0 ; err_tempu = 0 ; err_tempv = 0 - do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB + do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB if (CS%umask(I,J) == 1) then err_tempu = ABS(Au(I,J) - taudx(I,J)) if (err_tempu >= err_init) err_init = err_tempu @@ -1529,7 +1622,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i err_max = 0 - do J=G%jscB,G%jecB ; do I=G%jscB,G%iecB + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB if (CS%umask(I,J) == 1) then err_tempu = ABS(Au(I,J) - taudx(I,J)) if (err_tempu >= err_max) err_max = err_tempu @@ -1556,7 +1649,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (CS%vmask(I,J) == 1) then err_tempv = MAX(ABS(v_last(I,J)-v_shlf(I,J)), err_tempu) if (err_tempv >= err_max) err_max = err_tempv - tempv = SQRT(v_shlf(I,J)**2 + tempu**2) + tempv = SQRT((v_shlf(I,J)**2) + (tempu**2)) endif if (tempv >= max_vel) max_vel = tempv enddo ; enddo @@ -2540,8 +2633,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G),CS%visc_qps), & intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's - !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form - !! and units depend on the basal law exponent. + !! flow law [R L4 Z T-1 ~> kg m2 s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice !! shelf is floating: 0 if floating, 1 if not @@ -2608,53 +2700,53 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, qp = 2*(jq-1)+iq !current quad point - uq = (u_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - u_shlf(I,J) * (xquad(iq) * xquad(jq))) + & - (u_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - u_shlf(I-1,J) * (xquad(3-iq) * xquad(jq))) + uq = ((u_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (u_shlf(I,J) * (xquad(iq) * xquad(jq)))) + & + ((u_shlf(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (u_shlf(I-1,J) * (xquad(3-iq) * xquad(jq)))) - vq = (v_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - v_shlf(I,J) * (xquad(iq) * xquad(jq))) + & - (v_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - v_shlf(I-1,J) * (xquad(3-iq) * xquad(jq))) + vq = ((v_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (v_shlf(I,J) * (xquad(iq) * xquad(jq)))) + & + ((v_shlf(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (v_shlf(I-1,J) * (xquad(3-iq) * xquad(jq)))) - ux = (u_shlf(I-1,J-1) * Phi(1,qp,i,j) + & - u_shlf(I,J) * Phi(7,qp,i,j)) + & - (u_shlf(I,J-1) * Phi(3,qp,i,j) + & - u_shlf(I-1,J) * Phi(5,qp,i,j)) + ux = ((u_shlf(I-1,J-1) * Phi(1,qp,i,j)) + & + (u_shlf(I,J) * Phi(7,qp,i,j))) + & + ((u_shlf(I,J-1) * Phi(3,qp,i,j)) + & + (u_shlf(I-1,J) * Phi(5,qp,i,j))) - vx = (v_shlf(I-1,J-1) * Phi(1,qp,i,j) + & - v_shlf(I,J) * Phi(7,qp,i,j)) + & - (v_shlf(I,J-1) * Phi(3,qp,i,j) + & - v_shlf(I-1,J) * Phi(5,qp,i,j)) + vx = ((v_shlf(I-1,J-1) * Phi(1,qp,i,j)) + & + (v_shlf(I,J) * Phi(7,qp,i,j))) + & + ((v_shlf(I,J-1) * Phi(3,qp,i,j)) + & + (v_shlf(I-1,J) * Phi(5,qp,i,j))) - uy = (u_shlf(I-1,J-1) * Phi(2,qp,i,j) + & - u_shlf(I,J) * Phi(8,qp,i,j)) + & - (u_shlf(I,J-1) * Phi(4,qp,i,j) + & - u_shlf(I-1,J) * Phi(6,qp,i,j)) + uy = ((u_shlf(I-1,J-1) * Phi(2,qp,i,j)) + & + (u_shlf(I,J) * Phi(8,qp,i,j))) + & + ((u_shlf(I,J-1) * Phi(4,qp,i,j)) + & + (u_shlf(I-1,J) * Phi(6,qp,i,j))) - vy = (v_shlf(I-1,J-1) * Phi(2,qp,i,j) + & - v_shlf(I,J) * Phi(8,qp,i,j)) + & - (v_shlf(I,J-1) * Phi(4,qp,i,j) + & - v_shlf(I-1,J) * Phi(6,qp,i,j)) + vy = ((v_shlf(I-1,J-1) * Phi(2,qp,i,j)) + & + (v_shlf(I,J) * Phi(8,qp,i,j))) + & + ((v_shlf(I,J-1) * Phi(4,qp,i,j)) + & + (v_shlf(I-1,J) * Phi(6,qp,i,j))) if (visc_qp4) qpv = qp !current quad point for viscosity do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & - ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & - (uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) + (((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & + ((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & - ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & - (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) + (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & + ((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) if (float_cond(i,j) == 0) then ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) + ((basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) + ((basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))) endif enddo ; enddo enddo ; enddo @@ -2732,13 +2824,13 @@ subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, uloc_arr(:,:,:,:) = 0.0; vloc_arr(:,:,:,:)=0.0 do j=1,nsub ; do i=1,nsub; do qy=1,2 ; do qx=1,2 - hloc = (Phisub(qx,qy,i,j,1,1)*H(1,1) + Phisub(qx,qy,i,j,2,2)*H(2,2)) + & - (Phisub(qx,qy,i,j,1,2)*H(1,2) + Phisub(qx,qy,i,j,2,1)*H(2,1)) + hloc = ((Phisub(qx,qy,i,j,1,1)*H(1,1)) + (Phisub(qx,qy,i,j,2,2)*H(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*H(1,2)) + (Phisub(qx,qy,i,j,2,1)*H(2,1))) if (dens_ratio * hloc - bathyT > 0) then - uloc_arr(qx,qy,i,j) = ((Phisub(qx,qy,i,j,1,1) * U(1,1) + Phisub(qx,qy,i,j,2,2) * U(2,2)) + & - (Phisub(qx,qy,i,j,1,2) * U(1,2) + Phisub(qx,qy,i,j,2,1) * U(2,1))) - vloc_arr(qx,qy,i,j) = ((Phisub(qx,qy,i,j,1,1) * V(1,1) + Phisub(qx,qy,i,j,2,2) * V(2,2)) + & - (Phisub(qx,qy,i,j,1,2) * V(1,2) + Phisub(qx,qy,i,j,2,1) * V(2,1))) + uloc_arr(qx,qy,i,j) = (((Phisub(qx,qy,i,j,1,1) * U(1,1)) + (Phisub(qx,qy,i,j,2,2) * U(2,2))) + & + ((Phisub(qx,qy,i,j,1,2) * U(1,2)) + (Phisub(qx,qy,i,j,2,1) * U(2,1)))) + vloc_arr(qx,qy,i,j) = (((Phisub(qx,qy,i,j,1,1) * V(1,1)) + (Phisub(qx,qy,i,j,2,2) * V(2,2))) + & + ((Phisub(qx,qy,i,j,1,2) * V(1,2)) + (Phisub(qx,qy,i,j,2,1) * V(2,1)))) endif enddo; enddo ; enddo ; enddo @@ -2817,12 +2909,10 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, !! (corner) points [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G),CS%visc_qps), & intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's - !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form - !! and units depend on the basal law exponent. + !! flow law [R L4 Z T-1 ~> kg m2 s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear !! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1]. - real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf @@ -2891,8 +2981,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, vy = 0. u_diag_qp(iphi,jphi,qp) = & - ice_visc(i,j,qpv) * ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & - (uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) + ice_visc(i,j,qpv) * (((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & + ((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) if (float_cond(i,j) == 0) then uq = xquad(ilq) * xquad(jlq) @@ -2909,8 +2999,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, uy = 0. v_diag_qp(iphi,jphi,qp) = & - ice_visc(i,j,qpv) * ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & - (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) + ice_visc(i,j,qpv) * (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & + ((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) if (float_cond(i,j) == 0) then vq = xquad(ilq) * xquad(jlq) @@ -2941,15 +3031,15 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground) - if (CS%umask(I-1,J-1) == 1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + sub_ground(1,1) * basal_trac(i,j) - if (CS%umask(I-1,J ) == 1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + sub_ground(1,2) * basal_trac(i,j) - if (CS%umask(I ,J-1) == 1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + sub_ground(2,1) * basal_trac(i,j) - if (CS%umask(I ,J ) == 1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + sub_ground(2,2) * basal_trac(i,j) + if (CS%umask(I-1,J-1) == 1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + (sub_ground(1,1) * basal_trac(i,j)) + if (CS%umask(I-1,J ) == 1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + (sub_ground(1,2) * basal_trac(i,j)) + if (CS%umask(I ,J-1) == 1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + (sub_ground(2,1) * basal_trac(i,j)) + if (CS%umask(I ,J ) == 1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + (sub_ground(2,2) * basal_trac(i,j)) - if (CS%vmask(I-1,J-1) == 1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + sub_ground(1,1) * basal_trac(i,j) - if (CS%vmask(I-1,J ) == 1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + sub_ground(1,2) * basal_trac(i,j) - if (CS%vmask(I ,J-1) == 1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + sub_ground(2,1) * basal_trac(i,j) - if (CS%vmask(I ,J ) == 1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + sub_ground(2,2) * basal_trac(i,j) + if (CS%vmask(I-1,J-1) == 1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + (sub_ground(1,1) * basal_trac(i,j)) + if (CS%vmask(I-1,J ) == 1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + (sub_ground(1,2) * basal_trac(i,j)) + if (CS%vmask(I ,J-1) == 1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + (sub_ground(2,1) * basal_trac(i,j)) + if (CS%vmask(I ,J ) == 1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + (sub_ground(2,2) * basal_trac(i,j)) endif endif ; enddo ; enddo @@ -2986,8 +3076,8 @@ subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd grnd_stat(:,:,:,:)=0 do j=1,nsub ; do i=1,nsub; do qy=1,2 ; do qx=1,2 - hloc = (Phisub(qx,qy,i,j,1,1)*H_node(1,1) + Phisub(qx,qy,i,j,2,2)*H_node(2,2)) + & - (Phisub(qx,qy,i,j,1,2)*H_node(1,2) + Phisub(qx,qy,i,j,2,1)*H_node(2,1)) + hloc = ((Phisub(qx,qy,i,j,1,1)*H_node(1,1)) + (Phisub(qx,qy,i,j,2,2)*H_node(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*H_node(1,2)) + (Phisub(qx,qy,i,j,2,1)*H_node(2,1))) if (dens_ratio * hloc - bathyT > 0) grnd_stat(qx,qy,i,j) = 1 enddo; enddo ; enddo ; enddo @@ -3006,6 +3096,140 @@ subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd end subroutine CG_diagonal_subgrid_basal +!> Post_data calls related to ice-sheet flux divergence, strain-rate, and deviatoric stress +subroutine IS_dynamics_post_data_2(CS, ISS, G) + type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure + type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe + !! the ice-shelf state + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. + real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m]. + real, dimension(SZDIB_(G),SZDJB_(G)) :: Hu ! Ice shelf u_flux at corners [Z L T-1 ~> m2 s-1]. + real, dimension(SZDIB_(G),SZDJB_(G)) :: Hv ! Ice shelf v_flux at corners [Z L T-1 ~> m2 s-1]. + real, dimension(SZDI_(G),SZDJ_(G)) :: Hux ! Ice shelf d(u_flux)/dx at cell centers [Z T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)) :: Hvy ! Ice shelf d(v_flux)/dy at cell centers [Z T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)) :: flux_div ! horizontal flux divergence div(uH) [Z T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G),3) :: strain_rate ! strain-rate components xx,yy, and xy [T-1 ~> s-1] + real, dimension(SZDI_(G),SZDJ_(G),2) :: p_strain_rate ! horizontal principal strain-rates [T-1 ~> s-1] + real, dimension(SZDI_(G),SZDJ_(G),3) :: dev_stress ! deviatoric stress components xx,yy, and xy [R L Z T-2 ~> Pa] + real, dimension(SZDI_(G),SZDJ_(G),2) :: p_dev_stress ! horizontal principal deviatoric stress [R L Z T-2 ~> Pa] + real, dimension(SZDI_(G),SZDJ_(G)) :: ice_visc ! area-averaged ice viscosity [R L2 T-1 ~> Pa s] + real :: p1,p2 ! Used to calculate strain-rate principal components [T-1 ~> s-1] + integer :: i, j + + !Allocate the gradient basis functions for 1 cell-centered quadrature point per cell + if (.not. associated(CS%PhiC)) then + allocate(CS%PhiC(1:8,G%isc:G%iec,G%jsc:G%jec), source=0.0) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call bilinear_shape_fn_grid_1qp(G, i, j, CS%PhiC(:,i,j)) + enddo; enddo + endif + + !Calculate flux divergence and its components + if (CS%id_duHdx > 0 .or. CS%id_dvHdy > 0 .or. CS%id_fluxdiv > 0) then + call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node, CS%min_h_shelf) + + Hu(:,:) = 0.0; Hv(:,:) = 0.0; Hux(:,:) = 0.0 ; Hvy(:,:) = 0.0 ; flux_div(:,:) = 0.0 + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + if (CS%umask(I,J) > 0) then + Hu(I,J) = (H_node(I,J) * CS%u_shelf(I,J)) + endif + if (CS%vmask(I,J) > 0) then + Hv(I,J) = (H_node(I,J) * CS%v_shelf(I,J)) + endif + enddo; enddo + + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 3)) then + !components of flux divergence at cell centers + Hux(i,j) = (((Hu(I-1,J-1) * CS%PhiC(1,i,j)) + (Hu(I,J ) * CS%PhiC(7,i,j))) + & + ((Hu(I-1,J ) * CS%PhiC(5,i,j)) + (Hu(I,J-1) * CS%PhiC(3,i,j)))) + + Hvy(i,j) = (((Hv(I-1,J-1) * CS%PhiC(2,i,j)) + (Hv(I,J ) * CS%PhiC(8,i,j))) + & + ((Hv(I-1,J ) * CS%PhiC(6,i,j)) + (Hv(I,J-1) * CS%PhiC(4,i,j)))) + flux_div(i,j) = Hux(i,j) + Hvy(i,j) + endif + enddo ; enddo + + if (CS%id_duHdx > 0) call post_data(CS%id_duHdx, Hux, CS%diag) + if (CS%id_dvHdy > 0) call post_data(CS%id_dvHdy, Hvy, CS%diag) + if (CS%id_fluxdiv > 0) call post_data(CS%id_fluxdiv, flux_div, CS%diag) + endif + + if (CS%id_devstress_xx > 0 .or. CS%id_devstress_yy > 0 .or. CS%id_devstress_xy > 0 .or. & + CS%id_strainrate_xx > 0 .or. CS%id_strainrate_yy > 0 .or. CS%id_strainrate_xy > 0 .or. & + CS%id_pdevstress_1 > 0 .or. CS%id_pdevstress_2 > 0 .or. & + CS%id_pstrainrate_1 > 0 .or. CS%id_pstrainrate_2 > 0) then + + strain_rate(:,:,:) = 0.0 + do j=G%jsc,G%jec ; do i=G%isc,G%iec + !strain-rates at cell centers + if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 3)) then + !strain_rate(:,:,1) = strain_rate_xx(:,:) = ux(:,:) + strain_rate(i,j,1) = (((CS%u_shelf(I-1,J-1) * CS%PhiC(1,i,j)) + (CS%u_shelf(I,J ) * CS%PhiC(7,i,j))) + & + ((CS%u_shelf(I-1,J ) * CS%PhiC(5,i,j)) + (CS%u_shelf(I,J-1) * CS%PhiC(3,i,j)))) + !strain_rate(:,:,2) = strain_rate_yy(:,:) = uy(:,:) + strain_rate(i,j,2) = (((CS%v_shelf(I-1,J-1) * CS%PhiC(2,i,j)) + (CS%v_shelf(I,J ) * CS%PhiC(8,i,j))) + & + ((CS%v_shelf(I-1,J ) * CS%PhiC(6,i,j)) + (CS%v_shelf(I,J-1) * CS%PhiC(4,i,j)))) + !strain_rate(:,:,3) = strain_rate_xy(:,:) = 0.5 * (uy(:,:) + vy(:,:)) + strain_rate(i,j,3) = 0.5 * ((((CS%u_shelf(I-1,J-1) * CS%PhiC(2,i,j)) + (CS%u_shelf(I,J ) * CS%PhiC(8,i,j))) + & + ((CS%u_shelf(I-1,J ) * CS%PhiC(6,i,j)) + (CS%u_shelf(I,J-1) * CS%PhiC(4,i,j))))+ & + (((CS%v_shelf(I-1,J-1) * CS%PhiC(1,i,j)) + (CS%v_shelf(I,J ) * CS%PhiC(7,i,j))) + & + ((CS%v_shelf(I-1,J ) * CS%PhiC(5,i,j)) + (CS%v_shelf(I,J-1) * CS%PhiC(3,i,j))))) + endif + enddo ; enddo + + + if (CS%id_strainrate_xx > 0) call post_data(CS%id_strainrate_xx, strain_rate(:,:,1), CS%diag) + if (CS%id_strainrate_yy > 0) call post_data(CS%id_strainrate_yy, strain_rate(:,:,2), CS%diag) + if (CS%id_strainrate_xy > 0) call post_data(CS%id_strainrate_xy, strain_rate(:,:,3), CS%diag) + + if (CS%id_pstrainrate_1 > 0 .or. CS%id_pstrainrate_2 > 0 .or. & + CS%id_pdevstress_1 > 0 .or. CS%id_pdevstress_2 > 0) then + p_strain_rate(:,:,:) = 0.0 + do j=G%jsc,G%jec ; do i=G%isc,G%iec + p1 = 0.5*( strain_rate(i,j,1) + strain_rate(i,j,2)) + p2 = sqrt( (( 0.5 * (strain_rate(i,j,1) - strain_rate(i,j,2)) )**2) + (strain_rate(i,j,3)**2) ) + p_strain_rate(i,j,1) = p1+p2 !Max horizontal principal strain-rate + p_strain_rate(i,j,2) = p1-p2 !Min horizontal principal strain-rate + enddo ; enddo + + if (CS%id_pstrainrate_1 > 0) call post_data(CS%id_pstrainrate_1, p_strain_rate(:,:,1), CS%diag) + if (CS%id_pstrainrate_2 > 0) call post_data(CS%id_pstrainrate_2, p_strain_rate(:,:,2), CS%diag) + endif + + if (CS%id_devstress_xx > 0 .or. CS%id_devstress_yy > 0 .or. CS%id_devstress_xy > 0 .or. & + CS%id_pdevstress_1 > 0 .or. CS%id_pdevstress_2 > 0) then + + call ice_visc_diag(CS,G,ice_visc) + + if (CS%id_devstress_xx > 0 .or. CS%id_devstress_yy > 0 .or. CS%id_devstress_xy > 0) then + dev_stress(:,:,:)=0.0 + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (ISS%h_shelf(i,j)>0) then + dev_stress(i,j,1) = 2*ice_visc(i,j)*strain_rate(i,j,1)/ISS%h_shelf(i,j) !deviatoric stress xx + dev_stress(i,j,2) = 2*ice_visc(i,j)*strain_rate(i,j,2)/ISS%h_shelf(i,j) !deviatoric stress yy + dev_stress(i,j,3) = 2*ice_visc(i,j)*strain_rate(i,j,3)/ISS%h_shelf(i,j) !deviatoric stress xy + endif + enddo; enddo + if (CS%id_devstress_xx > 0) call post_data(CS%id_devstress_xx, dev_stress(:,:,1), CS%diag) + if (CS%id_devstress_yy > 0) call post_data(CS%id_devstress_yy, dev_stress(:,:,2), CS%diag) + if (CS%id_devstress_xy > 0) call post_data(CS%id_devstress_xy, dev_stress(:,:,3), CS%diag) + endif + + if (CS%id_pdevstress_1 > 0 .or. CS%id_pdevstress_2 > 0) then + p_dev_stress(:,:,:)=0.0 + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (ISS%h_shelf(i,j)>0) then + p_dev_stress(i,j,1) = 2*ice_visc(i,j)*p_strain_rate(i,j,1)/ISS%h_shelf(i,j) !max horiz principal dev stress + p_dev_stress(i,j,2) = 2*ice_visc(i,j)*p_strain_rate(i,j,2)/ISS%h_shelf(i,j) !min horiz principal dev stress + endif + enddo; enddo + if (CS%id_pdevstress_1 > 0) call post_data(CS%id_pdevstress_1, p_dev_stress(:,:,1), CS%diag) + if (CS%id_pdevstress_2 > 0) call post_data(CS%id_pdevstress_2, p_dev_stress(:,:,2), CS%diag) + endif + endif + endif +end subroutine IS_dynamics_post_data_2 !> Update depth integrated viscosity, based on horizontal strain rates subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) @@ -3061,8 +3285,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! constant viscocity for debugging elseif (trim(CS%ice_viscosity_compute) == "OBS") then if (CS%AGlen_visc(i,j) >0) then - CS%ice_visc(i,j,1) = max(CS%AGlen_visc(i,j) * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)),& - CS%min_ice_visc) + CS%ice_visc(i,j,1) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & + max(CS%AGlen_visc(i,j) ,CS%min_ice_visc) endif ! Here CS%Aglen_visc(i,j) is the ice viscosity [Pa s ~> R L2 T-1] computed from obs and read from a file elseif (model_qp1) then @@ -3071,28 +3295,29 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) Visc_coef = (CS%AGlen_visc(i,j))**(-1./n_g) ! Units of Aglen_visc [Pa-(n_g) s-1] - ux = (u_shlf(I-1,J-1) * CS%PhiC(1,i,j) + & - u_shlf(I,J) * CS%PhiC(7,i,j)) + & - (u_shlf(I-1,J) * CS%PhiC(5,i,j) + & - u_shlf(I,J-1) * CS%PhiC(3,i,j)) - - vx = (v_shlf(I-1,J-1) * CS%PhiC(1,i,j) + & - v_shlf(I,J) * CS%PhiC(7,i,j)) + & - (v_shlf(I-1,J) * CS%PhiC(5,i,j) + & - v_shlf(I,J-1) * CS%PhiC(3,i,j)) - - uy = (u_shlf(I-1,J-1) * CS%PhiC(2,i,j) + & - u_shlf(I,J) * CS%PhiC(8,i,j)) + & - (u_shlf(I-1,J) * CS%PhiC(6,i,j) + & - u_shlf(I,J-1) * CS%PhiC(4,i,j)) - - vy = (v_shlf(I-1,J-1) * CS%PhiC(2,i,j) + & - v_shlf(I,J) * CS%PhiC(8,i,j)) + & - (v_shlf(I-1,J) * CS%PhiC(6,i,j) + & - v_shlf(I,J-1) * CS%PhiC(4,i,j)) - - CS%ice_visc(i,j,1) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & - (US%s_to_T**2 * ((ux**2 + vy**2) + (ux*vy + 0.25*(uy+vx)**2) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & + ux = ((u_shlf(I-1,J-1) * CS%PhiC(1,i,j)) + & + (u_shlf(I,J) * CS%PhiC(7,i,j))) + & + ((u_shlf(I-1,J) * CS%PhiC(5,i,j)) + & + (u_shlf(I,J-1) * CS%PhiC(3,i,j))) + + vx = ((v_shlf(I-1,J-1) * CS%PhiC(1,i,j)) + & + (v_shlf(I,J) * CS%PhiC(7,i,j))) + & + ((v_shlf(I-1,J) * CS%PhiC(5,i,j)) + & + (v_shlf(I,J-1) * CS%PhiC(3,i,j))) + + uy = ((u_shlf(I-1,J-1) * CS%PhiC(2,i,j)) + & + (u_shlf(I,J) * CS%PhiC(8,i,j))) + & + ((u_shlf(I-1,J) * CS%PhiC(6,i,j)) + & + (u_shlf(I,J-1) * CS%PhiC(4,i,j))) + + vy = ((v_shlf(I-1,J-1) * CS%PhiC(2,i,j)) + & + (v_shlf(I,J) * CS%PhiC(8,i,j))) + & + ((v_shlf(I-1,J) * CS%PhiC(6,i,j)) + & + (v_shlf(I,J-1) * CS%PhiC(4,i,j))) + + CS%ice_visc(i,j,1) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & + max(0.5 * Visc_coef * & + (US%s_to_T**2 * (((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) elseif (model_qp4) then !calculate viscosity at 4 quadrature points per cell @@ -3101,28 +3326,29 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) do iq=1,2 ; do jq=1,2 - ux = (u_shlf(I-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j)) + & - (u_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j)) - - vx = (v_shlf(I-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j) + & - v_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j)) + & - (v_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j) + & - v_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j)) - - uy = (u_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j)) + & - (u_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) - - vy = (v_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j) + & - v_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j)) + & - (v_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & - v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) - - CS%ice_visc(i,j,2*(jq-1)+iq) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & - (US%s_to_T**2 * ((ux**2 + vy**2) + (ux*vy + 0.25*(uy+vx)**2) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & + ux = ((u_shlf(I-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j)) + & + (u_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j))) + & + ((u_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j)) + & + (u_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j))) + + vx = ((v_shlf(I-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j)) + & + (v_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j))) + & + ((v_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j)) + & + (v_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j))) + + uy = ((u_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j)) + & + (u_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j))) + & + ((u_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j)) + & + (u_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j))) + + vy = ((v_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j)) + & + (v_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j))) + & + ((v_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j)) + & + (v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j))) + + CS%ice_visc(i,j,2*(jq-1)+iq) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & + max(0.5 * Visc_coef * & + (US%s_to_T**2 * (((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) enddo; enddo endif @@ -3181,7 +3407,7 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 - unorm = US%L_T_to_m_s * sqrt( (umid**2 + vmid**2) + (eps_min**2 * (G%dxT(i,j)**2 + G%dyT(i,j)**2)) ) + unorm = US%L_T_to_m_s * sqrt( ((umid**2) + (vmid**2)) + (eps_min**2 * (G%dxT(i,j)**2 + G%dyT(i,j)**2)) ) !Coulomb friction (Schoof 2005, Gagliardini et al 2007) if (CS%CoulombFriction) then @@ -3358,10 +3584,10 @@ subroutine bilinear_shape_functions (X, Y, Phi, area) do qpoint=1,4 - a = -X(1)*(1-yquad(qpoint)) + X(2)*(1-yquad(qpoint)) - X(3)*yquad(qpoint) + X(4)*yquad(qpoint) ! d(x)/d(x*) - b = -Y(1)*(1-yquad(qpoint)) + Y(2)*(1-yquad(qpoint)) - Y(3)*yquad(qpoint) + Y(4)*yquad(qpoint) ! d(y)/d(x*) - c = -X(1)*(1-xquad(qpoint)) - X(2)*xquad(qpoint) + X(3)*(1-xquad(qpoint)) + X(4)*xquad(qpoint) ! d(x)/d(y*) - d = -Y(1)*(1-xquad(qpoint)) - Y(2)*xquad(qpoint) + Y(3)*(1-xquad(qpoint)) + Y(4)*xquad(qpoint) ! d(y)/d(y*) + a = ((-X(1)*(1-yquad(qpoint)))+(X(4)*yquad(qpoint))) + ((X(2)*(1-yquad(qpoint)))-(X(3)*yquad(qpoint))) !d(x)/d(x*) + b = ((-Y(1)*(1-yquad(qpoint)))+(Y(4)*yquad(qpoint))) + ((Y(2)*(1-yquad(qpoint)))-(Y(3)*yquad(qpoint))) !d(y)/d(x*) + c = ((-X(1)*(1-xquad(qpoint)))+(X(4)*xquad(qpoint))) + ((-X(2)*xquad(qpoint))+(X(3)*(1-xquad(qpoint))))!d(x)/d(y*) + d = ((-Y(1)*(1-xquad(qpoint)))+(Y(4)*xquad(qpoint))) + ((-Y(2)*xquad(qpoint))+(Y(3)*(1-xquad(qpoint))))!d(y)/d(y*) do node=1,4 @@ -3379,8 +3605,8 @@ subroutine bilinear_shape_functions (X, Y, Phi, area) xexp = xquad(qpoint) endif - Phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c) - Phi(2*node,qpoint) = (-c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c) + Phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / ((a*d)-(b*c)) + Phi(2*node,qpoint) = (-c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / ((a*d)-(b*c)) enddo enddo @@ -3420,12 +3646,12 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) do qpoint=1,4 if (J>1) then - a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) + a = (G%dxCv(i,J-1) * (1-yquad(qpoint))) + (G%dxCv(i,J) * yquad(qpoint)) ! d(x)/d(x*) else a = G%dxCv(i,J) !* yquad(qpoint) ! d(x)/d(x*) endif if (I>1) then - d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) + d = (G%dyCu(I-1,j) * (1-xquad(qpoint))) + (G%dyCu(I,j) * xquad(qpoint)) ! d(y)/d(y*) else d = G%dyCu(I,j) !* xquad(qpoint) endif @@ -3707,7 +3933,7 @@ end subroutine update_velocity_masks !> Interpolate the ice shelf thickness from tracer point to nodal points, !! subject to a mask. subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node, min_h_shelf) - type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf. + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)), & @@ -3942,8 +4168,8 @@ subroutine ice_shelf_advect_temp_x(CS, G, time_step, hmask, h0, h_after_uflux) elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) - flux_diff = flux_diff + ABS(u_face) * G%dyCu(I-1,j)* time_step / G%areaT(i,j) * & - (stencil(-1) - phi * (stencil(-1)-stencil(0))/2) + flux_diff = flux_diff + ((ABS(u_face) * G%dyCu(I-1,j)* time_step / G%areaT(i,j)) * & + (stencil(-1) - (phi * (stencil(-1)-stencil(0))/2))) else ! h(i-1) is valid ! (o.w. flux would most likely be out of cell) @@ -3956,8 +4182,8 @@ subroutine ice_shelf_advect_temp_x(CS, G, time_step, hmask, h0, h_after_uflux) elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0)) - flux_diff = flux_diff - ABS(u_face) * G%dyCu(I-1,j) * time_step / G%areaT(i,j) * & - (stencil(0) - phi * (stencil(0)-stencil(-1))/2) + flux_diff = flux_diff - ((ABS(u_face) * G%dyCu(I-1,j) * time_step / G%areaT(i,j)) * & + (stencil(0) - (phi * (stencil(0)-stencil(-1))/2))) else flux_diff = flux_diff - ABS(u_face) * G%dyCu(I-1,j) * time_step / G%areaT(i,j) * stencil(0) @@ -3987,8 +4213,8 @@ subroutine ice_shelf_advect_temp_x(CS, G, time_step, hmask, h0, h_after_uflux) elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1)) - flux_diff = flux_diff + ABS(u_face) * G%dyCu(I,j) * time_step / G%areaT(i,j) * & - (stencil(1) - phi * (stencil(1)-stencil(0))/2) + flux_diff = flux_diff + ((ABS(u_face) * G%dyCu(I,j) * time_step / G%areaT(i,j)) * & + (stencil(1) - (phi * (stencil(1)-stencil(0))/2))) else ! h(i+1) is valid ! (o.w. flux would most likely be out of cell) @@ -4003,8 +4229,8 @@ subroutine ice_shelf_advect_temp_x(CS, G, time_step, hmask, h0, h_after_uflux) if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0)) - flux_diff = flux_diff - ABS(u_face) * G%dyCu(I,j) * time_step / G%areaT(i,j) * & - (stencil(0) - phi * (stencil(0)-stencil(1))/2) + flux_diff = flux_diff - ((ABS(u_face) * G%dyCu(I,j) * time_step / G%areaT(i,j)) * & + (stencil(0) - (phi * (stencil(0)-stencil(1))/2))) else ! h(i+1) is valid (o.w. flux would most likely be out of cell) but h(i+2) is not @@ -4108,8 +4334,8 @@ subroutine ice_shelf_advect_temp_y(CS, G, time_step, hmask, h_after_uflux, h_aft elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1)) - flux_diff = flux_diff + ABS(v_face) * G%dxCv(i,J-1) * time_step / G%areaT(i,j) * & - (stencil(-1) - phi * (stencil(-1)-stencil(0))/2) + flux_diff = flux_diff + ((ABS(v_face) * G%dxCv(i,J-1) * time_step / G%areaT(i,j)) * & + (stencil(-1) - (phi * (stencil(-1)-stencil(0))/2))) else ! h(j-1) is valid ! (o.w. flux would most likely be out of cell) @@ -4121,8 +4347,8 @@ subroutine ice_shelf_advect_temp_y(CS, G, time_step, hmask, h_after_uflux, h_aft if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0)) - flux_diff = flux_diff - ABS(v_face) * G%dxCv(i,J-1) * time_step / G%areaT(i,j) * & - (stencil(0) - phi * (stencil(0)-stencil(-1))/2) + flux_diff = flux_diff - ((ABS(v_face) * G%dxCv(i,J-1) * time_step / G%areaT(i,j)) * & + (stencil(0) - (phi * (stencil(0)-stencil(-1))/2))) else flux_diff = flux_diff - ABS(v_face) * G%dxCv(i,J-1) * time_step / G%areaT(i,j) * stencil(0) endif @@ -4148,8 +4374,8 @@ subroutine ice_shelf_advect_temp_y(CS, G, time_step, hmask, h_after_uflux, h_aft flux_diff = flux_diff + ABS(v_face) * G%dxCv(i,J) * time_step * stencil(1) / G%areaT(i,j) elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid phi = slope_limiter (stencil(1)-stencil(2), stencil(0)-stencil(1)) - flux_diff = flux_diff + ABS(v_face) * G%dxCv(i,J) * time_step / G%areaT(i,j) * & - (stencil(1) - phi * (stencil(1)-stencil(0))/2) + flux_diff = flux_diff + ((ABS(v_face) * G%dxCv(i,J) * time_step / G%areaT(i,j)) * & + (stencil(1) - (phi * (stencil(1)-stencil(0))/2))) else ! h(j+1) is valid ! (o.w. flux would most likely be out of cell) ! but h(j+2) is not @@ -4160,8 +4386,8 @@ subroutine ice_shelf_advect_temp_y(CS, G, time_step, hmask, h_after_uflux, h_aft if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid phi = slope_limiter (stencil(0)-stencil(-1), stencil(1)-stencil(0)) - flux_diff = flux_diff - ABS(v_face) * G%dxCv(i,J) * time_step / G%areaT(i,j) * & - (stencil(0) - phi * (stencil(0)-stencil(1))/2) + flux_diff = flux_diff - ((ABS(v_face) * G%dxCv(i,J) * time_step / G%areaT(i,j)) * & + (stencil(0) - (phi * (stencil(0)-stencil(1))/2))) else ! h(j+1) is valid ! (o.w. flux would most likely be out of cell) ! but h(j+2) is not diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index ea45e0cc2f..f54eb8a638 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -146,6 +146,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) endif ! Read sub-grid scale topography parameters at velocity points used for porous barrier calculation + ! TODO: The following routine call may eventually be merged as one of the CHANNEL_CONFIG options call get_param(PF, mdl, "SUBGRID_TOPO_AT_VEL", read_porous_file, & "If true, use variables from TOPO_AT_VEL_FILE as parameters for porous barrier.", & default=.False.) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 08f24f3a5c..9ebe8ae734 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -55,6 +55,7 @@ module MOM_MEKE logical :: initialized = .false. !< True if this control structure has been initialized. ! Parameters real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim] + real :: MEKE_bhFrCoeff!< Efficiency of conversion of ME into MEKE by the biharmonic dissipation [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] real :: MEKE_damping !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1]. @@ -126,8 +127,10 @@ module MOM_MEKE type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output !>@{ Diagnostic handles integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1 + integer :: id_src_adv = -1, id_src_mom_K4 = -1, id_src_btm_drag = -1 + integer :: id_src_GM = -1, id_src_mom_lp = -1, id_src_mom_bh = -1 integer :: id_Ub = -1, id_Ut = -1 - integer :: id_GM_src = -1, id_mom_src = -1, id_GME_snk = -1, id_decay = -1 + integer :: id_GM_src = -1, id_mom_src = -1, id_mom_src_bh = -1, id_GME_snk = -1, id_decay = -1 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 @@ -192,6 +195,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h depth_tot, & ! The depth of the water column [H ~> m or kg m-2]. src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3). MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1]. + src_adv, & ! The MEKE source/tendency from the horizontal advection of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_K4, & ! The MEKE source/tendency from the bihamornic of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_btm_drag, & ! The MEKE source/tendency from the bottom drag acting on MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). + damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1]. + MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1] drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1]. del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2]. @@ -222,9 +233,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: cdrag2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6] real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1] real :: mass_neglect ! A negligible mass [R Z ~> kg m-2]. - real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). + real :: damp_step ! Size of damping timestep relative to sdt [nondim] + real :: damp_rate ! The MEKE damping rate [T-1 ~> s-1]. + real :: damping ! The net damping of a field after sdt_damp [nondim] logical :: use_drag_rate ! Flag to indicate drag_rate is finite integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features @@ -254,6 +267,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%debug) then if (allocated(MEKE%mom_src)) & call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + if (allocated(MEKE%mom_src_bh)) & + call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GME_snk)) & call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GM_src)) & @@ -272,7 +287,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! With a depth-dependent (and possibly strong) damping, it seems ! advisable to use Strang splitting between the damping and diffusion. - sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt + damp_step = 1. + if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5 + sdt_damp = sdt * damp_step ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg] if (CS%MEKE_advection_factor>0.) then @@ -387,12 +404,21 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = CS%MEKE_BGsrc + src_adv(i,j) = 0. + src_mom_K4(i,j) = 0. + src_btm_drag(i,j) = 0. + src_GM(i,j) = 0. + src_mom_lp(i,j) = 0. + src_mom_bh(i,j) = 0. enddo ; enddo if (allocated(MEKE%mom_src)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) + src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) & + - (CS%MEKE_bhFrCoeff-CS%MEKE_FrCoeff)*I_mass(i,j)*MEKE%mom_src_bh(i,j) + src_mom_lp(i,j) = - CS%MEKE_FrCoeff*I_mass(i,j)*(MEKE%mom_src(i,j)-MEKE%mom_src_bh(i,j)) + src_mom_bh(i,j) = - CS%MEKE_bhFrCoeff*I_mass(i,j)*MEKE%mom_src_bh(i,j) enddo ; enddo endif @@ -414,6 +440,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) + src_GM(i,j) = -CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) enddo ; enddo endif endif @@ -433,6 +460,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie + MEKE_current(i,j) = MEKE%MEKE(i,j) MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j) enddo ; enddo @@ -453,12 +481,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! First stage of Strang splitting !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. + damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + + damping = 1. / (1. + sdt_damp * damp_rate) + + ! NOTE: MEKE%MEKE should use `damping` but we must preserve the existing + ! expression for bit reproducibility + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate) + MEKE_decay(i,j) = damp_rate * G%mask2dT(i,j) + + src_GM(i,j) = src_GM(i,j) * damping + src_mom_lp(i,j) = src_mom_lp(i,j) * damping + src_mom_bh(i,j) = src_mom_bh(i,j) * damping + + src_btm_drag(i,j) = - MEKE_current(i,j) * ( & + damp_step * (damp_rate * damping) & + ) + + ! Store the effective damping rate if sdt is split + if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) & + damp_rate_s1(i,j) = damp_rate * damping enddo ; enddo if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then @@ -528,6 +573,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h del4MEKE(i,j) = (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) + src_mom_K4(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & + ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & + (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) enddo ; enddo endif ! @@ -595,6 +643,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) + src_adv(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & + ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & + (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) enddo ; enddo endif ! MEKE_KH>0 @@ -608,25 +659,38 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Second stage of Strang splitting if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.0) then - if (sdt>sdt_damp) then - ! Recalculate the drag rate, since MEKE has changed. - if (use_drag_rate) then - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & - cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) - enddo ; enddo - endif + ! Recalculate the drag rate, since MEKE has changed. + if (use_drag_rate) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. - ! notice that the above line ensures a damping only if MEKE is positive, - ! while leaving MEKE unchanged if it is negative - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & + cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo endif + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. + ! notice that the above line ensures a damping only if MEKE is positive, + ! while leaving MEKE unchanged if it is negative + + damping = 1. / (1. + sdt_damp * damp_rate) + + ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve + ! the existing expression for bit reproducibility. + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate) + MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j) + + src_GM(i,j) = src_GM(i,j) * damping + src_mom_lp(i,j) = src_mom_lp(i,j) * damping + src_mom_bh(i,j) = src_mom_bh(i,j) * damping + src_adv(i,j) = src_adv(i,j) * damping + src_mom_K4(i,j) = src_mom_K4(i,j) * damping + + src_btm_drag(i,j) = -MEKE_current(i,j) * ( & + damp_step * damping * (damp_rate + damp_rate_s1(i,j)) & + ) + enddo ; enddo endif ! MEKE_KH>=0 if (CS%debug) then @@ -727,9 +791,16 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag) if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag) if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag) + if (CS%id_src_adv>0) call post_data(CS%id_src_adv, src_adv, CS%diag) + if (CS%id_src_mom_K4>0) call post_data(CS%id_src_mom_K4, src_mom_K4, CS%diag) + if (CS%id_src_btm_drag>0) call post_data(CS%id_src_btm_drag, src_btm_drag, CS%diag) + if (CS%id_src_GM>0) call post_data(CS%id_src_GM, src_GM, CS%diag) + if (CS%id_src_mom_lp>0) call post_data(CS%id_src_mom_lp, src_mom_lp, CS%diag) + if (CS%id_src_mom_bh>0) call post_data(CS%id_src_mom_bh, src_mom_bh, CS%diag) if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag) if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag) if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag) + if (CS%id_mom_src_bh>0) call post_data(CS%id_mom_src_bh, MEKE%mom_src_bh, CS%diag) if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag) if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag) if (CS%id_gamma_b>0) then @@ -1210,6 +1281,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME "The efficiency of the conversion of mean energy into "//& "MEKE. If MEKE_FRCOEFF is negative, this conversion "//& "is not used or calculated.", units="nondim", default=-1.0) + call get_param(param_file, mdl, "MEKE_BHFRCOEFF", CS%MEKE_bhFrCoeff, & + "The efficiency of the conversion of mean energy into "//& + "MEKE by the biharmonic dissipation. If MEKE_bhFRCOEFF is negative, this conversion "//& + "is not used or calculated.", units="nondim", default=-1.0) call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, & "The efficiency of the conversion of MEKE into mean energy "//& "by GME. If MEKE_GMECOEFF is negative, this conversion "//& @@ -1399,6 +1474,20 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (.not. allocated(MEKE%MEKE)) CS%id_Ut = -1 CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, & 'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + !add diagnostics for the terms in the MEKE budget + CS%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, Time, & + 'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', diag%axesT1, Time, & + 'MEKE energy source from the biharmonic of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, Time, & + 'MEKE energy source from the bottom drag acting on MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', diag%axesT1, Time, & + 'MEKE energy source from the thickness mixing (GM scheme)', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', diag%axesT1, Time, & + 'MEKE energy source from the Laplacian of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', diag%axesT1, Time, & + 'MEKE energy source from the biharmonic of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + !end CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, & 'MEKE decay rate', 's-1', conversion=US%s_to_T) CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, & @@ -1409,6 +1498,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME 'MEKE energy available from momentum', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (.not. allocated(MEKE%mom_src)) CS%id_mom_src = -1 + CS%id_mom_src_bh = register_diag_field('ocean_model', 'MEKE_mom_src_bh',diag%axesT1, Time, & + 'MEKE energy available from the biharmonic dissipation of momentum', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + if (.not. allocated(MEKE%mom_src_bh)) CS%id_mom_src_bh = -1 CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, & 'MEKE energy lost to GME backscatter', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) @@ -1742,7 +1835,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct ! Local variables - real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim] + real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_bhFrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim] real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim] logical :: Use_KH_in_MEKE logical :: useMEKE @@ -1754,6 +1847,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) ! Read these parameters to determine what should be in the restarts MEKE_GMcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff) MEKE_FrCoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff) + MEKE_bhFrCoeff = -1. ; call read_param(param_file,"MEKE_bhFRCOEFF",MEKE_bhFrCoeff) MEKE_GMEcoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff) MEKE_KhCoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff) MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) @@ -1770,8 +1864,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2) if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0) - if (MEKE_FrCoeff>=0. .or. MEKE_GMECoeff>=0.) & + if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) then allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0) + allocate(MEKE%mom_src_bh(isd:ied,jsd:jed), source=0.0) + endif + if (MEKE_FrCoeff<0.) MEKE_FrCoeff = 0. + if (MEKE_bhFrCoeff<0.) MEKE_bhFrCoeff = 0. if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0) if (MEKE_KhCoeff>=0.) then allocate(MEKE%Kh(isd:ied,jsd:jed), source=0.0) @@ -1817,6 +1915,7 @@ subroutine MEKE_end(MEKE) if (allocated(MEKE%Kh)) deallocate(MEKE%Kh) if (allocated(MEKE%GME_snk)) deallocate(MEKE%GME_snk) if (allocated(MEKE%mom_src)) deallocate(MEKE%mom_src) + if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh) if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src) if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE) end subroutine MEKE_end diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index e51f558ce3..a95578848d 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -11,6 +11,8 @@ module MOM_MEKE_types real, allocatable :: GM_src(:,:) !< MEKE source due to thickness mixing (GM) [R Z L2 T-3 ~> W m-2]. real, allocatable :: mom_src(:,:) !< MEKE source from lateral friction in the !! momentum equations [R Z L2 T-3 ~> W m-2]. + real, allocatable :: mom_src_bh(:,:) !< MEKE source from the biharmonic part of the lateral friction in the + !! momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: GME_snk(:,:) !< MEKE sink from GME backscatter in the momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: Kh(:,:) !< The MEKE-derived lateral mixing coefficient [L2 T-1 ~> m2 s-1]. real, allocatable :: Kh_diff(:,:) !< Uses the non-MEKE-derived thickness diffusion coefficient to diffuse diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 52a9bbefe1..56a359857f 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -50,6 +50,9 @@ module MOM_hor_visc !! limited to guarantee stability. logical :: better_bound_Kh !< If true, use a more careful bounding of the !! Laplacian viscosity to guarantee stability. + logical :: EY24_EBT_BS !! If true, use an equivalent barotropic backscatter + !! with a stabilizing kill switch in MEKE, + !< developed by Yankovsky et al. 2024 logical :: bound_Ah !< If true, the biharmonic coefficient is locally !! limited to guarantee stability. logical :: better_bound_Ah !< If true, use a more careful bounding of the @@ -60,6 +63,9 @@ module MOM_hor_visc !! the viscosity bounds to the theoretical maximum !! for stability without considering other terms [nondim]. !! The default is 0.8. + real :: KS_coef !< A nondimensional coefficient on the biharmonic viscosity that sets the + !! kill switch for backscatter. Default is 1.0 [nondim]. + real :: KS_timescale !< A timescale for computing CFL limit for turning off backscatter [T ~> s]. logical :: backscatter_underbound !< If true, the bounds on the biharmonic viscosity are allowed !! to increase where the Laplacian viscosity is negative (due to !! backscatter parameterizations) beyond the largest timestep-dependent @@ -145,6 +151,7 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & Kh_Max_xx, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1]. Ah_Max_xx, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1]. + Ah_Max_xx_KS, & !< The maximum permitted biharmonic viscosity for kill switch [L4 T-1 ~> m4 s-1]. n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points [nondim] n1n1_m_n2n2_h, & !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points [nondim] grid_sp_h2, & !< Harmonic mean of the squares of the grid [L2 ~> m2] @@ -163,6 +170,7 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & Kh_Max_xy, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1]. Ah_Max_xy, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1]. + Ah_Max_xy_KS, & !< The maximum permitted biharmonic viscosity for kill switch [L4 T-1 ~> m4 s-1]. n1n2_q, & !< Factor n1*n2 in the anisotropic direction tensor at q-points [nondim] n1n1_m_n2n2_q !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points [nondim] @@ -230,8 +238,13 @@ module MOM_hor_visc integer :: id_vort_xy_q = -1, id_div_xx_h = -1 integer :: id_sh_xy_q = -1, id_sh_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 + integer :: id_FrictWork_bh = -1, id_FrictWorkIntz_bh = -1 integer :: id_FrictWork_GME = -1 integer :: id_normstress = -1, id_shearstress = -1 + integer :: id_visc_limit_h = -1, id_visc_limit_q = -1 + integer :: id_visc_limit_h_flag = -1, id_visc_limit_q_flag = -1 + integer :: id_visc_limit_h_frac = -1, id_visc_limit_q_frac = -1 + integer :: id_BS_coeff_h = -1, id_BS_coeff_q = -1 !>@} end type hor_visc_CS @@ -313,6 +326,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [L2 T-2 ~> m2 s-2] bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] + FrictWorkIntz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R L2 T-3 ~> W m-2] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] @@ -321,7 +335,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] m_leithy, & ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2] Ah_sq, & ! The square of the biharmonic viscosity [L8 T-2 ~> m8 s-2] - htot ! The total thickness of all layers [H ~> m or kg m-2] + htot, & ! The total thickness of all layers [H ~> m or kg m-2] + str_xx_BS ! The diagonal term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2] real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] real :: boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] @@ -346,7 +361,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] ! This form guarantees that hq/hu < 4. - GME_effic_q ! The filtered efficiency of the GME terms at q points [nondim] + GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim] + str_xy_BS ! The cross term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2] real :: grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] real :: boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] @@ -356,6 +372,10 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1] sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1] GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1] + visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter [nondim] + visc_limit_q_flag, & ! determines whether backscatter is shut off [nondim] + visc_limit_q_frac, & ! determines how close backscatter is to shutting off [nondim] + BS_coeff_q, & ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1] ShSt ! A diagnostic array of shear stress [T-1 ~> s-1]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & KH_u_GME, & !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1] @@ -368,14 +388,19 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1] dz, & ! Height change across layers [Z ~> m] FrictWork, & ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2] + FrictWork_bh, & ! work done by the biharmonic MKE dissipation mechanisms [R L2 T-3 ~> W m-2] FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2] div_xx_h, & ! horizontal divergence [T-1 ~> s-1] sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] - NoSt ! A diagnostic array of normal stress [T-1 ~> s-1]. + NoSt, & ! A diagnostic array of normal stress [T-1 ~> s-1]. + BS_coeff_h ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] - GME_coeff_h ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] + GME_coeff_h, & ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] + visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter [nondim] + visc_limit_h_flag, & ! determines whether backscatter is shut off [nondim] + visc_limit_h_frac ! determines how close backscatter is to shutting off [nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & @@ -426,6 +451,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim] + real :: tmp ! Fields evaluated on active layers, used for constructing 3D stress fields ! NOTE: The position of these declarations can impact performance, due to the @@ -436,6 +462,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, real, dimension(SZIB_(G),SZJB_(G)) :: & Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1] Kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1] + Kh_BS, & ! Laplacian antiviscosity [L2 T-1 ~> m2 s-1] Shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1] vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1] vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (h or q) [L-1 T-1 ~> m-1 s-1] @@ -452,6 +479,13 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + visc_limit_h(:,:,:) = 0. + visc_limit_q(:,:,:) = 0. + visc_limit_h_flag(:,:,:) = 0. + visc_limit_q_flag(:,:,:) = 0. + visc_limit_h_frac(:,:,:) = 0. + visc_limit_q_frac(:,:,:) = 0. + m_leithy(:,:) = 0.0 ! Initialize if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then @@ -633,12 +667,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & - !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & + !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_bh, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont & !$OMP ) & !$OMP private( & - !$OMP i, j, k, n, & + !$OMP i, j, k, n, tmp, & !$OMP dudx, dudy, dvdx, dvdy, sh_xx, sh_xy, h_u, h_v, & !$OMP Del2u, Del2v, DY_dxBu, DX_dyBu, sh_xx_bt, sh_xy_bt, & !$OMP str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, & @@ -654,7 +688,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, & !$OMP vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, & !$OMP sh_xx_smooth, sh_xy_smooth, & - !$OMP vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy & + !$OMP vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy, & + !$OMP Kh_BS, str_xx_bs, str_xy_bs, bs_coeff_h, bs_coeff_q & + !$OMP ) & + !$OMP firstprivate( & + !$OMP visc_limit_h, visc_limit_h_frac, visc_limit_h_flag, & + !$OMP visc_limit_q, visc_limit_q_frac, visc_limit_q_flag & !$OMP ) do k=1,nz @@ -1140,15 +1179,15 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) enddo ; enddo - if (use_MEKE_Ku) then + if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then ! *Add* the MEKE contribution (which might be negative) if (CS%res_scale_MEKE) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k) enddo ; enddo else do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) enddo ; enddo endif endif @@ -1350,6 +1389,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif + if (CS%EY24_EBT_BS) then + do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh + tmp = CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j) + visc_limit_h(i,j,k) = tmp + visc_limit_h_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j)) + if (Ah(i,j) >= tmp) then + visc_limit_h_flag(i,j,k) = 1. + endif + enddo ; enddo + endif + if ((CS%id_Ah_h>0) .or. CS%debug .or. CS%use_Leithy) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh Ah_h(i,j,k) = Ah(i,j) @@ -1365,7 +1415,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (CS%id_grid_Re_Ah>0) then + if (CS%id_grid_Re_Ah > 0) then do j=js,je ; do i=is,ie KE = 0.125 * (((u(I,j,k) + u(I-1,j,k))**2) + ((v(i,J,k) + v(i,J-1,k))**2)) grid_Ah = max(Ah(i,j), CS%min_grid_Ah) @@ -1387,6 +1437,31 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx + ! Backscatter using MEKE + if (CS%EY24_EBT_BS) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_h_flag(i,j,k) > 0) then + Kh_BS(i,j) = 0. + else + Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) + endif + enddo ; enddo + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx_BS(i,j) = -Kh_BS(i,j) * sh_xx(i,j) + enddo ; enddo + + if (CS%id_BS_coeff_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + BS_coeff_h(i,j,k) = Kh_BS(i,j) + enddo ; enddo + endif + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx(i,j) = str_xx(i,j) + str_xx_BS(i,j) + enddo ; enddo + endif ! Backscatter + if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term do J=js-1,Jeq ; do I=is-1,Ieq @@ -1541,10 +1616,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired. - if (use_MEKE_Ku) then + if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then ! *Add* the MEKE contribution (might be negative) - Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & - (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn + Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn endif if (CS%anisotropic) & @@ -1671,6 +1748,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif + if (CS%EY24_EBT_BS) then + do J=js-1,Jeq ; do I=is-1,Ieq + tmp = CS%KS_coef *hrat_min(I,J) * CS%Ah_Max_xy_KS(I,J) + visc_limit_q(I,J,k) = tmp + visc_limit_q_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xy_KS(i,j)) + if (Ah(I,J) >= tmp) then + visc_limit_q_flag(I,J,k) = 1. + endif + enddo ; enddo + endif + ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points if (CS%use_Leithy) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1695,6 +1783,34 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif ! Get Ah at q points and biharmonic part of str_xy + ! Backscatter using MEKE + if (CS%EY24_EBT_BS) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_q_flag(I,J,k) > 0) then + Kh_BS(I,J) = 0. + else + Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) + endif + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy_BS(I,J) = -Kh_BS(I,J) * (sh_xy(I,J)) + enddo ; enddo + + if (CS%id_BS_coeff_q>0) then + do J=js-1,Jeq ; do I=is-1,Ieq + BS_coeff_q(I,J,k) = Kh_BS(I,J) + enddo ; enddo + endif + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = str_xy(I,J) + str_xy_BS(I,J) + enddo ; enddo + endif ! Backscatter + if (CS%use_GME) then ! The wider halo here is to permit one pass of smoothing without a halo update. do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 @@ -1795,34 +1911,42 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif if (find_FrictWork) then - if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie - ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion - FrictWork(i,j,k) = GV%H_to_RZ * ( & - ((str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & - - (str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & - + 0.25*(( (str_xy(I,J) * & - (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & - + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & - + (str_xy(I-1,J-1) * & - (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & - + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & - + ( (str_xy(I-1,J) * & - (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & - + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & - + (str_xy(I,J-1) * & - (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & - + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) + if (CS%FrictWork_bug) then + do j=js,je ; do i=is,ie + ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork(i,j,k) = 0 + else + FrictWork(i,j,k) = GV%H_to_RZ * ( & + ((str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & + - (str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & + + 0.25*(( (str_xy(I,J) * & + (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & + + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & + + (str_xy(I-1,J-1) * & + (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & + + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & + + ( (str_xy(I-1,J) * & + (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & + + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & + + (str_xy(I,J-1) * & + (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & + + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) + endif enddo ; enddo else ; do j=js,je ; do i=is,ie - FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork(i,j,k) = 0 + else + FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((str_xx(i,j)*CS%dy2h(i,j) * ( & (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) & - (str_xx(i,j)*CS%dx2h(i,j) * ( & (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) & - + (0.25*(((str_xy(I,J)*( & + + (0.25*(((str_xy(I,J)*( & (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) & - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) & + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) & @@ -1842,10 +1966,75 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) + + endif + enddo ; enddo ; endif + endif + + if (CS%id_FrictWork_bh>0 .or. CS%id_FrictWorkIntz_bh > 0 .or. allocated(MEKE%mom_src_bh)) then + if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie + ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork_bh(i,j,k) = 0 + else + ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion !cyc + FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & + (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + + 0.25*((bhstr_xy(I,J) * & + ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & + + bhstr_xy(I-1,J-1) * & + ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & + + (bhstr_xy(I-1,J) * & + ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & + + bhstr_xy(I,J-1) * & + ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + endif + enddo ; enddo + else ; do j=js,je ; do i=is,ie + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork_bh(i,j,k) = 0 + else + ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) + FrictWork_bh(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + ((bhstr_xx(i,j)*CS%dy2h(i,j) * ( & + (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & + - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) & + - (bhstr_xx(i,j)*CS%dx2h(i,j) * ( & + (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & + - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) & + + (0.25*(((bhstr_xy(I,J)*( & + (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) & + - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) & + + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) & + - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) & + +(bhstr_xy(I-1,J-1)*( & + (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) & + - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) & + + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) & + - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) & + +((bhstr_xy(I-1,J)*( & + (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) & + - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) & + + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & + - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) & + +(bhstr_xy(I,J-1)*( & + (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & + - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & + - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) + endif enddo ; enddo ; endif endif + if (CS%use_GME) then if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) @@ -1905,6 +2094,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (k==1) then do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = 0. + MEKE%mom_src_bh(i,j) = 0. enddo ; enddo if (allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie @@ -1937,27 +2127,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + GV%H_to_RZ * ( & - (((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & - - ((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & - + 0.25*( (((str_xy(I,J)-RoScl*bhstr_xy(I,J)) * & - (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & - + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & - + ((str_xy(I-1,J-1)-RoScl*bhstr_xy(I-1,J-1)) * & - (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & - + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & - + (((str_xy(I-1,J)-RoScl*bhstr_xy(I-1,J)) * & - (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & - + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & - + ((str_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & - (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & - + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) + MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + (FrictWork(i,j,k) - RoScl*FrictWork_bh(i,j,k)) + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + & + (FrictWork_bh(i,j,k) - RoScl*FrictWork_bh(i,j,k)) enddo ; enddo - endif ! MEKE%backscatter_Ro_c + else do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) enddo ; enddo + endif ! MEKE%backscatter_Ro_c if (CS%use_GME .and. allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie @@ -1975,6 +2155,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag) if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag) if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag) + if (CS%id_FrictWork_bh>0) call post_data(CS%id_FrictWork_bh, FrictWork_bh, CS%diag) if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) @@ -1994,6 +2175,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_dudy_bt > 0) call post_data(CS%id_dudy_bt, dudy_bt, CS%diag) if (CS%id_dvdx_bt > 0) call post_data(CS%id_dvdx_bt, dvdx_bt, CS%diag) endif + if (CS%id_visc_limit_h>0) call post_data(CS%id_visc_limit_h, visc_limit_h, CS%diag) + if (CS%id_visc_limit_q>0) call post_data(CS%id_visc_limit_q, visc_limit_q, CS%diag) + if (CS%id_visc_limit_h_frac>0) call post_data(CS%id_visc_limit_h_frac, visc_limit_h_frac, CS%diag) + if (CS%id_visc_limit_q_frac>0) call post_data(CS%id_visc_limit_q_frac, visc_limit_q_frac, CS%diag) + if (CS%id_visc_limit_h_flag>0) call post_data(CS%id_visc_limit_h_flag, visc_limit_h_flag, CS%diag) + if (CS%id_visc_limit_q_flag>0) call post_data(CS%id_visc_limit_q_flag, visc_limit_q_flag, CS%diag) + + if (CS%EY24_EBT_BS) then + if (CS%id_BS_coeff_h>0) call post_data(CS%id_BS_coeff_h, BS_coeff_h, CS%diag) + if (CS%id_BS_coeff_q>0) call post_data(CS%id_BS_coeff_q, BS_coeff_q, CS%diag) + endif if (CS%debug) then if (CS%Laplacian) then @@ -2014,6 +2206,16 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag) endif + if (CS%id_FrictWorkIntz_bh > 0) then + do j=js,je + do i=is,ie ; FrictWorkIntz_bh(i,j) = FrictWork_bh(i,j,1) ; enddo + do k=2,nz ; do i=is,ie + FrictWorkIntz_bh(i,j) = FrictWorkIntz_bh(i,j) + FrictWork_bh(i,j,k) + enddo ; enddo + enddo + call post_data(CS%id_FrictWorkIntz_bh, FrictWorkIntz_bh, CS%diag) + endif + if (present(ADp)) then ! Diagnostics of the fractional thicknesses times momentum budget terms ! 3D diagnostics of hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option. @@ -2206,8 +2408,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, the Laplacian coefficient is locally limited "//& "to be stable with a better bounding than just BOUND_KH.", & default=CS%bound_Kh, do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "EY24_EBT_BS", CS%EY24_EBT_BS, & + "If true, use the the backscatter scheme (EBT mode with kill switch)"//& + "developed by Yankovsky et al. (2024). ", & + default=.false., do_not_log=.not.CS%Laplacian) if (.not.CS%Laplacian) CS%bound_Kh = .false. if (.not.CS%Laplacian) CS%better_bound_Kh = .false. + if (.not.(CS%Laplacian.and.use_MEKE)) CS%EY24_EBT_BS = .false. call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& "horizontal viscosity.", default=.false., & @@ -2358,6 +2565,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "viscosity bounds to the theoretical maximum for "//& "stability without considering other terms.", units="nondim", & default=0.8, do_not_log=.not.(CS%better_bound_Ah .or. CS%better_bound_Kh)) + call get_param(param_file, mdl, "KILL_SWITCH_COEF", CS%KS_coef, & + "A nondimensional coefficient on the biharmonic viscosity that "// & + "sets the kill switch for backscatter. Default is 1.0.", units="nondim", & + default=1.0, do_not_log=.not.(CS%EY24_EBT_BS)) call get_param(param_file, mdl, "NOSLIP", CS%no_slip, & "If true, no slip boundary conditions are used; otherwise "//& "free slip boundary conditions are assumed. The "//& @@ -2423,6 +2634,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) fail_if_missing=.true.) Idt = 1.0 / dt endif + call get_param(param_file, mdl, "KILL_SWITCH_TIMESCALE", CS%KS_timescale, & + "A timescale for computing the CFL limit for viscosity "// & + "that determines when backscatter is shut off. Default is DT.", & + default= dt , units="s", scale=US%s_to_T, do_not_log=.not.(CS%EY24_EBT_BS)) + if (CS%no_slip .and. CS%biharmonic) & call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// & "at the same time in MOM.") @@ -2446,11 +2662,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0 ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0 ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0 - if (CS%bound_Kh .or. CS%better_bound_Kh) then + if (CS%bound_Kh .or. CS%better_bound_Kh .or. CS%EY24_EBT_BS) then ALLOC_(CS%Kh_Max_xx(Isd:Ied,Jsd:Jed)) ; CS%Kh_Max_xx(:,:) = 0.0 ALLOC_(CS%Kh_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_Max_xy(:,:) = 0.0 endif - if (CS%Smagorinsky_Kh) then + if (CS%Smagorinsky_Kh .or. CS%EY24_EBT_BS) then ALLOC_(CS%Laplac2_const_xx(isd:ied,jsd:jed)) ; CS%Laplac2_const_xx(:,:) = 0.0 ALLOC_(CS%Laplac2_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac2_const_xy(:,:) = 0.0 endif @@ -2508,6 +2724,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%Ah_Max_xx(isd:ied,jsd:jed)) ; CS%Ah_Max_xx(:,:) = 0.0 ALLOC_(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy(:,:) = 0.0 endif + if (CS%EY24_EBT_BS) then + ALLOC_(CS%Ah_Max_xx_KS(isd:ied,jsd:jed)) ; CS%Ah_Max_xx_KS(:,:) = 0.0 + ALLOC_(CS%Ah_Max_xy_KS(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy_KS(:,:) = 0.0 + endif if (CS%Smagorinsky_Ah) then ALLOC_(CS%Biharm_const_xx(isd:ied,jsd:jed)) ; CS%Biharm_const_xx(:,:) = 0.0 ALLOC_(CS%Biharm_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const_xy(:,:) = 0.0 @@ -2766,8 +2986,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) (CS%DX_dyT(i,j)*((G%IdxCv(i,J)*v0v(i,J)) + (G%IdxCv(i,J-1)*v0v(i,J-1))))) * & max(G%IdxCv(i,J)*G%IareaCv(i,J), G%IdxCv(i,J-1)*G%IareaCv(i,J-1)) ) ) CS%Ah_Max_xx(I,J) = 0.0 - if (denom > 0.0) & + if (denom > 0.0) then CS%Ah_Max_xx(I,J) = CS%bound_coef * 0.5 * Idt / denom + if (CS%EY24_EBT_BS) then + CS%Ah_Max_xx_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom) + endif + endif + enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq denom = max( & @@ -2780,8 +3005,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) (CS%DY_dxBu(I,J)*((v0v(i+1,J)*G%IdyCv(i+1,J)) + (v0v(i,J)*G%IdyCv(i,J))))) * & max(G%IdyCv(i,J)*G%IareaCv(i,J), G%IdyCv(i+1,J)*G%IareaCv(i+1,J)) ) ) CS%Ah_Max_xy(I,J) = 0.0 - if (denom > 0.0) & + if (denom > 0.0) then CS%Ah_Max_xy(I,J) = CS%bound_coef * 0.5 * Idt / denom + if (CS%EY24_EBT_BS) then + CS%Ah_Max_xy_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom) + endif + endif + enddo ; enddo if (CS%debug) then call hchksum(CS%Ah_Max_xx, "Ah_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T) @@ -2880,6 +3110,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T) CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, & 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim') + CS%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim') + CS%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim') + CS%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') + CS%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') if (CS%id_grid_Re_Ah > 0) & ! Compute the smallest biharmonic viscosity capable of modifying the @@ -2931,6 +3173,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) endif + + if (CS%EY24_EBT_BS) then + CS%id_BS_coeff_h = register_diag_field('ocean_model', 'BS_coeff_h', diag%axesTL, Time, & + 'Backscatter coefficient at h points', 'm2 s-1') + CS%id_BS_coeff_q = register_diag_field('ocean_model', 'BS_coeff_q', diag%axesBL, Time, & + 'Backscatter coefficient at q points', 'm2 s-1') + endif + CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& 'Integral work done by lateral friction terms. If GME is turned on, this '//& 'includes the GME contribution.', & @@ -2941,6 +3191,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) cmor_field_name='dispkexyfo', & cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',& cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction') + CS%id_FrictWork_bh = register_diag_field('ocean_model','FrictWork_bh',diag%axesTL,Time,& + 'Integral work done by the biharmonic lateral friction terms.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + CS%id_FrictWorkIntz_bh = register_diag_field('ocean_model','FrictWorkIntz_bh',diag%axesT1,Time,& + 'Depth integrated work done by the biharmonic lateral friction', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) end subroutine hor_visc_init @@ -3171,6 +3427,9 @@ subroutine hor_visc_end(CS) if (CS%bound_Ah) then DEALLOC_(CS%Ah_Max_xx) ; DEALLOC_(CS%Ah_Max_xy) endif + if (CS%EY24_EBT_BS) then + DEALLOC_(CS%Ah_Max_xx_KS) ; DEALLOC_(CS%Ah_Max_xy_KS) + endif if (CS%Smagorinsky_Ah) then DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy) endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7ff7f187fc..e3979fe35a 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -113,6 +113,9 @@ module MOM_lateral_mixing_coeffs real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim] real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim] real, allocatable :: ebt_struct(:,:,:) !< Vertical structure function to scale diffusivities with [nondim] + real, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim] + real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0. + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & Laplac3_const_u !< Laplacian metric-dependent constants [L3 ~> m3] @@ -228,7 +231,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2] integer :: power_2 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - integer :: i, j + integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -238,7 +241,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) if (CS%calculate_cg1) then if (.not. allocated(CS%cg1)) call MOM_error(FATAL, & "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.") - if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then + if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) then if (.not. allocated(CS%ebt_struct)) call MOM_error(FATAL, & "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.") if (CS%Resoln_use_ebt) then @@ -258,6 +261,11 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) call do_group_pass(CS%pass_cg1, G%Domain) endif + if (CS%BS_EBT_power>0.) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%BS_struct(i,j,k) = CS%ebt_struct(i,j,k)**CS%BS_EBT_power + enddo ; enddo ; enddo + endif ! Calculate and store the ratio between deformation radius and grid-spacing ! at h-points [nondim]. @@ -1249,6 +1257,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, uses the equivalent barotropic wave speed instead "//& "of first baroclinic wave for calculating the resolution fn.",& default=.false.) + call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, & + "Power to raise EBT vertical structure to when backscatter "// & + "has vertical structure.", units="nondim", default=0.0) call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, & "If true, uses the equivalent barotropic structure "//& "as the vertical structure of thickness diffusivity.",& @@ -1274,7 +1285,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) default=1.0e-17, units="s-1", scale=US%T_to_s) call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, & default=.false., do_not_log=.true.) - CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct + CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct & + .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0. CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE ! Indicate whether to calculate the Eady growth rate CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.) @@ -1298,7 +1310,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.") endif - if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & + .or. CS%BS_EBT_power>0.) then in_use = .true. call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, & "The depth below which N2 is monotonized to avoid stratification "//& @@ -1307,6 +1320,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) units="m", default=-1.0, scale=GV%m_to_H) allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif + allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0) + CS%BS_struct(:,:,:) = 1.0 if (CS%use_stored_slopes) then if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then @@ -1655,8 +1670,9 @@ end subroutine VarMix_init subroutine VarMix_end(CS) type(VarMix_CS), intent(inout) :: CS - if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) & + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) & deallocate(CS%ebt_struct) + if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) if (CS%use_stored_slopes) then deallocate(CS%slope_x)