diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 760b195de4..f0335aaac4 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -283,23 +283,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, AD, G, CS) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke h_neglect = G%H_subroundoff +!$OMP parallel default(none) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,nz,RV,PV, & +!$OMP is,ie,js,je,Isq,Ieq,Jsq,Jeq,h_neglect) & +!$OMP private(relative_vorticity,absolute_vorticity,Ih,hArea_q,q, & +!$OMP abs_vort,Ih_q,fv1,fv2,fu1,fu2,max_fvq,max_fuq, & +!$OMP min_fvq,min_fuq,q2,a,b,c,d,ep_u,ep_v,Fe_m2,rat_lin, & +!$OMP min_Ihq,max_Ihq,rat_m1,AL_wt,Sad_wt,c1,c2,c3,slope, & +!$OMP uhc,uhm,uh_min,uh_max,vhc,vhm,vh_min,vh_max,KE,KEx, & +!$OMP KEy,temp1,temp2,Heff1,Heff2,Heff3,Heff4,VHeff, & +!$OMP QVHeff,max_fv,min_fv,UHeff,QUHeff,max_fu,min_fu) +!$OMP do do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+2 Area_h(i,j) = G%mask2dT(i,j) * G%areaT(i,j) enddo ; enddo +!$OMP do do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 Area_q(i,j) = (Area_h(i,j) + Area_h(i+1,j+1)) + & (Area_h(i+1,j) + Area_h(i,j+1)) enddo ; enddo -!$OMP parallel do default(none) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,nz,RV,PV, & -!$OMP is,ie,js,je,Isq,Ieq,Jsq,Jeq,h_neglect) & -!$OMP private(relative_vorticity,absolute_vorticity,Ih,hArea_q,q, & -!$OMP abs_vort,Ih_q,fv1,fv2,fu1,fu2,max_fvq,max_fuq, & -!$OMP min_fvq,min_fuq,q2,a,b,c,d,ep_u,ep_v,Fe_m2,rat_lin, & -!$OMP min_Ihq,max_Ihq,rat_m1,AL_wt,Sad_wt,c1,c2,c3,slope, & -!$OMP uhc,uhm,uh_min,uh_max,vhc,vhm,vh_min,vh_max,KE,KEx, & -!$OMP KEy,temp1,temp2,Heff1,Heff2,Heff3,Heff4,VHeff, & -!$OMP QVHeff,max_fv,min_fv,UHeff,QUHeff,max_fu,min_fu) +!$OMP do do k=1,nz ! Here the second order accurate layer potential vorticities, q, ! are calculated. hq is second order accurate in space. Relative @@ -727,6 +730,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, AD, G, CS) endif enddo ! k-loop. +!$OMP end parallel ! Here the various Coriolis-related derived quantities are offered ! for averaging. diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 13e01ab8e3..528159f86c 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -469,28 +469,34 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! and loading. This should really be based on bottom pressure anomalies, ! but that is not yet implemented, and the current form is correct for ! barotropic tides. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,1) = -1.0*G%bathyT(i,j) - enddo ; enddo - do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,1) = e(i,j,1) + h(i,j,k) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,h,G) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 ; e(i,j,1) = -1.0*G%bathyT(i,j) ; enddo + do k=1,nz ; do i=Isq,Ieq+1 + e(i,j,1) = e(i,j,1) + h(i,j,k) + enddo ; enddo + enddo call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp) endif ! Here layer interface heights, e, are calculated. +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,h,G,e_tidal,CS) if (CS%tides) then +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = -1.0*G%bathyT(i,j) - e_tidal(i,j) enddo ; enddo else +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = -1.0*G%bathyT(i,j) enddo ; enddo endif - do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP do + do j=Jsq,Jeq+1 ; do k=nz,1,-1 ; do i=Isq,Ieq+1 e(i,j,K) = e(i,j,K+1) + h(i,j,k) enddo ; enddo ; enddo +!$OMP end parallel if (use_EOS) then ! Calculate in-situ densities (rho_star). @@ -503,28 +509,34 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) if (nkmb>0) then tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp tv_tmp%eqn_of_state => tv%eqn_of_state - do k=1,nkmb ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - enddo ; enddo ; enddo + do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo - do j=Jsq,Jeq+1 +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref, & +!$OMP Rho_cv_BL,G) + do j=Jsq,Jeq+1 + do k=1,nkmb ; do i=Isq,Ieq+1 + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, & Rho_cv_BL(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state) + + do k=nkmb+1,nz ; do i=Isq,Ieq+1 + if (G%Rlay(k) < Rho_cv_BL(i,j)) then + tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) + else + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + endif + enddo ; enddo enddo - do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (G%Rlay(k) < Rho_cv_BL(i,j)) then - tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) - else - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - endif - enddo ; enddo ; enddo else tv_tmp%T => tv%T ; tv_tmp%S => tv%S tv_tmp%eqn_of_state => tv%eqn_of_state + do i=Isq,Ieq+1 ; p_ref(i) = 0.0 ; enddo endif ! This no longer includes any pressure dependency, since this routine ! will come down with a fatal error if there is any compressibility. +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,tv_tmp,p_ref,rho_star,tv,G_Rho0) do k=1,nz+1 ; do j=Jsq,Jeq+1 call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), & Isq,Ieq-Isq+2,tv%eqn_of_state) @@ -534,21 +546,28 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! Here the layer Montgomery potentials, M, are calculated. if (use_EOS) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1)) - if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 - enddo ; enddo - do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,CS,rho_star,e,use_p_atm, & +!$OMP p_atm,I_Rho0) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1)) + if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 + enddo + do k=2,nz ; do i=Isq,Ieq+1 + M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K) + enddo ; enddo + enddo else ! not use_EOS - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,1) = G%g_prime(1) * e(i,j,1) - if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 - enddo ; enddo - do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - M(i,j,k) = M(i,j,k-1) + G%g_prime(K) * e(i,j,K) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,M,G,e,use_p_atm,p_atm,I_Rho0) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + M(i,j,1) = G%g_prime(1) * e(i,j,1) + if (use_p_atm) M(i,j,1) = M(i,j,1) + p_atm(i,j) * I_Rho0 + enddo + do k=2,nz ; do i=Isq,Ieq+1 + M(i,j,k) = M(i,j,k-1) + G%g_prime(K) * e(i,j,K) + enddo ; enddo + enddo endif ! use_EOS if (present(pbce)) then @@ -559,6 +578,9 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! Calculate the pressure force. On a Cartesian grid, ! PFu = - dM/dx and PFv = - dM/dy. if (use_EOS) then +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,js,je,is,ie,nz,e,h_neglect, & +!$OMP rho_star,G,PFu,CS,PFv,M) & +!$OMP private(h_star,PFu_bc,PFv_bc) do k=1,nz do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 h_star(i,j) = (e(i,j,K) - e(i,j,K+1)) + h_neglect @@ -579,6 +601,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) enddo ; enddo enddo ! k-loop else ! .not. use_EOS +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,is,ie,js,je,nz,PFu,PFv,M,G) do k=1,nz do j=js,je ; do I=Isq,Ieq PFu(I,j,k) = -(M(i+1,j,k) - M(i,j,k)) * G%IdxCu(I,j) @@ -594,10 +617,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! eta is the sea surface height relative to a time-invariant geoid, for ! comparison with what is used for eta in btstep. See how e was calculated ! about 200 lines above. +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = e(i,j,1) + e_tidal(i,j) enddo ; enddo else +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = e(i,j,1) enddo ; enddo diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 386c3b4e56..31dd13e109 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -447,21 +447,27 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! astronomical sources and self-attraction and loading, in m. dM, & ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model, in m2 s-2. - pa, & ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the + pa ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the ! the interface atop a layer, in Pa. + real, dimension(SZI_(G)) :: & + Rho_cv_BL ! The coordinate potential density in the deepest variable + ! density near-surface layer, in kg m-3. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & dpa, & ! The change in pressure anomaly between the top and bottom ! of a layer, in Pa. - Rho_cv_BL,& ! The coordinate potential density in the deepest variable - ! density near-surface layer, in kg m-3. - intz_dpa ! The vertical integral in depth of the pressure anomaly less + intz_dpa, & ! The vertical integral in depth of the pressure anomaly less ! the pressure anomaly at the top of the layer, in m Pa. + pa_h_intz_dpa ! pa*h+intz_dpa + real, dimension(SZIB_(G),SZJ_(G)) :: & - intx_pa, & ! The zonal integral of the pressure anomaly along the interface + intx_pa ! The zonal integral of the pressure anomaly along the interface ! atop a layer, divided by the grid spacing, in Pa. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & intx_dpa ! The change in intx_pa through a layer, in Pa. real, dimension(SZI_(G),SZJB_(G)) :: & - inty_pa, & ! The meridional integral of the pressure anomaly along the + inty_pa ! The meridional integral of the pressure anomaly along the ! interface atop a layer, divided by the grid spacing, in Pa. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & inty_dpa ! The change in inty_pa through a layer, in Pa. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: & T_tmp, & ! Temporary array of temperatures where layers that are lighter @@ -488,6 +494,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, real, parameter :: C1_6 = 1.0/6.0 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k + integer :: PRScheme is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke nkmb=G%nk_rho_varies @@ -501,6 +508,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce: Module must be initialized before it is used.") + PRScheme = pressureReconstructionScheme(ALE_CSp) h_neglect = G%H_subroundoff I_Rho0 = 1.0/G%Rho0 G_Rho0 = G%g_Earth/G%Rho0 @@ -511,28 +519,36 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! and loading. This should really be based on bottom pressure anomalies, ! but that is not yet implemented, and the current form is correct for ! barotropic tides. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,1) = -1.0*G%bathyT(i,j) - enddo ; enddo - do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,1) = e(i,j,1) + G%H_to_m*h(i,j,k) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + e(i,j,1) = -1.0*G%bathyT(i,j) + enddo + do k=1,nz ; do i=Isq,Ieq+1 + e(i,j,1) = e(i,j,1) + G%H_to_m*h(i,j,k) + enddo ; enddo + enddo call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp) endif ! Here layer interface heights, e, are calculated. +!$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h,CS,e_tidal) if (CS%tides) then +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = -1.0*G%bathyT(i,j) - e_tidal(i,j) enddo ; enddo else +!$OM do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = -1.0*G%bathyT(i,j) enddo ; enddo endif - do k=nz,1,-1 ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 +!$OMP do + do j=Jsq,Jeq+1; do k=nz,1,-1 ; do i=Isq,Ieq+1 e(i,j,K) = e(i,j,K+1) + G%H_to_m*h(i,j,k) enddo ; enddo ; enddo +!$OMP end parallel if (use_EOS) then @@ -544,30 +560,38 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, if (nkmb>0) then tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp tv_tmp%eqn_of_state => tv%eqn_of_state - do k=1,nkmb ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - enddo ; enddo ; enddo + do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nkmb,nz,G,tv_tmp,tv,p_ref) & +!$OMP private(Rho_cv_BL) do j=Jsq,Jeq+1 + do k=1,nkmb ; do i=Isq,Ieq+1 + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, & - Rho_cv_BL(:,j), Isq, Ieq-Isq+2, tv%eqn_of_state) + Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state) + + do k=nkmb+1,nz ; do i=Isq,Ieq+1 + if (G%Rlay(k) < Rho_cv_BL(i)) then + tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) + else + tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) + endif + enddo ; enddo enddo - do k=nkmb+1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (G%Rlay(k) < Rho_cv_BL(i,j)) then - tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) - else - tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) - endif - enddo ; enddo ; enddo else tv_tmp%T => tv%T ; tv_tmp%S => tv%S tv_tmp%eqn_of_state => tv%eqn_of_state endif endif +!$OMP parallel default(none) shared(Jsq,Jeq,Isq,Ieq,tv_tmp,p_atm,rho_in_situ,tv, & +!$OMP p0,dM,CS,G_Rho0,e,use_p_atm,use_EOS,G,pa, & +!$OMP rho_ref,js,je,is,ie,intx_pa,inty_pa) if (CS%GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. if (use_EOS) then +!$OMP do do j=Jsq,Jeq+1 if (use_p_atm) then call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), & @@ -581,6 +605,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, enddo enddo else +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * G%Rlay(1)) * e(i,j,1) enddo ; enddo @@ -591,20 +616,25 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! integrals, assuming that the surface pressure anomaly varies linearly ! in x and y. if (use_p_atm) then +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) + p_atm(i,j) enddo ; enddo else +!$OMP do do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 pa(i,j) = (rho_ref*G%g_Earth)*e(i,j,1) enddo ; enddo endif +!$OMP do do j=js,je ; do I=Isq,Ieq intx_pa(I,j) = 0.5*(pa(i,j) + pa(i+1,j)) enddo ; enddo +!$OMP do do J=Jsq,Jeq ; do i=is,ie inty_pa(i,J) = 0.5*(pa(i,j) + pa(i,j+1)) enddo ; enddo +!$OMP end parallel ! Have checked that rho_0 drops out and that the 1-layer case is right. RWH. @@ -616,13 +646,19 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, use_ALE = .false. if (associated(ALE_CSp)) use_ALE = usePressureReconstruction(ALE_CSp) if ( use_ALE ) then - if ( pressureReconstructionScheme(ALE_CSp) == PRESSURE_RECONSTRUCTION_PLM ) then + if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then call pressure_gradient_plm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h ); - elseif ( pressureReconstructionScheme(ALE_CSp) == PRESSURE_RECONSTRUCTION_PPM ) then + elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then call pressure_gradient_ppm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h ); endif endif +!$OMP parallel default(none) shared(nz,Jsq,Jeq,Isq,Ieq,G,h,use_EOS,use_ALE,PRScheme,& +!$OMP T_t,T_b,S_t,S_b,e,rho_ref,CS,tv,js,je,is,ie, & +!$OMP dpa,intz_dpa,intx_dpa,inty_dpa,tv_tmp,I_Rho0,pa, & +!$OMP pa_h_intz_dpa,PFu,PFv,intx_pa,h_neglect,inty_pa) & +!$OMP private(dz) +!$OMP do do k=1,nz ! Calculate 4 integrals through the layer that are required in the ! subsequent calculation. @@ -636,67 +672,75 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! is used, whereby densities within each layer are constant no matter ! where the layers are located. if ( use_ALE ) then - if ( pressureReconstructionScheme(ALE_CSp) == PRESSURE_RECONSTRUCTION_PLM ) then + if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then call int_density_dz_generic_plm ( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), & S_b(:,:,k), & e(:,:,K), e(:,:,K+1), rho_ref, & CS%Rho0, G%g_Earth, & - G, tv%eqn_of_state, dpa, intz_dpa, & - intx_dpa, inty_dpa, & + G, tv%eqn_of_state, dpa(:,:,k), intz_dpa(:,:,k), & + intx_dpa(:,:,k), inty_dpa(:,:,k), & useMassWghtInterp = CS%useMassWghtInterp) - elseif ( pressureReconstructionScheme(ALE_CSp) == PRESSURE_RECONSTRUCTION_PPM ) then + elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then call int_density_dz_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), & e(:,:,K), e(:,:,K+1), rho_ref, & CS%Rho0, G%G_Earth, & - G, tv%eqn_of_state, dpa, intz_dpa, & - intx_dpa, inty_dpa) + G, tv%eqn_of_state, dpa(:,:,k), intz_dpa(:,:,k), & + intx_dpa(:,:,k), inty_dpa(:,:,k)) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, G%g_Earth, G, tv%eqn_of_state, & - dpa, intz_dpa, intx_dpa, inty_dpa) + dpa(:,:,k), intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k)) endif else do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dpa(i,j) = (G%Rlay(k) - rho_ref)*dz(i,j) - intz_dpa(i,j) = 0.5*(G%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k) + dpa(i,j,k) = (G%Rlay(k) - rho_ref)*dz(i,j) + intz_dpa(i,j,k) = 0.5*(G%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k) enddo ; enddo do j=js,je ; do I=Isq,Ieq - intx_dpa(I,j) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j)) + intx_dpa(I,j,k) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j)) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - inty_dpa(i,J) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1)) + inty_dpa(i,J,k) = 0.5*(G%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1)) enddo ; enddo endif + enddo - ! Compute pressure gradient in x direction - do j=js,je ; do I=Isq,Ieq - PFu(I,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - & - (pa(i+1,j)*h(i+1,j,k) + intz_dpa(i+1,j))) + & +!$OMP do + do j=jsq,Jeq+1; do k=1,nz; do i=Isq,Ieq+1 + pa_h_intz_dpa(i,j,k) = pa(i,j)*h(i,j,k) + intz_dpa(i,j,k) + pa(i,j) = pa(i,j)+dpa(i,j,k) + enddo; enddo; enddo + + ! Compute pressure gradient in x direction +!$OMP do + do j =js,je + do k = 1,nz; do I=Isq,Ieq + PFu(I,j,k) = ((pa_h_intz_dpa(i,j,k) - pa_h_intz_dpa(i+1,j,k)) + & ((h(i+1,j,k) - h(i,j,k)) * intx_pa(I,j) - & - (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa(I,j))) * & + (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa(I,j,k))) * & ((2.0*I_Rho0*G%IdxCu(I,j)) / & ((h(i,j,k) + h(i+1,j,k)) + h_neglect)) - intx_pa(I,j) = intx_pa(I,j) + intx_dpa(I,j) + intx_pa(I,j) = intx_pa(I,j) + intx_dpa(I,j,k) enddo ; enddo + enddo ! Compute pressure gradient in y direction - do J=Jsq,Jeq ; do i=is,ie - PFv(i,J,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - & - (pa(i,j+1)*h(i,j+1,k) + intz_dpa(i,j+1))) + & +!$OMP do + do J=Jsq,Jeq + do k=1,nz ; do i=is,ie + PFv(i,J,k) = ((pa_h_intz_dpa(i,j,k) - pa_h_intz_dpa(i,j+1,k)) + & ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,J) - & - (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa(i,J))) * & + (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa(i,J,k))) * & ((2.0*I_Rho0*G%IdyCv(i,J)) / & ((h(i,j,k) + h(i,j+1,k)) + h_neglect)) - inty_pa(i,J) = inty_pa(i,J) + inty_dpa(i,J) - enddo ; enddo - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j) = pa(i,j) + dpa(i,j) + inty_pa(i,J) = inty_pa(i,J) + inty_dpa(i,J,k) enddo ; enddo enddo - +!$OMP end parallel if (CS%GFS_scale < 1.0) then +!$OMP parallel do default(none) shared(is,ie,js,je,nz,Isq,Ieq,Jsq,Jeq,PFu,PFv,dM,G) do k=1,nz do j=js,je ; do I=Isq,Ieq PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j) @@ -716,10 +760,12 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! eta is the sea surface height relative to a time-invariant geoid, for ! comparison with what is used for eta in btstep. See how e was calculated ! about 200 lines above. +!$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = e(i,j,1) + e_tidal(i,j) enddo ; enddo else +!$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = e(i,j,1) enddo ; enddo diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index d18003daaa..e09759b718 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -690,6 +690,10 @@ subroutine zonal_face_thickness(u, h, hL, hR, h_u, dt, G, LB, vol_CFL, visc_rem_ integer :: i, j, k, ish, ieh, jsh, jeh, nz ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke +!$OMP parallel default(none) shared(ish,ieh,jsh,jeh,nz,u,vol_CFL,dt,G, & +!$OMP hL,hR,h,h_u,visc_rem_u) & +!$OMP private(CFL,curv_3,h_marg) +!$OMP do do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh if (u(I,j,k) > 0.0) then if (vol_CFL) then ; CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j)) @@ -714,8 +718,14 @@ subroutine zonal_face_thickness(u, h, hL, hR, h_u, dt, G, LB, vol_CFL, visc_rem_ endif ! This could perhaps be replaced by h_avg. h_u(I,j,k) = h_marg - if (present(visc_rem_u)) h_u(I,j,k) = h_u(I,j,k) * visc_rem_u(I,j,k) - enddo ; enddo ; enddo + enddo; enddo ; enddo + if (present(visc_rem_u)) then +!$OMP do + do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh + h_u(I,j,k) = h_u(I,j,k) * visc_rem_u(I,j,k) + enddo ; enddo ; enddo + endif +!$OMP end parallel end subroutine zonal_face_thickness @@ -1426,6 +1436,10 @@ subroutine merid_face_thickness(v, h, hL, hR, h_v, dt, G, LB, vol_CFL, visc_rem_ integer :: i, j, k, ish, ieh, jsh, jeh, nz ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke +!$OMP parallel default(none) shared(ish,ieh,jsh,jeh,nz,v,vol_CFL,dt,G, & +!$OMP hL,hR,h,h_v,visc_rem_v) & +!$OMP private(CFL,curv_3,h_marg) +!$OMP do do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh if (v(i,J,k) > 0.0) then if (vol_CFL) then ; CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) @@ -1451,8 +1465,15 @@ subroutine merid_face_thickness(v, h, hL, hR, h_v, dt, G, LB, vol_CFL, visc_rem_ endif ! This could perhaps be replaced by h_avg. h_v(i,J,k) = h_marg - if (present(visc_rem_v)) h_v(i,J,k) = h_v(i,J,k) * visc_rem_v(i,J,k) enddo ; enddo ; enddo + + if (present(visc_rem_v)) then +!$OMP do + do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh + h_v(i,J,k) = h_v(i,J,k) * visc_rem_v(i,J,k) + enddo ; enddo ; enddo + endif +!$OMP end parallel end subroutine merid_face_thickness diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 3819e0cb04..1989faed3d 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -1069,15 +1069,22 @@ subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, CS) call continuity(u, v, h, h_temp, uh_temp, vh_temp, dt, G, & CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) - - do j=js,je ; do I=is-1,ie ; CS%uhbt_in(I,j) = uh_temp(I,j,1) ; enddo ; enddo - do k=2,nz ; do j=js,je ; do I=is-1,ie - CS%uhbt_in(I,j) = CS%uhbt_in(I,j) + uh_temp(I,j,k) - enddo ; enddo ; enddo - do J=js-1,je ; do i=is,ie ; CS%vhbt_in(i,J) = vh_temp(i,J,1) ; enddo ; enddo - do k=2,nz ; do J=js-1,je ; do i=is,ie - CS%vhbt_in(i,J) = CS%vhbt_in(i,J) + vh_temp(i,J,k) - enddo ; enddo ; enddo +!$OMP parallel default(none) shared(is,ie,js,je,nz,CS,uh_temp,vh_temp) +!$OMP do + do j=js,je + do I=is-1,ie ; CS%uhbt_in(I,j) = uh_temp(I,j,1) ; enddo + do k=2,nz ; do I=is-1,ie + CS%uhbt_in(I,j) = CS%uhbt_in(I,j) + uh_temp(I,j,k) + enddo ; enddo + enddo +!$OMP do + do J=js-1,je + do i=is,ie ; CS%vhbt_in(i,J) = vh_temp(i,J,1) ; enddo + do k=2,nz ; do i=is,ie + CS%vhbt_in(i,J) = CS%vhbt_in(i,J) + vh_temp(i,J,k) + enddo ; enddo + enddo +!$OMP end parallel CS%readjust_velocity = .true. endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index e47c9816b4..ad598b2fc4 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -370,7 +370,13 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90") - up(:,:,:) = 0.0 ; vp(:,:,:) = 0.0 ; hp(:,:,:) = h(:,:,:) + +!$OMP parallel do default(none) shared(nz,G,up,vp,hp,h) + do k = 1, nz + do j=G%jsd,G%jed ; do i=G%isdB,G%iedB ; up(i,j,k) = 0.0 ; enddo ; enddo + do j=G%jsdB,G%jedB ; do i=G%isd,G%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo + do j=G%jsd,G%jed ; do i=G%isd,G%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo + enddo ! Update CFL truncation value as function of time call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp) @@ -459,6 +465,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & else Pa_to_eta = 1.0 / G%H_to_Pa endif +!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta_PF_start,CS,Pa_to_eta,p_surf_begin,p_surf_end) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta_PF_start(i,j) = CS%eta_PF(i,j) - Pa_to_eta * & (p_surf_begin(i,j) - p_surf_end(i,j)) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 2d5b27868b..0a3763cfb0 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -388,6 +388,9 @@ subroutine extractFluxes2d(G, fluxes, optics, nsw, dt, ! heat flux associated with net mass fluxes into the ocean. integer :: j +!$OMP parallel do default(none) shared(G,fluxes, optics, nsw,dt,DepthBeforeScalingFluxes, & +!$OMP useRiverHeatContent, useCalvingHeatContent, & +!$OMP h,T,Net_H,Net_heat,Net_salt,Pen_SW_bnd,tv) do j=G%jsc, G%jec call extractFluxes1d(G, fluxes, optics, nsw, j, dt, & DepthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, & diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index e31484c5c6..1cb060be56 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -91,8 +91,8 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) ! (in,opt) halo_size - The width of halo points on which to calculate eta. real :: p(SZI_(G),SZJ_(G),SZK_(G)+1) - real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height - ! across a layer, in m2 s-2. + real :: dz_geo(SZI_(G),SZJ_(G),SZK_(G)) ! The change in geopotential height + ! across a layer, in m2 s-2. real :: dilate(SZI_(G)), htot(SZI_(G)) real :: I_gEarth integer i, j, k, isv, iev, jsv, jev, nz, halo @@ -107,15 +107,21 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) I_gEarth = 1.0 / G_Earth +!$OMP parallel default(none) shared(isv,iev,jsv,jev,nz,eta,G,h,eta_bt,tv,p, & +!$OMP G_Earth,dz_geo,halo,I_gEarth) & +!$OMP private(dilate,htot) +!$OMP do do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -G%bathyT(i,j) ; enddo ; enddo if (G%Boussinesq) then - do k=nz,1,-1 ; do j=jsv,jev ; do i=isv,iev +!$OMP do + do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + h(i,j,k) enddo ; enddo ; enddo if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! that is used for the dynamics. +!$OMP do do j=jsv,jev do i=isv,iev dilate(i) = (eta_bt(i,j)+G%bathyT(i,j)) / (eta(i,j,1)+G%bathyT(i,j)) @@ -128,26 +134,34 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) else if (associated(tv%eqn_of_state)) then ! ### THIS SHOULD BE P_SURF, IF AVAILABLE. - do j=jsv,jev ; do i=isv,iev ; p(i,j,1) = 0.0 ; enddo ; enddo - do k=1,nz ; do j=jsv,jev ; do i=isv,iev - p(i,j,K+1) = p(i,j,K) + G_Earth*G%H_to_kg_m2*h(i,j,k) - enddo ; enddo ; enddo - - do k=nz,1,-1 +!$OMP do + do j=jsv,jev + do i=isv,iev ; p(i,j,1) = 0.0 ; enddo + do k=1,nz ; do i=isv,iev + p(i,j,K+1) = p(i,j,K) + G_Earth*G%H_to_kg_m2*h(i,j,k) + enddo ; enddo + enddo +!$OMP do + do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & - 0.0, G, tv%eqn_of_state, dz_geo, halo_size=halo) - do j=jsv,jev ; do i=isv,iev - eta(i,j,K) = eta(i,j,K+1) + I_gEarth * dz_geo(i,j) + 0.0, G, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo) + enddo +!$OMP do + do j=jsv,jev + do k=nz,1,-1 ; do i=isv,iev + eta(i,j,K) = eta(i,j,K+1) + I_gEarth * dz_geo(i,j,k) enddo ; enddo enddo else - do k=nz,1,-1 ; do j=jsv,jev ; do i=isv,iev +!$OMP do + do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + G%H_to_kg_m2*h(i,j,k)/G%Rlay(k) enddo ; enddo ; enddo endif if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! from the time-averaged barotropic solution. +!$OMP do do j=jsv,jev do i=isv,iev ; htot(i) = G%H_subroundoff ; enddo do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo @@ -158,6 +172,7 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) enddo endif endif +!$OMP end parallel end subroutine find_eta_3d @@ -184,9 +199,9 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size) ! mass per unit aread (non-Boussinesq), in m or kg m-2. ! (in,opt) halo_size - The width of halo points on which to calculate eta. - real, dimension(SZI_(G),SZJ_(G)) :: & - p_top, & ! The pressure at the interface above a layer, in Pa. - p_bot, & ! The pressure at the interface below a layer, in Pa. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: & + p ! The pressure in Pa. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & dz_geo ! The change in geopotential height across a layer, in m2 s-2. real :: htot(SZI_(G)) ! The sum of all layers' thicknesses, in kg m-2 or m. real :: I_gEarth @@ -198,42 +213,53 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size) I_gEarth = 1.0 / G_Earth +!$OMP parallel default(none) shared(is,ie,js,je,nz,eta,G,eta_bt,h,tv,p, & +!$OMP G_Earth,dz_geo,halo,I_gEarth) & +!$OMP private(htot) +!$OMP do do j=js,je ; do i=is,ie ; eta(i,j) = -G%bathyT(i,j) ; enddo ; enddo if (G%Boussinesq) then if (present(eta_bt)) then +!$OMP do do j=js,je ; do i=is,ie eta(i,j) = eta_bt(i,j) enddo ; enddo else - do k=1,nz ; do j=js,je ; do i=is,ie +!$OMP do + do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + h(i,j,k) enddo ; enddo ; enddo endif else if (associated(tv%eqn_of_state)) then - do j=js,je ; do i=is,ie - p_top(i,j) = 0.0 ; p_bot(i,j) = 0.0 - enddo ; enddo - do k=1,nz - do j=js,je ; do i=is,ie - p_top(i,j) = p_bot(i,j) - p_bot(i,j) = p_top(i,j) + G_Earth*G%H_to_kg_m2*h(i,j,k) - enddo ; enddo - call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, 0.0, & - G, tv%eqn_of_state, dz_geo, halo_size=halo) - do j=js,je ; do i=is,ie - eta(i,j) = eta(i,j) + I_gEarth * dz_geo(i,j) +!$OMP do + do j=js,je + do i=is,ie ; p(i,j,1) = 0.0 ; enddo + + do k=1,nz ; do i=is,ie + p(i,j,k+1) = p(i,j,k) + G_Earth*G%H_to_kg_m2*h(i,j,k) enddo ; enddo enddo +!$OMP do + do k = 1, nz + call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, & + G, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo) + enddo +!$OMP do + do j=js,je ; do k=1,nz ; do i=is,ie + eta(i,j) = eta(i,j) + I_gEarth * dz_geo(i,j,k) + enddo ; enddo ; enddo else - do k=1,nz ; do j=js,je ; do i=is,ie +!$OMP do + do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + G%H_to_kg_m2*h(i,j,k)/G%Rlay(k) enddo ; enddo ; enddo endif if (present(eta_bt)) then ! Dilate the water column to agree with the the time-averaged column ! mass from the barotropic solution. +!$OMP do do j=js,je do i=is,ie ; htot(i) = G%H_subroundoff ; enddo do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo @@ -244,7 +270,8 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size) enddo endif endif - +!$OMP end parallel + end subroutine find_eta_2d end module MOM_interface_heights diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index fd01dc095e..7a94aabc19 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -494,7 +494,7 @@ subroutine set_grid_derived_metrics(G, param_file) G%IareaBu(I,J) = G%IdxBu(I,J) * G%IdyBu(I,J) enddo ; enddo -68 FORMAT ("WARNING: PE ",I4," ",a3,"(",I4,",",I4,") = ",ES10.4, & +68 FORMAT ("WARNING: PE ",I4," ",a4,"(",I4,",",I4,") = ",ES10.4, & " is being changed to ",ES10.4,".") call callTree_leave("set_grid_derived_metrics()") diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 80c3a3c0ac..3866fd4d36 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -93,39 +93,59 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping + Rho0 = G%H_to_kg_m2 * G%m_to_H + mass_neglect = G%H_to_kg_m2 * G%H_subroundoff if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_MEKE: Module must be initialized before it is used.") if (.not.associated(MEKE)) call MOM_error(FATAL, & "MOM_MEKE: MEKE must be initialized before it is used.") + ! Only integerate the MEKE equations if MEKE is required. if (associated(MEKE%MEKE)) then sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping - if (CS%id_src > 0) then ; do j=js,je ; do i=is,ie - src(i,j) = CS%MEKE_BGsrc - enddo ; enddo ; endif - Rho0 = G%H_to_kg_m2 * G%m_to_H mass_neglect = G%H_to_kg_m2 * G%H_subroundoff + cdrag2 = CS%MEKE_Cd_scale * CS%cdrag**2 + ! 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) sdt_damp = 0.5*sdt - do j=js-1,je+1 ; do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo ; enddo - do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 - mass(i,j) = mass(i,j) + G%H_to_kg_m2 * h(i,j,k) - enddo ; enddo ; enddo - do j=js-1,je+1 ; do i=is-1,ie+1 - I_mass(i,j) = 0.0 - if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) - if (mass(i,j) < 0.0) mass(i,j) = 0.0 ! Should be unneccesary? - enddo ; enddo +!$OMP parallel default(none) shared(MEKE,CS,is,ie,js,je,nz,src,mass,G,h,I_mass, & +!$OMP sdt,drag_vel_u,visc,drag_vel_v,drag_rate_visc, & +!$OMP drag_rate,Rho0,MEKE_decay,sdt_damp,cdrag2) & +!$OMP private(ldamping) + if (CS%id_src > 0) then +!$OMP do + do j=js,je ; do i=is,ie + src(i,j) = CS%MEKE_BGsrc + enddo ; enddo + endif + +!$OMP do + do j=js-1,je+1 + do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo + do k=1,nz ; do i=is-1,ie+1 + mass(i,j) = mass(i,j) + G%H_to_kg_m2 * h(i,j,k) + enddo ; enddo + do i=is-1,ie+1 + I_mass(i,j) = 0.0 + if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) + if (mass(i,j) < 0.0) mass(i,j) = 0.0 ! Should be unneccesary? + enddo + enddo if (associated(MEKE%mom_src)) then +!$OMP do do j=js,je ; do i=is,ie MEKE%MEKE(i,j) = MEKE%MEKE(i,j) - I_mass(i,j) * & (sdt*CS%MEKE_FrCoeff)*MEKE%mom_src(i,j) enddo ; enddo if (CS%id_src > 0) then +!$OMP do 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) enddo ; enddo @@ -133,17 +153,20 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) endif if (associated(MEKE%GM_src)) then +!$OMP do do j=js,je ; do i=is,ie MEKE%MEKE(i,j) = MEKE%MEKE(i,j) - I_mass(i,j) * & (sdt*CS%MEKE_GMcoeff)*MEKE%GM_src(i,j) enddo ; enddo if (CS%id_src > 0) then +!$OMP do 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) enddo ; enddo endif endif +!$OMP do do j=js,je ; do i=is,ie MEKE%MEKE(i,j) = max(0.0,MEKE%MEKE(i,j) + (sdt*CS%MEKE_BGsrc)*G%mask2dT(i,j)) enddo ; enddo @@ -153,18 +176,20 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0) sdt_damp = 0.5*sdt if (CS%visc_drag) then +!$OMP do do j=js,je ; do I=is-1,ie drag_vel_u(I,j) = 0.0 if ((G%mask2dCu(I,j) > 0.0) .and. (visc%bbl_thick_u(I,j) > 0.0)) & drag_vel_u(I,j) = visc%kv_bbl_u(I,j) / visc%bbl_thick_u(I,j) enddo ; enddo +!$OMP do do J=js-1,je ; do i=is,ie drag_vel_v(i,J) = 0.0 if ((G%mask2dCu(i,J) > 0.0) .and. (visc%bbl_thick_v(i,J) > 0.0)) & drag_vel_v(i,J) = visc%kv_bbl_v(i,J) / visc%bbl_thick_v(i,J) enddo ; enddo - cdrag2 = CS%MEKE_Cd_scale * CS%cdrag**2 +!$OMP do do j=js,je ; do i=is,ie drag_rate_visc(i,j) = (0.25*G%IareaT(i,j) * & ((G%areaCu(I-1,j)*drag_vel_u(I-1,j) + & @@ -175,6 +200,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) sqrt(drag_rate_visc(i,j)**2 + 2.0*cdrag2 * max(0.0,MEKE%MEKE(i,j))) enddo ; enddo elseif (CS%MEKE_Cd_scale >= 0.0) then +!$OMP do do j=js,je ; do i=is,ie drag_rate(i,j) = (Rho0 * I_mass(i,j)) * & (CS%cdrag*(2.0*sqrt(max(0.0,MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2))) @@ -182,6 +208,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) endif if (CS%MEKE_damping + CS%MEKE_Cd_scale > 0.0) then +!$OMP do do j=js,je ; do i=is,ie ldamping = CS%MEKE_damping if (CS%MEKE_Cd_scale > 0.0) & @@ -191,7 +218,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) enddo ; enddo endif - +!$OMP end parallel if (CS%MEKE_KH >= 0.0) then call cpu_clock_begin(id_clock_pass) call pass_var(MEKE%MEKE, G%Domain) @@ -199,6 +226,12 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) ! ### More elaborate prescriptions for Kh could be used here. Kh_here = CS%MEKE_Kh +!$OMP parallel default(none) shared(is,ie,js,je,MEKE,CS,sdt,G,Kh_u,MEKE_uflux, & +!$OMP mass,mass_neglect,Kh_v,MEKE_vflux,I_mass, & +!$OMP sdt_damp,drag_rate,Rho0,drag_rate_visc, & +!$OMP cdrag2,MEKE_decay) & +!$OMP private(Kh_here,Inv_Kh_max,ldamping) +!$OMP do do j=js,je ; do I=is-1,ie ! Limit Kh to avoid CFL violations. if (associated(MEKE%Kh)) & @@ -212,6 +245,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * & (MEKE%MEKE(i,j) - MEKE%MEKE(i+1,j)) enddo ; enddo +!$OMP do do J=js-1,je ; do i=is,ie if (associated(MEKE%Kh)) & Kh_here = CS%MEKE_Kh + CS%KhMEKE_Fac*0.5*(MEKE%Kh(i,j)+MEKE%Kh(i,j+1)) @@ -224,6 +258,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * & (MEKE%MEKE(i,j) - MEKE%MEKE(i,j+1)) enddo ; enddo +!$OMP do do j=js,je ; do i=is,ie 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)) + & @@ -233,17 +268,19 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) if ((CS%MEKE_damping + CS%MEKE_Cd_scale > 0.0) .and. (sdt>sdt_damp)) then ! Recalculate the drag rate, since MEKE has changed. if (CS%visc_drag) then +!$OMP do do j=js,je ; do i=is,ie drag_rate(i,j) = (Rho0 * I_mass(i,j)) * & sqrt(drag_rate_visc(i,j)**2 + cdrag2 * 2.0*max(0.0,MEKE%MEKE(i,j))) enddo ; enddo elseif (CS%MEKE_Cd_scale >= 0.0) then +!$OMP do do j=js,je ; do i=is,ie drag_rate(i,j) = (Rho0 * I_mass(i,j)) * & (CS%cdrag*2.0*(sqrt(max(0.0,MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2))) enddo ; enddo endif - +!$OMP do do j=js,je ; do i=is,ie ldamping = CS%MEKE_damping if (CS%MEKE_Cd_scale > 0.0) & @@ -253,6 +290,7 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) MEKE_decay(i,j) = 0.5 * G%mask2dT(i,j) * (MEKE_decay(i,j) + ldamping) enddo ; enddo endif +!$OMP end parallel endif call cpu_clock_begin(id_clock_pass) @@ -262,11 +300,13 @@ subroutine step_forward_MEKE(MEKE, h, visc, dt, G, CS) ! Calculate diffusivity for main model to use if (CS%MEKE_KhCoeff>0.) then if (CS%Rd_as_max_scale) then +!$OMP parallel do default(none) shared(is,ie,js,je,MEKE,CS,G) do j=js-1,je+1 ; do i=is-1,ie+1 MEKE%Kh(i,j) = (CS%MEKE_KhCoeff*sqrt(2.*max(0.,MEKE%MEKE(i,j))*G%areaT(i,j))) * & min(MEKE%Rd_dx_h(i,j), 1.0) enddo ; enddo else +!$OMP parallel do default(none) shared(is,ie,js,je,MEKE,CS,G) do j=js-1,je+1 ; do i=is-1,ie+1 MEKE%Kh(i,j) = CS%MEKE_KhCoeff*sqrt(2.*max(0.,MEKE%MEKE(i,j))*G%areaT(i,j)) enddo ; enddo diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 8149203ffc..63d7d6893c 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -551,10 +551,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, CS, OBC) enddo ! end of k loop if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then - do j=js,je ; do i=is,ie ; MEKE%mom_src(i,j) = FrictWork(i,j,1) ; enddo ; enddo - do k=2,nz ; do j=js,je ; do i=is,ie - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) - enddo ; enddo ; enddo +!$OMP parallel do default(none) shared(is,ie,js,je,nz,MEKE,FrictWork) + do j=js,je + do i=is,ie ; MEKE%mom_src(i,j) = FrictWork(i,j,1) ; enddo + do k=2,nz ; do i=is,ie + MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) + enddo ; enddo + enddo endif ; endif ! Offer fields for averaging. diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index eb9ab5fe6e..493a0cf21e 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -145,7 +145,11 @@ subroutine calc_resoln_function(h, tv, G, CS) ! Do this calculation on the extent used in MOM_hor_visc.F90, and ! MOM_tracer.F90 so that no halo update is needed. mod_power_2 = mod(CS%Res_fn_power, 2) + +!$OMP parallel default(none) shared(is,ie,js,je,Ieq,Jeq,CS,mod_power_2) & +!$OMP private(dx_term,cg1_q,power_2) if (CS%Res_fn_power >= 100) then +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 dx_term = CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j) if ((CS%Res_coef * CS%cg1(i,j))**2 > dx_term) then @@ -154,7 +158,7 @@ subroutine calc_resoln_function(h, tv, G, CS) CS%Res_fn_h(i,j) = 1.0 endif enddo ; enddo - +!$OMP do do J=js-1,Jeq ; do I=is-1,Ieq cg1_q = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + & (CS%cg1(i+1,j) + CS%cg1(i,j+1))) @@ -166,11 +170,12 @@ subroutine calc_resoln_function(h, tv, G, CS) endif enddo ; enddo elseif (CS%Res_fn_power == 2) then +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 dx_term = CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j) CS%Res_fn_h(i,j) = dx_term / (dx_term + (CS%Res_coef * CS%cg1(i,j))**2) enddo ; enddo - +!$OMP do do J=js-1,Jeq ; do I=is-1,Ieq cg1_q = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + & (CS%cg1(i+1,j) + CS%cg1(i,j+1))) @@ -179,12 +184,13 @@ subroutine calc_resoln_function(h, tv, G, CS) enddo ; enddo elseif (mod_power_2 == 0) then power_2 = CS%Res_fn_power / 2 +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 dx_term = (CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j))**power_2 CS%Res_fn_h(i,j) = dx_term / & (dx_term + (CS%Res_coef * CS%cg1(i,j))**CS%Res_fn_power) enddo ; enddo - +!$OMP do do J=js-1,Jeq ; do I=is-1,Ieq cg1_q = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + & (CS%cg1(i+1,j) + CS%cg1(i,j+1))) @@ -193,6 +199,7 @@ subroutine calc_resoln_function(h, tv, G, CS) (dx_term + (CS%Res_coef * cg1_q)**CS%Res_fn_power) enddo ; enddo else +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 dx_term = (sqrt(CS%f2_dx2_h(i,j) + & CS%cg1(i,j)*CS%beta_dx2_h(i,j)))**CS%Res_fn_power @@ -200,7 +207,7 @@ subroutine calc_resoln_function(h, tv, G, CS) (dx_term + (CS%Res_coef * CS%cg1(i,j))**CS%Res_fn_power) enddo ; enddo - +!$OMP do do J=js-1,Jeq ; do I=is-1,Ieq cg1_q = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + & (CS%cg1(i+1,j) + CS%cg1(i,j+1))) @@ -213,10 +220,12 @@ subroutine calc_resoln_function(h, tv, G, CS) ! Calculate and store the ratio between deformation radius and grid-spacing ! at h-points (non-dimensional). +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 CS%Rd_dx_h(i,j) = CS%cg1(i,j) / & (sqrt(CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j))) enddo ; enddo +!$OMP end parallel if (query_averaging_enabled(CS%diag)) then if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 0fa843adc5..98674f04d1 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -315,7 +315,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, CS, & int_slope_u, int_slope_v) - +!$OMP parallel do default(none) shared(is,ie,js,je,nz,uhtr,uhD,dt,vhtr,CDp,vhD,h,G) do k=1,nz do j=js,je ; do I=is-1,ie uhtr(I,j,k) = uhtr(I,j,k) + uhD(I,j,k)*dt @@ -334,6 +334,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then +!$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix) do j=js,je ; do i=is,ie MEKE%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j) enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index c9b0bef716..f4ca070976 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -440,18 +440,29 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & reset_diags = .false. ! This is the second call to mixedlayer. if (reset_diags) then - if (CS%TKE_diagnostics) then ; do j=js,je ; do i=is,ie - CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_RiBulk(i,j) = 0.0 - CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_pen_SW(i,j) = 0.0 - CS%diag_TKE_mixing(i,j) = 0.0 ; CS%diag_TKE_mech_decay(i,j) = 0.0 - CS%diag_TKE_conv_decay(i,j) = 0.0 ; CS%diag_TKE_conv_s2(i,j) = 0.0 - enddo ; enddo ; endif - if (ALLOCATED(CS%diag_PE_detrain)) then ; do j=js,je ; do i=is,ie - CS%diag_PE_detrain(i,j) = 0.0 - enddo ; enddo ; endif - if (ALLOCATED(CS%diag_PE_detrain2)) then ; do j=js,je ; do i=is,ie - CS%diag_PE_detrain2(i,j) = 0.0 - enddo ; enddo ; endif +!$OMP parallel default(none) shared(is,ie,js,je,CS) + if (CS%TKE_diagnostics) then +!$OMP do + do j=js,je ; do i=is,ie + CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_RiBulk(i,j) = 0.0 + CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_pen_SW(i,j) = 0.0 + CS%diag_TKE_mixing(i,j) = 0.0 ; CS%diag_TKE_mech_decay(i,j) = 0.0 + CS%diag_TKE_conv_decay(i,j) = 0.0 ; CS%diag_TKE_conv_s2(i,j) = 0.0 + enddo ; enddo + endif + if (ALLOCATED(CS%diag_PE_detrain)) then +!$OMP do + do j=js,je ; do i=is,ie + CS%diag_PE_detrain(i,j) = 0.0 + enddo ; enddo + endif + if (ALLOCATED(CS%diag_PE_detrain2)) then +!$OMP do + do j=js,je ; do i=is,ie + CS%diag_PE_detrain2(i,j) = 0.0 + enddo ; enddo + endif +!$OMP end parallel endif if (CS%ML_resort) then diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 24f9cccda9..03f54ed030 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -313,7 +313,6 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! interpolated into depth space. logical :: showCallTree ! If true, show the call tree integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb - real, pointer :: T(:,:,:), S(:,:,:) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB nkmb = G%nk_rho_varies @@ -323,9 +322,6 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") - T => tv%T - S => tv%S - ! This sets the equivalence between the same bits of memory for these arrays. eaml => eatr ; ebml => ebtr @@ -337,7 +333,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debug) then call MOM_state_chksum("Start of diabatic ", u(:,:,:), v(:,:,:), h(:,:,:), G) endif - if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, G) call cpu_clock_begin(id_clock_set_diffusivity) call set_BBL_diffusivity(u, v, h, fluxes, visc, G, CS%set_diff_CSp) @@ -346,11 +342,11 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! Frazil formation keeps the temperature above the freezing point. ! make_frazil is deliberately called at both the beginning and at ! the end of the diabatic processes. - if (ASSOCIATED(T) .AND. ASSOCIATED(tv%frazil)) then + if (ASSOCIATED(tv%T) .AND. ASSOCIATED(tv%frazil)) then call make_frazil(h,tv,G,CS) if (showCallTree) call callTree_waypoint("done with 1st make_frazil (diabatic)") endif - if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then !$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml) @@ -363,7 +359,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call geothermal(h, tv, dt, eaml, ebml, G, CS%geothermal_CSp) call cpu_clock_end(id_clock_geothermal) if (showCallTree) call callTree_waypoint("geothermal (diabatic)") - if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) endif ! Set_opacity estimates the optical properties of the water column. @@ -396,12 +392,12 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! Keep salinity from falling below a small but positive threshold. ! This occurs when the ice model attempts to extract more salt than ! is actually present in the ocean. - if (ASSOCIATED(S) .and. ASSOCIATED(tv%salt_deficit)) & + if (ASSOCIATED(tv%S) .and. ASSOCIATED(tv%salt_deficit)) & call adjust_salt(h, tv, G, CS) call cpu_clock_end(id_clock_mixedlayer) if (CS%debug) call MOM_state_chksum("After mixedlayer ", u, v, h, G) if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif endif @@ -557,7 +553,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call KPP_applyNonLocalTransport(CS%KPP_CSp, G, h, CS%KPP_NLTscalar, CS%netSalt, dt, tv%S, isSalt=.true.) call cpu_clock_end(id_clock_kpp) if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)") - if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, G) if (CS%debug) then call MOM_state_chksum("after KPP_applyNLT ", u(:,:,:), v(:,:,:), h(:,:,:), G) @@ -568,13 +564,13 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! When used with KPP this needs to provide a diffusivity and happen before KPP ????? if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. & - associated(T)) then + associated(tv%T)) then call cpu_clock_begin(id_clock_differential_diff) ! Changes: tv%T, tv%S call differential_diffuse_T_S(h, tv, visc, dt, G) call cpu_clock_end(id_clock_differential_diff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") - if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) endif ! This block sets ea, eb from Kd or Kd_int. @@ -630,7 +626,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call hchksum(G%H_to_m*ea, "after applyBoundaryFluxes ea",G,haloshift=0) endif if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") - if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) endif ! Update h according to divergence of the difference between @@ -665,28 +661,28 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call MOM_thermovar_chksum("after negative check ", tv, G) endif if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)") - if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, hold, T, S, G) + if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, hold, tv%T, tv%S, G) ! Here, T and S are updated according to ea and eb. ! If using the bulk mixed layer, T and S are also updated ! by surface fluxes (in fluxes%*). if (CS%bulkmixedlayer) then - if (ASSOCIATED(T)) then + if (ASSOCIATED(tv%T)) then call cpu_clock_begin(id_clock_tridiag) ! Temperature and salinity (as state variables) are treated slightly ! differently from other tracers to insure that massless layers that ! are lighter than the mixed layer have temperatures and salinities ! that correspond to their prescribed densities. if (CS%massless_match_targets) then -!$OMP parallel do default (none) shared(is,ie,js,je,nkmb,hold,h_neglect,eb,T,ea,S,nz,kb) & +!$OMP parallel do default (none) shared(is,ie,js,je,nkmb,hold,h_neglect,eb,ea,nz,kb,tv) & !$OMP private(h_tr,b1,d1,c1,b_denom_1) do j=js,je do i=is,ie h_tr = hold(i,j,1) + h_neglect b1(i) = 1.0 / (h_tr + eb(i,j,1)) d1(i) = h_tr * b1(i) - T(i,j,1) = b1(i) * (h_tr*T(i,j,1)) - S(i,j,1) = b1(i) * (h_tr*S(i,j,1)) + tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1)) + tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1)) enddo do k=2,nkmb ; do i=is,ie c1(i,k) = eb(i,j,k-1) * b1(i) @@ -694,8 +690,8 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) b_denom_1 = h_tr + d1(i)*ea(i,j,k) b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) if (k kb(i,j)) then c1(i,k) = eb(i,j,k-1) * b1(i) h_tr = hold(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i)*ea(i,j,k) b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) d1(i) = b_denom_1 * b1(i) - T(i,j,k) = b1(i) * (h_tr*T(i,j,k) + ea(i,j,k)*T(i,j,k-1)) - S(i,j,k) = b1(i) * (h_tr*S(i,j,k) + ea(i,j,k)*S(i,j,k-1)) + tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1)) + tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1)) elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j)) ! The bottommost buffer layer might entrain all the mass from some ! of the interior layers that are thin and lighter in the coordinate ! density than that buffer layer. The T and S of these newly ! massless interior layers are unchanged. - T(i,j,nkmb) = T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * T(i,j,k) - S(i,j,nkmb) = S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * S(i,j,k) + tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k) + tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k) endif enddo ; enddo do k=nz-1,nkmb,-1 ; do i=is,ie if (k >= kb(i,j)) then - T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1) - S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1) + tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) + tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) endif enddo ; enddo do i=is,ie ; if (kb(i,j) <= nz) then - T(i,j,nkmb) = T(i,j,nkmb) + c1(i,kb(i,j))*T(i,j,kb(i,j)) - S(i,j,nkmb) = S(i,j,nkmb) + c1(i,kb(i,j))*S(i,j,kb(i,j)) + tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j)) + tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j)) endif ; enddo do k=nkmb-1,1,-1 ; do i=is,ie - T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1) - S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1) + tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) + tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) enddo ; enddo enddo ! end of j loop else ! .not. massless_match_targets ! This simpler form allows T & S to be too dense for the layers ! between the buffer layers and the interior. ! Changes: T, S - call triDiagTS(G, is, ie, js, je, hold, ea, eb, T, S) + call triDiagTS(G, is, ie, js, je, hold, ea, eb, tv%T, tv%S) endif ! massless_match_targets call cpu_clock_end(id_clock_tridiag) - if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G) endif ! end of ASSOCIATED(T) if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then @@ -788,28 +784,28 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! Keep salinity from falling below a small but positive threshold. ! This occurs when the ice model attempts to extract more salt than ! is actually present in the ocean. - if (ASSOCIATED(S) .and. ASSOCIATED(tv%salt_deficit)) & + if (ASSOCIATED(tv%S) .and. ASSOCIATED(tv%salt_deficit)) & call adjust_salt(h, tv, G, CS) call cpu_clock_end(id_clock_mixedlayer) if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif else ! Not BULKMIXEDLAYER. ! Calculate the change in temperature & salinity due to entrainment. - if (ASSOCIATED(T)) then + if (ASSOCIATED(tv%T)) then if (CS%debug) then call hchksum(G%H_to_m*ea, "before triDiagTS ea ",G,haloshift=0) call hchksum(G%H_to_m*eb, "before triDiagTS eb ",G,haloshift=0) endif call cpu_clock_begin(id_clock_tridiag) ! Changes: T, S - call triDiagTS(G, is, ie, js, je, hold, ea, eb, T, S) + call triDiagTS(G, is, ie, js, je, hold, ea, eb, tv%T, tv%S) call cpu_clock_end(id_clock_tridiag) if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") - if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) endif endif ! end BULKMIXEDLAYER @@ -823,7 +819,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call regularize_layers(h, tv, dt, ea, eb, G, CS%regularize_layers_CSp) call cpu_clock_end(id_clock_remap) if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)") - if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) endif if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & @@ -832,12 +828,12 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) Tdif_flx(i,j,1) = 0.0 ; Tdif_flx(i,j,nz+1) = 0.0 Tadv_flx(i,j,1) = 0.0 ; Tadv_flx(i,j,nz+1) = 0.0 enddo ; enddo -!$OMP parallel do default(none) shared(is,ie,js,je,nz,Tdif_flx,Idt,ea,eb,T,Tadv_flx) +!$OMP parallel do default(none) shared(is,ie,js,je,nz,Tdif_flx,Idt,ea,eb,Tadv_flx,tv) do K=2,nz ; do j=js,je ; do i=is,ie Tdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & - (T(i,j,k-1) - T(i,j,k)) + (tv%T(i,j,k-1) - tv%T(i,j,k)) Tadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & - 0.5*(T(i,j,k-1) + T(i,j,k)) + 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k)) enddo ; enddo ; enddo endif if ((CS%id_Sdif > 0) .or. (CS%id_Sdif_z > 0) .or. & @@ -846,12 +842,12 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) Sdif_flx(i,j,1) = 0.0 ; Sdif_flx(i,j,nz+1) = 0.0 Sadv_flx(i,j,1) = 0.0 ; Sadv_flx(i,j,nz+1) = 0.0 enddo ; enddo -!$OMP parallel do default(none) shared(is,ie,js,je,nz,Sdif_flx,Idt,ea,eb,S,Sadv_flx) +!$OMP parallel do default(none) shared(is,ie,js,je,nz,Sdif_flx,Idt,ea,eb,Sadv_flx,tv) do K=2,nz ; do j=js,je ; do i=is,ie Sdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & - (S(i,j,k-1) - S(i,j,k)) + (tv%S(i,j,k-1) - tv%S(i,j,k)) Sadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & - 0.5*(S(i,j,k-1) + S(i,j,k)) + 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) enddo ; enddo ; enddo endif @@ -936,9 +932,9 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) call cpu_clock_begin(id_clock_sponge) if (CS%bulkmixedlayer .and. ASSOCIATED(tv%eqn_of_state)) then do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo -!$OMP parallel do default(none) shared(js,je,T,S,p_ref_cv,Rcv_ml,is,ie,tv) +!$OMP parallel do default(none) shared(js,je,p_ref_cv,Rcv_ml,is,ie,tv) do j=js,je - call calculate_density(T(:,j,1), S(:,j,1), p_ref_cv, Rcv_ml(:,j), & + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & is, ie-is+1, tv%eqn_of_state) enddo call apply_sponge(h, dt, G, ea, eb, CS%sponge_CSp, Rcv_ml) @@ -1089,10 +1085,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) ! Frazil formation keeps the temperature above the freezing point. ! make_frazil is deliberately called at both the beginning and at ! the end of the diabatic processes. - if (ASSOCIATED(T) .AND. ASSOCIATED(tv%frazil)) then + if (ASSOCIATED(tv%T) .AND. ASSOCIATED(tv%frazil)) then call make_frazil(h,tv,G,CS) if (showCallTree) call callTree_waypoint("done with 2nd make_frazil (diabatic)") - if (CS%debugConservation) call MOM_state_stats('2nd make_frazil', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, G) endif if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) @@ -1140,7 +1136,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (num_z_diags > 0) & call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, CS%diag_to_Z_CSp) - if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, T, S, G) + if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G) if (showCallTree) call callTree_leave("diabatic()") end subroutine diabatic @@ -1198,8 +1194,8 @@ subroutine make_frazil(h, tv, G, CS) do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo endif !$OMP parallel do default(none) shared(is,ie,js,je,CS,G,h,nz,tv) & -!$OMP firstprivate(pressure) & -!$OMP private(fraz_col,T_fr_set,T_freeze,hc) +!$OMP private(fraz_col,T_fr_set,T_freeze,hc) & +!$OMP firstprivate(pressure) do j=js,je do i=is,ie ; fraz_col(:) = 0.0 ; enddo diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 5763fe23e1..46ef1ea21c 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -298,7 +298,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, CS) if ((G%nkml>0) .and. .not.use_BBL_EOS) then do i=G%IscB,G%IecB+1 ; p_ref(i) = tv%P_ref ; enddo - do k=1,nkmb ; do j=Jsq,Jeq+1 +!$OMP parallel do default(none) shared(Jsq,Jeq,Isq,Ieq,nkmb,tv,p_ref,Rml) + do j=Jsq,Jeq+1 ; do k=1,nkmb call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, & Rml(:,j,k), Isq, Ieq-Isq+2, tv%eqn_of_state) enddo ; enddo