diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 9fa37f9dbb..9448553a0f 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -702,8 +702,8 @@ subroutine ALE_remap_scalar(CS, G, nk_src, h_src, s_src, h_dst, s_dst, all_cells if (present(old_remap)) use_remapping_core_w = old_remap n_points = nk_src -!$OMP parallel default(none) shared(CS,G,h_src,s_src,h_dst,s_dst,use_remapping_core_w & -!$OMP ,ignore_vanished_layers, nk_src,dx ) & +!$OMP parallel default(none) shared(CS,G,h_src,s_src,h_dst,s_dst & +!$OMP ,ignore_vanished_layers, use_remapping_core_w, nk_src,dx ) & !$OMP firstprivate(n_points) !$OMP do do j = G%jsc,G%jec diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 641bc8e072..c5432934b6 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -133,11 +133,11 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & !$OMP parallel do default(none) shared(nz,is,ie,js,je,use_EOS,G,GV,pres,T,S, & !$OMP IsdB,tv,h,h_neglect,e,dz_neglect, & -!$OMP h_neglect2,present_N2_u,G_Rho0,N2_u) & +!$OMP h_neglect2,present_N2_u,G_Rho0,N2_u,slope_x) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & -!$OMP drdx,mag_grad2,Slope,slope2_Ratio,slope_x) +!$OMP drdx,mag_grad2,Slope,slope2_Ratio) do j=js,je ; do K=nz,2,-1 if (.not.(use_EOS)) then drdiA = 0.0 ; drdiB = 0.0 @@ -220,11 +220,11 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & ! Calculate the meridional isopycnal slope. !$OMP parallel do default(none) shared(nz,is,ie,js,je,use_EOS,G,GV,pres,T,S, & !$OMP IsdB,tv,h,h_neglect,e,dz_neglect, & -!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v) & +!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & -!$OMP drdy,mag_grad2,Slope,slope2_Ratio,slope_y) +!$OMP drdy,mag_grad2,Slope,slope2_Ratio) do j=js-1,je ; do K=nz,2,-1 if (.not.(use_EOS)) then drdjA = 0.0 ; drdjB = 0.0 diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 5069788f4e..eb1d15fced 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -374,9 +374,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & reset_diags = .false. ! This is the second call to mixedlayer. if (reset_diags) then -!$OMP parallel default(none) shared(is,ie,js,je,CS) +!!OMP parallel default(none) shared(is,ie,js,je,CS) if (CS%TKE_diagnostics) then -!$OMP do +!!OMP do do j=js,je ; do i=is,ie CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_MKE(i,j) = 0.0 CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_forcing(i,j) = 0.0 @@ -384,37 +384,37 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & CS%diag_TKE_conv_decay(i,j) = 0.0 enddo ; enddo endif -!$OMP end parallel +!!OMP end parallel endif -!$OMP parallel do default(none) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, & -!$OMP CS,G,GV,sfc_connected,fluxes,IdtdR0, & -!$OMP TKE_forced,debug,H_neglect,dSV_dT, & -!$OMP dSV_dS,I_dtmrho,C1_3,h_tt_min,vonKar, & -!$OMP max_itt,Kd_int) & -!$OMP private(h,u,v,T,S,Kd,mech_TKE_k,conv_PErel_k, & -!$OMP U_Star,absf,mech_TKE,conv_PErel,nstar_k, & -!$OMP h_sum,I_hs,h_bot,hb_hs,T0,S0,num_itts, & -!$OMP pres,dMass,dPres,dT_to_dPE,dS_to_dPE, & -!$OMP dT_to_dColHt,dS_to_dColHt,Kddt_h, & -!$OMP b_den_1,dT_to_dPE_a,dT_to_dColHt_a, & -!$OMP dS_to_dColHt_a,htot,uhtot,vhtot, & -!$OMP Idecay_len_TKE,exp_kh,nstar_FC,tot_TKE, & -!$OMP TKE_reduc,dTe_t2,dSe_t2,dTe,dSe,dt_h, & -!$OMP Convectively_stable,sfc_disconnect,b1, & -!$OMP c1,dT_km1_t2,dS_km1_t2,dTe_term, & -!$OMP dSe_term,MKE2_Hharm,vstar,h_tt, & -!$OMP Kd_guess0,Kddt_h_g0,dPEc_dKd_Kd0, & -!$OMP PE_chg_max,dPEa_dKd_g0,PE_chg_g0, & -!$OMP MKE_src,dPE_conv,Kddt_h_max,Kddt_h_min, & -!$OMP TKE_left_max,TKE_left_min,Kddt_h_guess, & -!$OMP TKE_left_itt,dPEa_dKd_itt,PE_chg_itt, & -!$OMP MKE_src_itt,Kddt_h_itt,dPEc_dKd,PE_chg, & -!$OMP dMKE_src_dK,TKE_left,use_Newt, & -!$OMP dKddt_h_Newt,Kddt_h_Newt,Kddt_h_next, & -!$OMP dKddt_h,Te,Se,Hsfc_used,dS_to_dPE_a, & -!$OMP dMKE_max,TKE_here) +!!OMP parallel do default(none) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, & +!!OMP CS,G,GV,fluxes,IdtdR0, & +!!OMP TKE_forced,debug,H_neglect,dSV_dT, & +!!OMP dSV_dS,I_dtmrho,C1_3,h_tt_min,vonKar, & +!!OMP max_itt,Kd_int) & +!!OMP private(i,j,k,h,u,v,T,S,Kd,mech_TKE_k,conv_PErel_k, & +!!OMP U_Star,absf,mech_TKE,conv_PErel,nstar_k, & +!!OMP h_sum,I_hs,h_bot,hb_hs,T0,S0,num_itts, & +!!OMP pres,dMass,dPres,dT_to_dPE,dS_to_dPE, & +!!OMP dT_to_dColHt,dS_to_dColHt,Kddt_h, & +!!OMP b_den_1,dT_to_dPE_a,dT_to_dColHt_a, & +!!OMP dS_to_dColHt_a,htot,uhtot,vhtot, & +!!OMP Idecay_len_TKE,exp_kh,nstar_FC,tot_TKE, & +!!OMP TKE_reduc,dTe_t2,dSe_t2,dTe,dSe,dt_h, & +!!OMP Convectively_stable,sfc_disconnect,b1, & +!!OMP c1,dT_km1_t2,dS_km1_t2,dTe_term, & +!!OMP dSe_term,MKE2_Hharm,vstar,h_tt, & +!!OMP Kd_guess0,Kddt_h_g0,dPEc_dKd_Kd0, & +!!OMP PE_chg_max,dPEa_dKd_g0,PE_chg_g0, & +!!OMP MKE_src,dPE_conv,Kddt_h_max,Kddt_h_min, & +!!OMP TKE_left_max,TKE_left_min,Kddt_h_guess, & +!!OMP TKE_left_itt,dPEa_dKd_itt,PE_chg_itt, & +!!OMP MKE_src_itt,Kddt_h_itt,dPEc_dKd,PE_chg, & +!!OMP dMKE_src_dK,TKE_left,use_Newt, & +!!OMP dKddt_h_Newt,Kddt_h_Newt,Kddt_h_next, & +!!OMP dKddt_h,Te,Se,Hsfc_used,dS_to_dPE_a, & +!!OMP dMKE_max,sfc_connected,TKE_here) do j=js,je ! Copy the thicknesses and other fields to 2-d arrays. do k=1,nz ; do i=is,ie diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 3be335cc75..7462d1c920 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -151,6 +151,10 @@ module MOM_vert_friction ! to calculate the viscous coupling between layers ! except near the bottom. Otherwise the arithmetic ! mean thickness is used except near the bottom. + real :: harm_BL_val ! A scale to determine when water is in the boundary + ! layers based solely on harmonic mean thicknesses + ! for the purpose of determining the extent to which + ! the thicknesses used in the viscosities are upwinded. logical :: direct_stress ! If true, the wind stress is distributed over the ! topmost Hmix_stress of fluid and KVML may be very small. logical :: dynamic_viscous_ML ! If true, use the results from a dynamic @@ -626,11 +630,15 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) real :: topfn ! A function which goes from 1 at the top to 0 much more ! than Htbl into the interior. real :: z2 ! The distance from the bottom, normalized by Hbbl, nondim. + real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected, in H. real :: H_to_m, m_to_H ! Unit conversion factors. real :: h_arith ! The arithmetic mean thickness, in m or kg m-2. + real :: I_valBL ! The inverse of a scaling factor determining when water is + ! still within the boundary layer, as determined by the sum + ! of the harmonic mean thicknesses. logical, dimension(SZIB_(G)) :: do_i, do_i_shelf logical :: do_any_shelf @@ -644,6 +652,7 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) h_neglect = GV%H_subroundoff H_to_m = GV%H_to_m ; m_to_H = GV%m_to_H I_Hbbl(:) = 1.0 / (CS%Hbbl * GV%m_to_H + h_neglect) + I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val if (CS%debug .or. (CS%id_hML_u > 0)) then allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; hML_u(:,:) = 0.0 @@ -662,11 +671,11 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) endif !$OMP parallel do default(none) shared(G,GV,CS,visc,Isq,ieq,nz,u,h,fluxes,hML_u, & -!$OMP h_neglect,dt,m_to_H) & +!$OMP h_neglect,dt,m_to_H,I_valBL) & !$OMP firstprivate(i_hbbl) & !$OMP private(do_i,kv_bbl,bbl_thick,z_i,h_harm,h_arith,hvel,z2, & !$OMP botfn,zh,Dmin,zcol,a,do_any_shelf,do_i_shelf, & -!$OMP a_shelf,Ztop_min,I_HTbl,hvel_shelf,topfn,h_ml) +!$OMP a_shelf,Ztop_min,I_HTbl,hvel_shelf,topfn,h_ml,z2_wt) do j=G%Jsc,G%Jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo @@ -710,9 +719,15 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) hvel(I,k) = h_arith if (u(I,j,k) * (h(i+1,j,k)-h(i,j,k)) > 0) then - z2 = max(zh(I), max(zcol(i),zcol(i+1)) + Dmin(I)) * I_Hbbl(I) - botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) - hvel(I,k) = (1.0-botfn)*h_arith + botfn*h_harm(I,k) + if (zh(I) * I_Hbbl(I) < CS%harm_BL_val) then + hvel(I,k) = h_harm(I,k) + else + z2_wt = 1.0 ; if (zh(I) * I_Hbbl(I) < 2.0*CS%harm_BL_val) & + z2_wt = max(0.0, min(1.0, zh(I) * I_Hbbl(I) * I_valBL - 1.0)) + z2 = z2_wt * (max(zh(I), max(zcol(i),zcol(i+1)) + Dmin(I)) * I_Hbbl(I)) + botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) + hvel(I,k) = (1.0-botfn)*h_arith + botfn*h_harm(I,k) + endif endif endif ; enddo ! i loop @@ -750,9 +765,15 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) hvel_shelf(I,k) = hvel(I,k) if (u(I,j,k) * (h(i+1,j,k)-h(i,j,k)) > 0) then - z2 = max(zh(I), Ztop_min(I) - min(zcol(i),zcol(i+1))) * I_HTbl(I) - topfn = 1.0 / (1.0 + 0.09*z2**6) - hvel_shelf(I,k) = min(hvel(I,k), (1.0-topfn)*h_arith + topfn*h_harm(I,k)) + if (zh(I) * I_HTbl(I) < CS%harm_BL_val) then + hvel_shelf(I,k) = min(hvel(I,k), h_harm(I,k)) + else + z2_wt = 1.0 ; if (zh(I) * I_HTbl(I) < 2.0*CS%harm_BL_val) & + z2_wt = max(0.0, min(1.0, zh(I) * I_HTbl(I) * I_valBL - 1.0)) + z2 = z2_wt * (max(zh(I), Ztop_min(I) - min(zcol(i),zcol(i+1))) * I_HTbl(I)) + topfn = 1.0 / (1.0 + 0.09*z2**6) + hvel_shelf(I,k) = min(hvel(I,k), (1.0-topfn)*h_arith + topfn*h_harm(I,k)) + endif endif endif ; enddo enddo @@ -766,8 +787,11 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) if (do_any_shelf) then do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i_shelf(I)) then - CS%a_u(I,j,K) = fluxes%frac_shelf_u(I,j) * a_shelf(I,k) + & + CS%a_u(I,j,K) = fluxes%frac_shelf_u(I,j) * a_shelf(I,K) + & (1.0-fluxes%frac_shelf_u(I,j)) * a(I,K) +! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH +! CS%a_u(I,j,K) = fluxes%frac_shelf_u(I,j) * max(a_shelf(I,K), a(I,K)) + & +! (1.0-fluxes%frac_shelf_u(I,j)) * a(I,K) elseif (do_i(I)) then CS%a_u(I,j,K) = a(I,K) endif ; enddo ; enddo @@ -788,11 +812,11 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) ! Now work on v-points. !$OMP parallel do default(none) shared(G,GV,CS,visc,is,ie,Jsq,Jeq,nz,v,h,fluxes,hML_v, & -!$OMP h_neglect,dt,m_to_H) & +!$OMP h_neglect,dt,m_to_H,I_valBL) & !$OMP firstprivate(i_hbbl) & !$OMP private(do_i,kv_bbl,bbl_thick,z_i,h_harm,h_arith,hvel,z2, & !$OMP botfn,zh,Dmin,zcol1,zcol2,a,do_any_shelf,do_i_shelf, & -!$OMP a_shelf,Ztop_min,I_HTbl,hvel_shelf,topfn,h_ml) +!$OMP a_shelf,Ztop_min,I_HTbl,hvel_shelf,topfn,h_ml,z2_wt) do J=Jsq,Jeq do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo @@ -837,10 +861,16 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) hvel(i,k) = h_arith if (v(i,J,k) * (h(i,j+1,k)-h(i,j,k)) > 0) then - z2 = max(zh(i), max(zcol1(i),zcol2(i)) + Dmin(i)) * I_Hbbl(i) - botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) - hvel(i,k) = (1.0-botfn)*h_arith + botfn*h_harm(i,k) - endif + if (zh(i) * I_Hbbl(i) < CS%harm_BL_val) then + hvel(i,k) = h_harm(i,k) + else + z2_wt = 1.0 ; if (zh(i) * I_Hbbl(i) < 2.0*CS%harm_BL_val) & + z2_wt = max(0.0, min(1.0, zh(i) * I_Hbbl(i) * I_valBL - 1.0)) + z2 = z2_wt * (max(zh(i), max(zcol1(i),zcol2(i)) + Dmin(i)) * I_Hbbl(i)) + botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) + hvel(i,k) = (1.0-botfn)*h_arith + botfn*h_harm(i,k) + endif + endif endif ; enddo ; enddo ! i & k loops endif @@ -876,10 +906,16 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) hvel_shelf(i,k) = hvel(i,k) if (v(i,j,k) * (h(i,j+1,k)-h(i,j,k)) > 0) then - z2 = max(zh(i), Ztop_min(i) - min(zcol1(i),zcol2(i))) * I_HTbl(i) - topfn = 1.0 / (1.0 + 0.09*z2**6) - hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith + topfn*h_harm(i,k)) - endif + if (zh(i) * I_HTbl(i) < CS%harm_BL_val) then + hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k)) + else + z2_wt = 1.0 ; if (zh(i) * I_HTbl(i) < 2.0*CS%harm_BL_val) & + z2_wt = max(0.0, min(1.0, zh(i) * I_HTbl(i) * I_valBL - 1.0)) + z2 = z2_wt * (max(zh(i), Ztop_min(i) - min(zcol1(i),zcol2(i))) * I_HTbl(i)) + topfn = 1.0 / (1.0 + 0.09*z2**6) + hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith + topfn*h_harm(i,k)) + endif + endif endif ; enddo enddo call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, & @@ -894,6 +930,9 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, GV, CS) do K=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then CS%a_v(i,J,K) = fluxes%frac_shelf_v(i,J) * a_shelf(i,k) + & (1.0-fluxes%frac_shelf_v(i,J)) * a(i,K) +! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH +! CS%a_v(i,J,K) = fluxes%frac_shelf_v(i,J) * max(a_shelf(i,K), a(i,K)) + & +! (1.0-fluxes%frac_shelf_v(i,J)) * a(i,K) elseif (do_i(i)) then CS%a_v(i,J,K) = a(i,K) endif ; enddo ; enddo @@ -1435,6 +1474,12 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, & call get_param(param_file, mod, "HARMONIC_VISC", CS%harmonic_visc, & "If true, use the harmonic mean thicknesses for \n"//& "calculating the vertical viscosity.", default=.false.) + call get_param(param_file, mod, "HARMONIC_BL_SCALE", CS%harm_BL_val, & + "A scale to determine when water is in the boundary \n"//& + "layers based solely on harmonic mean thicknesses for \n"//& + "the purpose of determining the extent to which the \n"//& + "thicknesses used in the viscosities are upwinded.", & + default=0.0, units="nondim") call get_param(param_file, mod, "DEBUG", CS%debug, default=.false.) if (GV%nkml < 1) &