From af232757aa0958805e4a77c96c790c848381f5f0 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 11 May 2016 13:00:38 -0400 Subject: [PATCH 1/3] Fix answer mismatch between regular and openmp runs - global_ALE/z was not reproducing ocean.stats when compiled with -openmp - This branch fixes issue # 290 , partially. - The remaining issue is that OMP directives in MOM_energetic_PBL.F90 still cause answer change and I cannot find anything wrong with them. In this update I just commented out the OMP directives in MOM_energetic_PBL.F90 and hopefully myself or Zhi can come back to it later to find the issue. --- src/ALE/MOM_ALE.F90 | 4 +- src/core/MOM_isopycnal_slopes.F90 | 8 +-- .../vertical/MOM_energetic_PBL.F90 | 60 +++++++++---------- 3 files changed, 36 insertions(+), 36 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 23c4159a95..ae0eca9a27 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -673,8 +673,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 From b5a1a139028845c620c6204e6f3cc29aacddea48 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 11 May 2016 18:36:23 -0400 Subject: [PATCH 2/3] +Added parameter to thicken viscous boundary layers Added a new run-time parameter, HARMONIC_BL_SCALE, that forces the model to treat any layers that are within this fraction of the bottom boundary layer thickess of the bottom as though they are at the bottom for the purpose of upwinding estimates of the thickness at velocity points. The default value, 0, reproduces the previous answers bitwise identically. --- .../vertical/MOM_vert_friction.F90 | 75 +++++++++++++++---- 1 file changed, 60 insertions(+), 15 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 3be335cc75..784951a8a3 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 @@ -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 @@ -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) & From 2f542864d1173f66011600d4a338d75eabaf3611 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 12 May 2016 09:54:51 -0400 Subject: [PATCH 3/3] Openmp fixes for MOM_vert_friction.F90 - Added missing variables to directives --- src/parameterizations/vertical/MOM_vert_friction.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 784951a8a3..7462d1c920 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -671,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 @@ -812,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