From 94689e8cb6b08e7fc80afc6fbdc2e45965676bcb Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 11:46:03 -0400 Subject: [PATCH 01/11] Fix that is,ie,js,je might not be initialized --- config_src/solo_driver/MOM_surface_forcing.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 8119df43dd..52033c8137 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -1461,7 +1461,7 @@ subroutine buoyancy_forcing_const(state, fluxes, day, dt, G, CS) ! (in) G - The ocean's grid structure. ! (in) CS - A pointer to the control structure returned by a previous ! call to surface_forcing_init. - integer :: i, j, is, ie, js, je + integer :: i, j if ( CS%use_temperature ) then if (.not.associated(fluxes%evap)) then @@ -1510,7 +1510,7 @@ subroutine buoyancy_forcing_const(state, fluxes, day, dt, G, CS) endif if (associated(fluxes%p_surf)) then - do j=js,je ; do i=is,ie + do j=G%jsc,G%jec ; do i=G%isc,G%iec fluxes%p_surf(i,j) = 0.0 enddo ; enddo endif From d0ced28c74f3daf6fac754b6d9ba25e4e48ac0a3 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 15:56:17 -0400 Subject: [PATCH 02/11] Fix segmentation fault when compiling with -openmp and -O2. Possible due to Intel compiler bug --- src/equation_of_state/MOM_EOS_Wright.F90 | 44 ++++++++++++++++-------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 7b33516239..66e4bdf21d 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -405,22 +405,38 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, G, dza, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - do j=jsh,jeh ; do i=ish,ieh - al0_2d(i,j) = a0 + a1*T(i,j) + a2*S(i,j) - p0_2d(i,j) = b0 + b4*S(i,j) + T(i,j) * (b1 + T(i,j)*(b2 + b3*T(i,j)) + b5*S(i,j)) - lambda_2d(i,j) = c0 +c4*S(i,j) + T(i,j) * (c1 + T(i,j)*(c2 + c3*T(i,j)) + c5*S(i,j)) + if (present(intp_dza)) then + do j=jsh,jeh ; do i=ish,ieh + al0_2d(i,j) = a0 + a1*T(i,j) + a2*S(i,j) + p0_2d(i,j) = b0 + b4*S(i,j) + T(i,j) * (b1 + T(i,j)*(b2 + b3*T(i,j)) + b5*S(i,j)) + lambda_2d(i,j) = c0 +c4*S(i,j) + T(i,j) * (c1 + T(i,j)*(c2 + c3*T(i,j)) + c5*S(i,j)) + + al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) + dp = p_b(i,j) - p_t(i,j) + p_ave = 0.5*(p_t(i,j)+p_b(i,j)) - al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) - dp = p_b(i,j) - p_t(i,j) - p_ave = 0.5*(p_t(i,j)+p_b(i,j)) - - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - alpha_anom = al0 + lambda / (p0 + p_ave) - alpha_ref - rem = lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) - dza(i,j) = alpha_anom*dp + 2.0*eps*rem - if (present(intp_dza)) & + eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps + alpha_anom = al0 + lambda / (p0 + p_ave) - alpha_ref + rem = lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) + dza(i,j) = alpha_anom*dp + 2.0*eps*rem intp_dza(i,j) = 0.5*alpha_anom*dp**2 - dp*(1.0-eps)*rem - enddo ; enddo + enddo ; enddo + else + do j=jsh,jeh ; do i=ish,ieh + al0_2d(i,j) = a0 + a1*T(i,j) + a2*S(i,j) + p0_2d(i,j) = b0 + b4*S(i,j) + T(i,j) * (b1 + T(i,j)*(b2 + b3*T(i,j)) + b5*S(i,j)) + lambda_2d(i,j) = c0 +c4*S(i,j) + T(i,j) * (c1 + T(i,j)*(c2 + c3*T(i,j)) + c5*S(i,j)) + + al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) + dp = p_b(i,j) - p_t(i,j) + p_ave = 0.5*(p_t(i,j)+p_b(i,j)) + + eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps + alpha_anom = al0 + lambda / (p0 + p_ave) - alpha_ref + rem = lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) + dza(i,j) = alpha_anom*dp + 2.0*eps*rem + enddo ; enddo + endif if (present(intx_dza)) then ; do j=G%jsc,G%jec ; do I=Isq,Ieq intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) From 4a0d3e12b10f83f924cb7be291a6f140799c1c68 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:03:22 -0400 Subject: [PATCH 03/11] Initialize isym in routine calc_shelf_visc_triangular --- src/ice_shelf/MOM_ice_shelf.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 1685b3c95b..c2679c04b3 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -4887,6 +4887,12 @@ subroutine calc_shelf_visc_triangular (CS,u,v) G => CS%grid + if (G%symmetric) then + isym = 1 + else + isym = 0 + endif + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB gjsd = G%jsd_global ; gisd = G%isd_global From c386ac15fd48e1ac2cfd5b0f1fb6a0e6f33dc759 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:04:27 -0400 Subject: [PATCH 04/11] fix. hrat_min need to be private --- src/parameterizations/lateral/MOM_hor_visc.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 47ac85a814..28f67244d2 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -298,7 +298,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, CS, OBC) !$OMP parallel do default(shared) private(i, j, k, u0, v0, sh_xx, str_xx, & !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, & -!$OMP Shear_mag, huq, hvq, hq, Kh_scale) +!$OMP Shear_mag, huq, hvq, hq, Kh_scale, hrat_min) do k=1,nz ! This code uses boundary conditions that are consistent with From b7ef981c8cf52f3a918e69bb84539acd3d87f1ce Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:06:40 -0400 Subject: [PATCH 05/11] Fix that caused the model to crash when compiling with -openmp and -O2 or -O3. Implement openmp in this module --- .../lateral/MOM_thickness_diffuse.F90 | 180 +++++++++++++----- 1 file changed, 130 insertions(+), 50 deletions(-) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index a6340c7842..daa11a939b 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -100,7 +100,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) type(thermo_var_ptrs), intent(in) :: tv real, intent(in) :: dt type(ocean_grid_type), intent(in) :: G - type(MEKE_type), pointer :: MEKE + type(MEKE_type), intent(inout) :: MEKE type(VarMix_CS), pointer :: VarMix type(cont_diag_ptrs), intent(inout) :: CDp type(thickness_diffuse_CS), pointer :: CS @@ -117,7 +117,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) ! (in) G - The ocean's grid structure. ! (in) VarMix - A structure containing a number of fields related to ! variable lateral mixing. -! (in) MEKE - A structure containing information about the Mesoscale Eddy +! (inout) MEKE - A structure containing information about the Mesoscale Eddy ! Kinetic Energy parameterization; this might be unassociated. ! (inout) CDp - A structure with pointers to various terms in the continuity ! equations. @@ -142,22 +142,26 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) ! weighting of the interface slopes to that calculated also ! using density gradients at v points. The physically correct ! slopes occur at 0, while 1 is used for numerical closures. - real :: KH_u_CFL(SZIB_(G), SZJ_(G)) ! The maximum stable interface height - real :: KH_v_CFL(SZI_(G), SZJB_(G)) ! diffusivities at u & v grid points, - ! in m2 s-1. - real :: Khth_Loc ! Locally calculated thickness mixing coefficient (m2/s) + real, allocatable, save :: KH_u_CFL(:,:) ! The maximum stable interface height + real, allocatable, save :: KH_v_CFL(:,:) ! diffusivities at u & v grid points, + ! in m2 s-1. + logical, save :: first_call = .TRUE. + real :: Khth_Loc_u(SZIB_(G), SZJ_(G)) + real :: Khth_Loc(SZIB_(G), SZJB_(G)) ! Locally calculated thickness mixing coefficient (m2/s) logical :: use_VarMix, Resoln_scaled integer :: i, j, k, is, ie, js, je, nz - + logical :: MEKE_not_null if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse:"// & "Module must be initialized before it is used.") + MEKE_not_null = (LOC(MEKE) .NE. 0) + if ((.not.CS%thickness_diffuse) .or. & - .not.( CS%Khth > 0.0 .or. associated(VarMix) .or. associated(MEKE) ) ) return + .not.( CS%Khth > 0.0 .or. associated(VarMix) .or. MEKE_not_null ) ) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - if (ASSOCIATED(MEKE)) then + if (MEKE_not_null) then if (ASSOCIATED(MEKE%GM_src)) then do j=js,je ; do i=is,ie ; MEKE%GM_src(i,j) = 0. ; enddo ; enddo endif @@ -169,55 +173,124 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) Resoln_scaled = VarMix%Resoln_scaled_KhTh endif + if(first_call) then + allocate(KH_u_CFL(SZIB_(G), SZJ_(G)) ) + allocate(KH_v_CFL(SZI_(G), SZJB_(G)) ) +!$OMP parallel do default(shared) + do j = js,je; do I=is-1,ie + KH_u_CFL(I,j) = 0.2/(dt*(G%IdxCu(I,j)*G%IdxCu(I,j) + G%IdyCu(I,j)*G%IdyCu(I,j))) + enddo; enddo +!$OMP parallel do default(shared) + do j = js-1,je; do I=is,ie + KH_v_CFL(i,J) = 0.2/(dt*(G%IdxCv(i,J)*G%IdxCv(i,J) + G%IdyCv(i,J)*G%IdyCv(i,J))) + enddo; enddo + first_call = .false. + endif + call find_eta(h, tv, G%g_Earth, G, e, halo_size=1) ! Set the diffusivities. - do j=js,je ; do I=is-1,ie - Khth_Loc = CS%Khth - if (use_VarMix) & - Khth_Loc = Khth_Loc + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) - if (ASSOCIATED(MEKE)) then - if (associated(MEKE%Kh)) & - Khth_Loc = Khth_Loc + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) - endif - if (Resoln_scaled) & - Khth_Loc = Khth_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) - if (CS%Khth_Max > 0) then - Khth_Loc = max(CS%Khth_min, min(Khth_Loc,CS%Khth_Max)) - else - Khth_Loc = max(CS%Khth_min, Khth_Loc) - endif - KH_u_CFL(I,j) = 0.2/(dt*(G%IdxCu(I,j)*G%IdxCu(I,j) + G%IdyCu(I,j)*G%IdyCu(I,j))) - KH_u(I,j,1) = min(KH_u_CFL(I,j), Khth_Loc) +!$OMP parallel default(shared) +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = CS%Khth + enddo; enddo + + if(use_VarMix) then +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = Khth_Loc_u(i,j) + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j) + enddo; enddo + endif + + if(MEKE_not_null) then; + if (associated(MEKE%Kh)) then +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = Khth_Loc_u(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j)) + enddo; enddo + endif + endif + + if (Resoln_scaled) then +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = Khth_Loc_u(i,j) * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j)) + enddo; enddo + endif + + if (CS%Khth_Max > 0) then +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = max(CS%Khth_min, min(Khth_Loc_u(i,j),CS%Khth_Max)) + enddo; enddo + else +!$OMP do + do j = js,je; do I=is-1,ie + Khth_Loc_u(i,j) = max(CS%Khth_min, Khth_Loc_u(i,j)) + enddo; enddo + endif +!$OMP do + do j = js,je; do I=is-1,ie + KH_u(I,j,1) = min(KH_u_CFL(I,j), Khth_Loc_u(i,j)) enddo ; enddo + +!$OMP do do K=2,nz+1 ; do j=js,je ; do I=is-1,ie KH_u(I,j,K) = KH_u(I,j,1) enddo ; enddo ; enddo +!$OMP do do J=js-1,je ; do i=is,ie - Khth_Loc = CS%Khth - if (use_VarMix) & - Khth_Loc = Khth_Loc + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) - if (ASSOCIATED(MEKE)) then - if (associated(MEKE%Kh)) & - Khth_Loc = Khth_Loc + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) - endif - if (Resoln_scaled) & - Khth_Loc = Khth_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) - if (CS%Khth_Max > 0) then - Khth_Loc = max(CS%Khth_min, min(Khth_Loc,CS%Khth_Max)) - else - Khth_Loc = max(CS%Khth_min, Khth_Loc) + Khth_Loc(i,j) = CS%Khth + enddo; enddo + + if (use_VarMix) then +!$OMP do + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = Khth_Loc(i,j) + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J) + enddo; enddo + endif + if (MEKE_not_null) then + if (associated(MEKE%Kh)) then +!$OMP do + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = Khth_Loc(i,j) + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1)) + enddo; enddo endif - KH_v_CFL(i,J) = 0.2/(dt*(G%IdxCv(i,J)*G%IdxCv(i,J) + G%IdyCv(i,J)*G%IdyCv(i,J))) - KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_Loc) + endif + + if (Resoln_scaled) then +!$OMP do + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = Khth_Loc(i,j) * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1)) + enddo; enddo + endif + + if (CS%Khth_Max > 0) then +!$OMP do + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = max(CS%Khth_min, min(Khth_Loc(i,j),CS%Khth_Max)) + enddo; enddo + else +!$OMP do + do J=js-1,je ; do i=is,ie + Khth_Loc(i,j) = max(CS%Khth_min, Khth_Loc(i,j)) + enddo; enddo + endif +!$OMP do + do J=js-1,je ; do i=is,ie + KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_Loc(i,j)) enddo ; enddo +!$OMP do do K=2,nz+1 ; do J=js-1,je ; do i=is,ie KH_v(i,J,K) = KH_v(i,J,1) enddo ; enddo ; enddo - +!$OMP do do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo +!$OMP do do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo +!$OMP end parallel if (CS%detangle_interfaces) then call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, CS, & @@ -252,7 +325,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) enddo ; enddo enddo - if (ASSOCIATED(MEKE) .and. ASSOCIATED(VarMix)) then + if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then do j=js,je ; do i=is,ie MEKE%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j) @@ -292,7 +365,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & real, dimension(NIMEM_,NJMEMB_,NKMEM_), intent(out) :: vhD real, intent(in) :: dt type(ocean_grid_type), intent(in) :: G - type(MEKE_type), pointer :: MEKE + type(MEKE_type), intent(inout) :: MEKE type(thickness_diffuse_CS), pointer :: CS real, dimension(NIMEMB_,NJMEM_,NK_INTERFACE_), optional, intent(in) :: int_slope_u real, dimension(NIMEM_,NJMEMB_,NK_INTERFACE_), optional, intent(in) :: int_slope_v @@ -391,7 +464,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & ! Diagnostics that should be eliminated altogether later... ! real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: sfn_x, sfn_slope_x ! real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: sfn_y, sfn_slope_y - + logical :: MEKE_not_null integer :: is, ie, js, je, nz, IsdB integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; IsdB = G%IsdB @@ -403,18 +476,21 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & dz_neglect = G%H_subroundoff*G%H_to_m use_EOS = associated(tv%eqn_of_state) + MEKE_not_null = (LOC(MEKE) .NE. 0) nk_linear = max(G%nkml, 1) find_work = .false. - if (ASSOCIATED(MEKE)) find_work = ASSOCIATED(MEKE%GM_src) + if (MEKE_not_null) find_work = ASSOCIATED(MEKE%GM_src) find_work = (ASSOCIATED(CS%GMwork) .or. find_work) if (use_EOS) then call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth, dt, T, S, G, 1) endif +!$OMP parallel default(shared) ! Find the maximum and minimum permitted streamfunction. +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 h_avail_rsum(i,j,1) = 0.0 pres(i,j,1) = 0.0 ! ### This should be atmospheric pressure. @@ -424,8 +500,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & h_frac(i,j,1) = 1.0 pres(i,j,2) = pres(i,j,1) + G%H_to_Pa*h(i,j,1) enddo ; enddo - do k=2,nz - do j=js-1,je+1 ; do i=is-1,ie+1 +!$OMP do + do j=js-1,je+1 + do k=2,nz ; do i=is-1,ie+1 h_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-G%Angstrom),0.0) h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k) h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) & @@ -433,18 +510,21 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & pres(i,j,K+1) = pres(i,j,K) + G%H_to_Pa*h(i,j,k) enddo ; enddo enddo - +!$OMP do do j=js,je ; do I=is-1,ie uhtot(I,j) = 0.0 ; Work_u(I,j) = 0.0 ! sfn_x(I,j,1) = 0.0 ; sfn_slope_x(I,j,1) = 0.0 ! sfn_x(I,j,nz+1) = 0.0 ; sfn_slope_x(I,j,nz+1) = 0.0 enddo ; enddo +!$OMP do do J=js-1,je ; do i=is,ie vhtot(i,J) = 0.0 ; Work_v(i,J) = 0.0 ! sfn_y(i,J,1) = 0.0 ; sfn_slope_y(i,J,1) = 0.0 ! sfn_y(i,J,nz+1) = 0.0 ; sfn_slope_y(i,J,nz+1) = 0.0 enddo ; enddo +!$OMP end parallel +!!$OMP parallel do default(none) do K=nz,2,-1 if (find_work .and. .not.(use_EOS)) then drdiA = 0.0 ; drdiB = 0.0 ; drdjA = 0.0 ; drdjB = 0.0 @@ -822,7 +902,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, & Work_h = 0.5 * G%IareaT(i,j) * & ((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J))) if (ASSOCIATED(CS%GMwork)) CS%GMwork(i,j) = Work_h - if (ASSOCIATED(MEKE)) then ; if (ASSOCIATED(MEKE%GM_src)) then + if (MEKE_not_null) then ; if (ASSOCIATED(MEKE%GM_src)) then MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h endif ; endif enddo ; enddo ; endif From 204e261f8753d29fa0067d04408dbc69cebc27eb Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:29:03 -0400 Subject: [PATCH 06/11] Implement openmp. Add parameter allow_clocks_in_omp_loops to control clock inside openmp do loop. --- .../vertical/MOM_bulk_mixed_layer.F90 | 103 +++++++++++------- 1 file changed, 65 insertions(+), 38 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index cd9941e0ef..0f7d11e31d 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -172,7 +172,9 @@ module MOM_bulk_mixed_layer diag_PE_detrain2 ! The spurious source of potential ! energy due to mixed layer only ! detrainment, W m-2. - + logical :: allow_clocks_in_omp_loops ! If true, clocks can be called + ! from inside loops that can be threaded. + ! To run with multiple threads, set to False. integer :: id_ML_depth = -1, id_TKE_wind = -1, id_TKE_mixing = -1 integer :: id_TKE_RiBulk = -1, id_TKE_conv = -1, id_TKE_pen_SW = -1 integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1, id_TKE_conv_s2 = -1 @@ -180,8 +182,8 @@ module MOM_bulk_mixed_layer integer :: id_Hsfc_used = -1, id_Hsfc_max = -1, id_Hsfc_min = -1 end type bulkmixedlayer_CS -integer :: id_clock_detrain, id_clock_mech, id_clock_conv, id_clock_adjustment -integer :: id_clock_EOS, id_clock_resort, id_clock_pass +integer :: id_clock_detrain=0, id_clock_mech=0, id_clock_conv=0, id_clock_adjustment=0 +integer :: id_clock_EOS=0, id_clock_resort=0, id_clock_pass=0 integer :: num_msg = 0, max_msg = 2 @@ -407,21 +409,24 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & nsw = 0 if (ASSOCIATED(optics)) nsw = optics%nbands - allocate(Pen_SW_bnd(max(nsw,1),SZI_(G))) - allocate(opacity_band(max(nsw,1),SZI_(G),SZK_(G))) if (CS%limit_det .or. (CS%id_Hsfc_min > 0)) then +!$OMP parallel default(shared) +!$OMP do do j=js-1,je+1 ; do i=is-1,ie+1 h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0 enddo ; enddo - - do k=1,nkmb ; do j=js-1,je+1 ; do i=is-1,ie+1 - h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k) - hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k) - enddo ; enddo ; enddo - do k=nkmb+1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 - h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k) - enddo ; enddo ; enddo +!$OMP do + do j=js-1,je+1 + do k=1,nkmb ; do i=is-1,ie+1 + h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k) + hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k) + enddo ; enddo + do k=nkmb+1,nz ; do i=is-1,ie+1 + h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k) + enddo ; enddo + enddo +!$OMP end parallel call cpu_clock_begin(id_clock_pass) call pass_var(h_sum,G%Domain,complete=.false.) @@ -433,6 +438,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & reset_diags = .true. if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) & 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 @@ -454,6 +460,16 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & endif max_BL_det(:) = -1 +!$OMP parallel default(private) shared(is,ie,js,je,nz,h_3d,u_3d,v_3d,nkmb,G,nsw,optics, & +!$OMP CS,tv,fluxes,Irho0,dt,Idt_diag,Ih,write_diags, & +!$OMP hmbl_prev,h_sum,Hsfc_min,Hsfc_max,dt__diag, & +!$OMP Hsfc_used,Inkmlm1,Inkml,ea,eb,h_miss, & +!$OMP id_clock_EOS,id_clock_resort,id_clock_adjustment, & +!$OMP id_clock_conv,id_clock_mech,id_clock_detrain ) & +!$OMP firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) + allocate(Pen_SW_bnd(max(nsw,1),SZI_(G))) + allocate(opacity_band(max(nsw,1),SZI_(G),SZK_(G))) +!$OMP do do j=js,je ! Copy the thicknesses and other fields to 2-d arrays. do k=1,nz ; do i=is,ie @@ -470,7 +486,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 enddo ; enddo - call cpu_clock_begin(id_clock_EOS) + if(id_clock_EOS>0) call cpu_clock_begin(id_clock_EOS) ! Calculate an estimate of the mid-mixed layer pressure (in Pa) do i=is,ie ; p_ref(i) = 0.0 ; enddo do k=1,CS%nkml ; do i=is,ie @@ -486,21 +502,21 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), is, & ie-is+1, tv%eqn_of_state) enddo - call cpu_clock_end(id_clock_EOS) + if(id_clock_EOS>0) call cpu_clock_end(id_clock_EOS) if (CS%ML_resort) then - call cpu_clock_begin(id_clock_resort) + if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort) if (CS%ML_presort_nz_conv_adj > 0) & call convective_adjustment(h(:,1:), u, v, R0(:,1:), Rcv(:,1:), T(:,1:), & S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, CS, & CS%ML_presort_nz_conv_adj) call sort_ML(h(:,1:), R0(:,1:), eps, G, CS, ksort) - call cpu_clock_end(id_clock_resort) + if(id_clock_resort>0) call cpu_clock_end(id_clock_resort) else do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo - call cpu_clock_begin(id_clock_adjustment) + if(id_clock_adjustment>0) call cpu_clock_begin(id_clock_adjustment) ! Undergo instantaneous entrainment into the buffer layers and mixed layers ! to remove hydrostatic instabilities. Any water that is lighter than ! currently in the mixed or buffer layer is entrained. @@ -508,7 +524,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, CS) do i=is,ie ; h_CA(i) = h(i,1) ; enddo - call cpu_clock_end(id_clock_adjustment) + if(id_clock_adjustment>0) call cpu_clock_end(id_clock_adjustment) endif if (associated(fluxes%liq_runoff) .and. CS%do_rivermix) then @@ -533,7 +549,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & endif - call cpu_clock_begin(id_clock_conv) + if(id_clock_conv>0) call cpu_clock_begin(id_clock_conv) ! The surface forcing is contained in the fluxes type. ! Here, we "unpack" it and aggregate all the thermodynamic forcing into @@ -552,7 +568,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & nsw, Pen_SW_bnd, opacity_band, Conv_en, & dKE_FC, j, ksort, G, CS, tv) - call cpu_clock_end(id_clock_conv) + if(id_clock_conv>0) call cpu_clock_end(id_clock_conv) ! Now the mixed layer undergoes mechanically forced entrainment. ! The mixed layer may entrain down to the Monin-Obukhov depth if the @@ -560,7 +576,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & ! First the TKE at the depth of free convection that is available ! to drive mixing is calculated. - call cpu_clock_begin(id_clock_mech) + if(id_clock_mech>0) call cpu_clock_begin(id_clock_mech) call find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & @@ -579,7 +595,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & if (CS%TKE_diagnostics) then ; do i=is,ie CS%diag_TKE_mech_decay(i,j) = CS%diag_TKE_mech_decay(i,j) - Idt_diag*TKE(i) enddo ; endif - call cpu_clock_end(id_clock_mech) + if(id_clock_mech>0) call cpu_clock_end(id_clock_mech) ! Calculate the homogeneous mixed layer properties and store them in layer 0. do i=is,ie ; if (htot(i) > 0.0) then @@ -609,10 +625,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & ! these unused layers (but not currently in the code). if (CS%ML_resort) then - call cpu_clock_begin(id_clock_resort) + if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort) call resort_ML(h(:,0:), T(:,0:), S(:,0:), R0(:,0:), Rcv(:,0:), eps, & d_ea, d_eb, ksort, G, CS, dR0_dT, dR0_dS, dRcv_dT, dRcv_dS) - call cpu_clock_end(id_clock_resort) + if(id_clock_resort>0) call cpu_clock_end(id_clock_resort) endif if (CS%limit_det .or. (CS%id_Hsfc_max > 0) .or. (CS%id_Hsfc_min > 0)) then @@ -643,7 +659,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & ! Move water left in the former mixed layer into the buffer layer and ! from the buffer layer into the interior. These steps might best be ! treated in conjuction. - call cpu_clock_begin(id_clock_detrain) + if(id_clock_detrain>0) call cpu_clock_begin(id_clock_detrain) if (CS%nkbl == 1) then call mixedlayer_detrain_1(h(:,0:), T(:,0:), S(:,0:), R0(:,0:), Rcv(:,0:), & dt, dt__diag, d_ea, d_eb, j, G, CS, & @@ -656,7 +672,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & ! This code only works with 1 or 2 buffer layers. call MOM_error(FATAL, "MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.") endif - call cpu_clock_end(id_clock_detrain) + if(id_clock_detrain>0) call cpu_clock_end(id_clock_detrain) if (CS%id_Hsfc_used > 0) then @@ -764,6 +780,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & deallocate(Pen_SW_bnd) deallocate(opacity_band) +!$OMP end parallel if (write_diags) then if (CS%id_ML_depth > 0) & @@ -3384,6 +3401,12 @@ subroutine bulkmixedlayer_init(Time, G, param_file, diag, CS) "If true, use the fluxes%calving_Hflx field to set the \n"//& "heat carried by runoff, instead of using SST*CP*froz_runoff.", & default=.false.) + + call get_param(param_file, mod, "ALLOW_CLOCKS_IN_OMP_LOOPS", & + CS%allow_clocks_in_omp_loops, & + "If true, clocks can be called from inside loops that can \n"//& + "be threaded. To run with multiple threads, set to False.", & + default=.true.) CS%id_ML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, & Time, 'Surface mixed layer depth', 'meter') @@ -3445,20 +3468,24 @@ subroutine bulkmixedlayer_init(Time, G, param_file, diag, CS) if (CS%id_PE_detrain2 > 0) call safe_alloc_alloc(CS%diag_PE_detrain2, isd, ied, jsd, jed) if (CS%id_ML_depth > 0) call safe_alloc_alloc(CS%ML_depth, isd, ied, jsd, jed) - id_clock_detrain = cpu_clock_id('(Ocean mixed layer detrain)', grain=CLOCK_ROUTINE) - id_clock_mech = cpu_clock_id('(Ocean mixed layer mechanical entrainment)', grain=CLOCK_ROUTINE) - id_clock_conv = cpu_clock_id('(Ocean mixed layer convection)', grain=CLOCK_ROUTINE) - if (CS%ML_resort) then - id_clock_resort = cpu_clock_id('(Ocean mixed layer resorting)', grain=CLOCK_ROUTINE) - else - id_clock_adjustment = cpu_clock_id('(Ocean mixed layer convective adjustment)', grain=CLOCK_ROUTINE) + if(CS%allow_clocks_in_omp_loops) then + id_clock_detrain = cpu_clock_id('(Ocean mixed layer detrain)', grain=CLOCK_ROUTINE) + id_clock_mech = cpu_clock_id('(Ocean mixed layer mechanical entrainment)', grain=CLOCK_ROUTINE) + id_clock_conv = cpu_clock_id('(Ocean mixed layer convection)', grain=CLOCK_ROUTINE) + if (CS%ML_resort) then + id_clock_resort = cpu_clock_id('(Ocean mixed layer resorting)', grain=CLOCK_ROUTINE) + else + id_clock_adjustment = cpu_clock_id('(Ocean mixed layer convective adjustment)', grain=CLOCK_ROUTINE) + endif + id_clock_EOS = cpu_clock_id('(Ocean mixed layer EOS)', grain=CLOCK_ROUTINE) endif - id_clock_EOS = cpu_clock_id('(Ocean mixed layer EOS)', grain=CLOCK_ROUTINE) + if (CS%limit_det .or. (CS%id_Hsfc_min > 0)) & - id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=CLOCK_ROUTINE) + id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=CLOCK_ROUTINE) - if (CS%limit_det) then - endif + +! if (CS%limit_det) then +! endif end subroutine bulkmixedlayer_init From bf6cde4158bb86c12ff971284bdddfb18545e636 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:30:00 -0400 Subject: [PATCH 07/11] implement openmp --- .../vertical/MOM_entrain_diffusive.F90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index d886216e7d..14f6ccf5c8 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -279,7 +279,7 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, CS, ea, eb, & endif tolerance = G%m_to_H * CS%Tolerance_Ent - + g_2dt = 0.5 * G%g_Earth / dt kmb = G%nk_rho_varies K2 = max(kmb+1,2) ; kb_min = K2 if (.not. CS%bulkmixedlayer) then @@ -287,14 +287,26 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, CS, ea, eb, & ! These lines fill in values that are arbitrary, but needed because ! they are used to normalize the buoyancy flux in layer nz. do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ; enddo + else + kb(:) = 0 + do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo endif if (CS%id_diff_work > 0) allocate(diff_work(G%isd:G%ied,G%jsd:G%jed,nz+1)) if (CS%id_Kd > 0) allocate(Kd_eff(G%isd:G%ied,G%jsd:G%jed,nz)) correct_density = (CS%correct_density .and. associated(tv%eqn_of_state)) - if (correct_density) pres(:) = tv%P_Ref + if (correct_density) then + pres(:) = tv%P_Ref + else + pres(:) = 0.0 + endif +!$OMP parallel do default(private) shared(is,ie,js,je,nz,Kd_Lay,G,dt,Kd_int,CS,h,tv, & +!$OMP kmb,Angstrom,fluxes,K2,h_neglect,tolerance, & +!$OMP ea,eb,correct_density,Kd_eff,diff_work, & +!$OMP g_2dt, kb_out) & +!$OMP firstprivate(kb,ds_dsp1,dsp1_ds,pres,kb_min) do j=js,je do i=is,ie ; kb(i) = 1 ; enddo @@ -845,7 +857,6 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, CS, ea, eb, & endif if (CS%id_diff_work > 0) then - g_2dt = 0.5 * G%g_Earth / dt do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo if (associated(tv%eqn_of_state)) then if (associated(fluxes%p_surf)) then From 45c17c59fd1f2929927198d3dcda7154340e1136 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:30:51 -0400 Subject: [PATCH 08/11] openmp update --- src/parameterizations/vertical/MOM_kappa_shear.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 816c91cd80..ab89d6c0ac 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -323,7 +323,11 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, & !$OMP tv,G,CS,kappa_io,dz_massless,k0dt,p_surf,gR0,g_R0,dt, & -!$OMP tol_dksrc,tol_dksrc_low,tol2,Ri_crit,dt_refinements) +!$OMP tol_dksrc,tol_dksrc_low,tol2,Ri_crit,dt_refinements, & +#ifdef ADD_DIAGNOSTICS +!$OMP I_Ld2_3d,dz_Int_3d, & +#endif +!$OMP tke_io,kv_io) do j=js,je do k=1,nz ; do i=is,ie h_2d(i,k) = h(i,j,k)*G%H_to_m From 73e1a9ad012aa0aef8bf250cc3a6095348eb4c20 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:31:32 -0400 Subject: [PATCH 09/11] implement openmp --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index f0f019b360..b2d332ec06 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -546,6 +546,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, G, C if (CS%debug) call hchksum(Kd_sfc,"Kd_sfc",G,haloshift=0) +!$OMP parallel do default(private) shared(is,ie,js,je,nz,G,CS,h,tv,T_f,S_f,fluxes,dd, & +!$OMP Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc, & +!$OMP Kd_int,dt,u,v,Omega2) do j=js,je ! Set up variables related to the stratification. call find_N2(h, tv, T_f, S_f, fluxes, j, G, CS, dRho_int, N2_lay, N2_int, N2_bot) From dc4e11d58f908e63ed33c0289c1dd378580e91e6 Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Tue, 1 Apr 2014 16:32:01 -0400 Subject: [PATCH 10/11] openmp fix --- src/parameterizations/vertical/MOM_vert_friction.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index dc03c7f390..5347e42be5 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -956,6 +956,7 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m dz_neglect = G%H_subroundoff*G%H_to_m do_shelf = .false. ; if (present(shelf)) do_shelf = shelf + h_ml(:) = 0.0 ! The following loop calculates the vertical average velocity and ! surface mixed layer contributions to the vertical viscosity. From b57b2db81981800c235adb2e15a3c2b4ce8778aa Mon Sep 17 00:00:00 2001 From: Zhi Liang Date: Thu, 3 Apr 2014 10:41:38 -0400 Subject: [PATCH 11/11] fix minor bug --- src/tracer/MOM_tracer_registry.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index f1d7152111..0edef9262e 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -117,7 +117,7 @@ subroutine register_tracer(tr1, name, param_file, Reg, ad_x, ad_y, & ! (in,opt) df_2d_y - An array in which the vertically summed meridional diffusive ! fluxes are stored in units of CONC m3 s-1. - integer :: m, ntr + integer :: ntr type(tracer_type) :: temp character(len=256) :: mesg ! Message for error messages. @@ -143,10 +143,10 @@ subroutine register_tracer(tr1, name, param_file, Reg, ad_x, ad_y, & Reg%Tr(ntr)%OBC_in_u => OBC_in_u ; endif if (present(OBC_in_v)) then ; if (associated(OBC_in_v)) & Reg%Tr(ntr)%OBC_in_v => OBC_in_v ; endif - if (present(ad_2d_x)) then ; if (associated(ad_2d_x)) Reg%Tr(m)%ad2d_x => ad_2d_x ; endif - if (present(ad_2d_y)) then ; if (associated(ad_2d_y)) Reg%Tr(m)%ad2d_y => ad_2d_y ; endif - if (present(df_2d_x)) then ; if (associated(df_2d_x)) Reg%Tr(m)%df2d_x => df_2d_x ; endif - if (present(df_2d_y)) then ; if (associated(df_2d_y)) Reg%Tr(m)%df2d_y => df_2d_y ; endif + if (present(ad_2d_x)) then ; if (associated(ad_2d_x)) Reg%Tr(ntr)%ad2d_x => ad_2d_x ; endif + if (present(ad_2d_y)) then ; if (associated(ad_2d_y)) Reg%Tr(ntr)%ad2d_y => ad_2d_y ; endif + if (present(df_2d_x)) then ; if (associated(df_2d_x)) Reg%Tr(ntr)%df2d_x => df_2d_x ; endif + if (present(df_2d_y)) then ; if (associated(df_2d_y)) Reg%Tr(ntr)%df2d_y => df_2d_y ; endif end subroutine register_tracer