From d7012046d9838a6398c9b382daac2930f30268f7 Mon Sep 17 00:00:00 2001 From: Nic Hannah Date: Thu, 9 Jul 2015 11:07:19 -0400 Subject: [PATCH] Unit conversions to support H vertical unit independant of m. No results change. Closes #189 --- src/core/MOM.F90 | 25 ++++++++++++++----- src/core/MOM_PressureForce_Montgomery.F90 | 22 ++++++++-------- src/core/MOM_PressureForce_analytic_FV.F90 | 21 ++++++++-------- src/core/MOM_barotropic.F90 | 16 ++++++------ src/core/MOM_dynamics_split_RK2.F90 | 13 +++++----- src/core/MOM_forcing_type.F90 | 6 ++--- src/core/MOM_grid.F90 | 16 +++++++----- src/core/MOM_interface_heights.F90 | 10 ++++---- src/diagnostics/MOM_sum_output.F90 | 2 +- .../MOM_state_initialization.F90 | 14 ++++++++++- src/parameterizations/vertical/MOM_KPP.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 6 ++--- .../vertical/MOM_set_diffusivity.F90 | 2 +- .../vertical/MOM_shortwave_abs.F90 | 2 +- 14 files changed, 94 insertions(+), 63 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c6017c7465..d0a7da8695 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -958,7 +958,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call check_redundant("Pre-ALE 1 ", u, v, G) endif call cpu_clock_begin(id_clock_ALE) + ! Switch thickness units from H to m for remapping + h = h*G%H_to_m call ALE_main(G, h, u, v, CS%tv, CS%ALE_CSp) + h = h*G%m_to_H call cpu_clock_end(id_clock_ALE) if (CS%debug) then call MOM_state_chksum("Post-ALE 1 ", u, v, h, CS%uh, CS%vh, G) @@ -1257,7 +1260,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call check_redundant("Pre-ALE ", u, v, G) endif call cpu_clock_begin(id_clock_ALE) + ! Switch thickness units from H to m for remapping + h = h*G%H_to_m call ALE_main(G, h, u, v, CS%tv, CS%ALE_CSp) + h = h*G%m_to_H call cpu_clock_end(id_clock_ALE) if (CS%debug) then call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G) @@ -1854,7 +1860,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! Allocate and initialize space for primary MOM variables. ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz)) ; CS%u(:,:,:) = 0.0 ALLOC_(CS%v(isd:ied,JsdB:JedB,nz)) ; CS%v(:,:,:) = 0.0 - ALLOC_(CS%h(isd:ied,jsd:jed,nz)) ; CS%h(:,:,:) = G%Angstrom + ALLOC_(CS%h(isd:ied,jsd:jed,nz)) ; CS%h(:,:,:) = G%Angstrom_z ALLOC_(CS%uh(IsdB:IedB,jsd:jed,nz)) ; CS%uh(:,:,:) = 0.0 ALLOC_(CS%vh(isd:ied,JsdB:JedB,nz)) ; CS%vh(:,:,:) = 0.0 if (CS%use_temperature) then @@ -1988,18 +1994,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) if (CS%debug) then call uchksum(CS%u,"Pre initialize_ALE u", G, haloshift=1) call vchksum(CS%v,"Pre initialize_ALE v", G, haloshift=1) - call hchksum(CS%h,"Pre initialize_ALE h", G, haloshift=1) + call hchksum(CS%h*G%H_to_m,"Pre initialize_ALE h", G, haloshift=1) endif + ! Switch thickness units from H to m for regridding + CS%h = CS%h*G%H_to_m + if (.not. query_initialized(CS%h,"h",CS%restart_CSp)) then ! This is a not a restart so we do the following... call adjustGridForIntegrity(CS%ALE_CSp, G, CS%h ) call ALE_main( G, CS%h, CS%u, CS%v, CS%tv, CS%ALE_CSp ) endif + + ! Switch thickness units back to H + CS%h = CS%h*G%m_to_H + call ALE_updateVerticalGridType( CS%ALE_CSp, G%GV ) if (CS%debug) then call uchksum(CS%u,"Post initialize_ALE u", G, haloshift=1) call vchksum(CS%v,"Post initialize_ALE v", G, haloshift=1) - call hchksum(CS%h, "Post initialize_ALE h", G, haloshift=1) + call hchksum(CS%h*G%H_to_m, "Post initialize_ALE h", G, haloshift=1) endif endif @@ -2652,8 +2665,8 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm) enddo do k=1,nz ; do i=is,ie - if (depth(i) + h(i,j,k) < depth_ml) then - dh = h(i,j,k) + if (depth(i) + h(i,j,k)*G%H_to_m < depth_ml) then + dh = h(i,j,k)*G%H_to_m elseif (depth(i) < depth_ml) then dh = depth_ml - depth(i) else @@ -2669,7 +2682,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm) enddo ; enddo ! Calculate the average properties of the mixed layer depth. do i=is,ie - if (depth(i) < G%H_subroundoff) depth(i) = G%H_subroundoff + if (depth(i) < G%H_subroundoff*G%H_to_m) depth(i) = G%H_subroundoff*G%H_to_m if (CS%use_temperature) then state%SST(i,j) = state%SST(i,j) / depth(i) state%SSS(i,j) = state%SSS(i,j) / depth(i) diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index bb862464f9..c7f1444f7b 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -516,7 +516,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) 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) + e(i,j,1) = e(i,j,1) + h(i,j,k)*G%H_to_m enddo ; enddo enddo call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp) @@ -537,7 +537,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) endif !$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) + e(i,j,K) = e(i,j,K+1) + h(i,j,k)*G%H_to_m enddo ; enddo ; enddo !$OMP end parallel if (use_EOS) then @@ -662,12 +662,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta) ! 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) + eta(i,j) = e(i,j,1)*G%m_to_H + e_tidal(i,j)*G%m_to_H 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) + eta(i,j) = e(i,j,1)*G%m_to_H enddo ; enddo endif endif @@ -729,7 +729,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star) Rho0xG = Rho0*g_Earth G_Rho0 = g_Earth/Rho0 use_EOS = associated(tv%eqn_of_state) - h_neglect = G%H_subroundoff * G%H_to_m + h_neglect = G%H_subroundoff*G%H_to_m if (use_EOS) then if (present(rho_star)) then @@ -737,8 +737,8 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star) !$OMP private(Ihtot) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect) - pbce(i,j,1) = GFS_scale * rho_star(i,j,1) + Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H) + pbce(i,j,1) = GFS_scale * rho_star(i,j,1) * (1.0 / G%m_to_H) enddo do k=2,nz ; do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * & @@ -750,13 +750,13 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star) !$OMP private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect) + Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H) press(i) = -Rho0xG*e(i,j,1) enddo call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, & Isq, Ieq-Isq+2, tv%eqn_of_state) do i=Isq,Ieq+1 - pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) + pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * (1.0 / G%m_to_H) enddo do k=2,nz do i=Isq,Ieq+1 @@ -779,8 +779,8 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star) !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h_neglect,pbce) private(Ihtot) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect) - pbce(i,j,1) = G%g_prime(1) + Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H) + pbce(i,j,1) = G%g_prime(1) * (1.0 / G%m_to_H) enddo do k=2,nz ; do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k-1) + & diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index eff2c2c09f..026ed7c415 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -576,13 +576,13 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! but that is not yet implemented, and the current form is correct for ! barotropic tides. !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h) - do j=Jsq,Jeq+1 + 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 + e(i,j,1) = e(i,j,1) + h(i,j,k)*G%H_to_m + enddo ; enddo enddo call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp) endif @@ -602,7 +602,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, endif !$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) + e(i,j,K) = e(i,j,K+1) + h(i,j,k)*G%H_to_m enddo ; enddo ; enddo !$OMP end parallel @@ -677,9 +677,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! of freedeom needed to know the linear profile). if ( use_ALE ) then if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then - call pressure_gradient_plm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h ); + call pressure_gradient_plm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h*G%H_to_m ); elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then - call pressure_gradient_ppm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h ); + call pressure_gradient_ppm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h*G%H_to_m ); endif endif @@ -754,6 +754,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, rho_ref, CS%Rho0, G%g_Earth, G%Block(n), tv%eqn_of_state, & dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk ) endif + intz_dpa_bk = intz_dpa_bk*G%m_to_H else do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk @@ -775,7 +776,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, PFu(I,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - & (pa_bk(ib+1,jb)*h(i+1,j,k) + intz_dpa_bk(ib+1,jb))) + & ((h(i+1,j,k) - h(i,j,k)) * intx_pa_bk(Ib,jb) - & - (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa_bk(Ib,jb))) * & + (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa_bk(Ib,jb) * G%m_to_H)) * & ((2.0*I_Rho0*G%IdxCu(I,j)) / & ((h(i,j,k) + h(i+1,j,k)) + h_neglect)) intx_pa_bk(Ib,jb) = intx_pa_bk(Ib,jb) + intx_dpa_bk(Ib,jb) @@ -786,7 +787,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, PFv(i,J,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - & (pa_bk(ib,jb+1)*h(i,j+1,k) + intz_dpa_bk(ib,jb+1))) + & ((h(i,j+1,k) - h(i,j,k)) * inty_pa_bk(ib,Jb) - & - (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa_bk(ib,Jb))) * & + (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa_bk(ib,Jb) * G%m_to_H)) * & ((2.0*I_Rho0*G%IdyCv(i,J)) / & ((h(i,j,k) + h(i,j+1,k)) + h_neglect)) inty_pa_bk(ib,Jb) = inty_pa_bk(ib,Jb) + inty_dpa_bk(ib,Jb) @@ -819,12 +820,12 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce, ! 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) + eta(i,j) = e(i,j,1)*G%m_to_H + e_tidal(i,j)*G%m_to_H 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) + eta(i,j) = e(i,j,1)*G%m_to_H enddo ; enddo endif endif diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 08a7d62cb6..146c90f298 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1375,7 +1375,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, & ! Limit the sink (inward) correction to the amount of mass that is already ! inside the cell, plus any mass added by eta_source. Htot = eta(i,j) - if (G%Boussinesq) Htot = CS%bathyT(i,j) + eta(i,j) + if (G%Boussinesq) Htot = CS%bathyT(i,j)*G%m_to_H + eta(i,j) CS%eta_cor(i,j) = max(CS%eta_cor(i,j), -max(0.0,Htot + dt*CS%eta_source(i,j))) endif @@ -2657,8 +2657,8 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, MS, halo, use_BT_cont, Datu, Datv, BT_OBC%Cg_u(I,j) = SQRT(G%g_prime(1)*(0.5* & (G%bathyT(i,j) + G%bathyT(i+1,j)))) if (G%Boussinesq) then - BT_OBC%H_u(I,j) = 0.5*((G%bathyT(i,j) + eta(i,j)) + & - (G%bathyT(i+1,j) + eta(i+1,j))) + BT_OBC%H_u(I,j) = 0.5*((G%bathyT(i,j)*G%m_to_H + eta(i,j)) + & + (G%bathyT(i+1,j)*G%m_to_H + eta(i+1,j))) else BT_OBC%H_u(I,j) = 0.5*(eta(i,j) + eta(i+1,j)) endif @@ -2693,8 +2693,8 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, MS, halo, use_BT_cont, Datu, Datv, BT_OBC%Cg_v(i,J) = SQRT(G%g_prime(1)*(0.5* & (G%bathyT(i,j) + G%bathyT(i,j+1)))) if (G%Boussinesq) then - BT_OBC%H_v(i,J) = 0.5*((G%bathyT(i,j) + eta(i,j)) + & - (G%bathyT(i,j+1) + eta(i,j+1))) + BT_OBC%H_v(i,J) = 0.5*((G%bathyT(i,j)*G%m_to_H + eta(i,j)) + & + (G%bathyT(i,j+1)*G%m_to_H + eta(i,j+1))) else BT_OBC%H_v(i,J) = 0.5*(eta(i,j) + eta(i,j+1)) endif @@ -3487,14 +3487,14 @@ subroutine find_face_areas(Datu, Datv, G, CS, MS, eta, halo, add_max) if (G%Boussinesq) then !$OMP do do j=js-hs,je+hs ; do I=is-1-hs,ie+hs - H1 = CS%bathyT(i,j) + eta(i,j) ; H2 = CS%bathyT(i+1,j) + eta(i+1,j) + H1 = CS%bathyT(i,j)*G%m_to_H + eta(i,j) ; H2 = CS%bathyT(i+1,j)*G%m_to_H + eta(i+1,j) Datu(I,j) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) & Datu(I,j) = CS%dy_Cu(I,j) * (2.0 * H1 * H2) / (H1 + H2) ! Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (H1 + H2) enddo ; enddo !$OMP do do J=js-1-hs,je+hs ; do i=is-hs,ie+hs - H1 = CS%bathyT(i,j) + eta(i,j) ; H2 = CS%bathyT(i,j+1) + eta(i,j+1) + H1 = CS%bathyT(i,j)*G%m_to_H + eta(i,j) ; H2 = CS%bathyT(i,j+1)*G%m_to_H + eta(i,j+1) Datv(i,J) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) & Datv(i,J) = CS%dx_Cv(i,J) * (2.0 * H1 * H2) / (H1 + H2) ! Datv(i,J) = CS%dy_v(i,J) * 0.5 * (H1 + H2) @@ -3597,7 +3597,7 @@ subroutine bt_mass_source(h, eta, fluxes, set_cor, dt_therm, dt_since_therm, & do j=js,je do i=is,ie ; h_tot(i) = h(i,j,1) ; enddo if (G%Boussinesq) then - do i=is,ie ; eta_h(i) = h(i,j,1) - G%bathyT(i,j) ; enddo + do i=is,ie ; eta_h(i) = h(i,j,1) - G%bathyT(i,j)*G%m_to_H ; enddo else do i=is,ie ; eta_h(i) = h(i,j,1) ; enddo endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 570424e17f..d5d35161ee 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -496,12 +496,13 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_begin(id_clock_pres) call PressureForce(h, tv, CS%PFu, CS%PFv, G, CS%PressureForce_CSp, & CS%ALE_CSp, p_surf, CS%pbce, CS%eta_PF) + if (CS%debug) then + call uchksum(CS%PFu,"After PressureForce CS%PFu",G,haloshift=0) + call vchksum(CS%PFv,"After PressureForce CS%PFv",G,haloshift=0) + endif + if (dyn_p_surf) then - if (G%Boussinesq) then - Pa_to_eta = 1.0 / (G%Rho0*G%g_Earth) - else - Pa_to_eta = 1.0 / G%H_to_Pa - endif + Pa_to_eta = 1.0 / G%H_to_Pa !$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 * & @@ -1242,7 +1243,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, & ! dimensions as h, either m or kg m-3. ! CS%eta(:,:) = 0.0 already from initialization. if (G%Boussinesq) then - do j=js,je ; do i=is,ie ; CS%eta(i,j) = -G%bathyT(i,j) ; enddo ; enddo + do j=js,je ; do i=is,ie ; CS%eta(i,j) = -G%bathyT(i,j) * G%m_to_H ; enddo ; enddo endif do k=1,nz ; do j=js,je ; do i=is,ie CS%eta(i,j) = CS%eta(i,j) + h(i,j,k) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 022dda86af..6556d1dc6b 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -687,7 +687,7 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, & useRiverHeatContent = .False. useCalvingHeatContent = .False. - depthBeforeScalingFluxes = max( G%Angstrom, 1.e-30 ) + depthBeforeScalingFluxes = max( G%Angstrom_z, 1.e-30 ) pressure(:) = 0. ! Ignore atmospheric pressure GoRho = G%g_Earth / G%Rho0 start = 1 + G%isc - G%isd @@ -724,10 +724,10 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, & ! Convert to a buoyancy flux, excluding penetrating SW heating buoyancyFlux(G%isc:G%iec,1) = - GoRho * ( dRhodS(G%isc:G%iec) * netSalt(G%isc:G%iec) + & - dRhodT(G%isc:G%iec) * netHeat(G%isc:G%iec) ) ! m^2/s^3 + dRhodT(G%isc:G%iec) * netHeat(G%isc:G%iec) ) * G%H_to_m ! m^2/s^3 ! We also have a penetrative buoyancy flux associated with penetrative SW do k=2, G%ke+1 - buoyancyFlux(G%isc:G%iec,k) = - GoRho * ( dRhodT(G%isc:G%iec) * netPen(G%isc:G%iec,k) ) ! m^2/s^3 + buoyancyFlux(G%isc:G%iec,k) = - GoRho * ( dRhodT(G%isc:G%iec) * netPen(G%isc:G%iec,k) ) * G%H_to_m ! m^2/s^3 enddo end subroutine calculateBuoyancyFlux1d diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index eee7f3732c..6027d2c7a1 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -236,11 +236,16 @@ subroutine MOM_grid_init(G, param_file) call get_param(param_file, "MOM_grid", "ANGSTROM", G%Angstrom_z, & "The minumum layer thickness, usually one-Angstrom.", & units="m", default=1.0e-10) - if (.not.G%Boussinesq) & + + if (.not.G%Boussinesq) then call get_param(param_file, "MOM_grid", "H_TO_KG_M2", G%H_to_kg_m2,& "A constant that translates thicknesses from the model's \n"//& "internal units of thickness to kg m-2.", units="kg m-2 H-1", & default=1.0) + else + call get_param(param_file, "MOM_grid", "H_TO_M", G%H_to_m, default=1.0) + endif + call get_param(param_file, "MOM_grid", "BATHYMETRY_AT_VEL", G%bathymetry_at_vel, & "If true, there are separate values for the basin depths \n"//& "at velocity points. Otherwise the effects of of \n"//& @@ -306,11 +311,10 @@ subroutine MOM_grid_init(G, param_file) endif if (G%Boussinesq) then - G%H_to_kg_m2 = G%Rho0 - G%kg_m2_to_H = 1.0/G%Rho0 - G%m_to_H = 1.0 - G%H_to_m = 1.0 - G%Angstrom = G%Angstrom_z + G%H_to_kg_m2 = G%Rho0 * G%H_to_m + G%kg_m2_to_H = 1.0 / G%H_to_kg_m2 + G%m_to_H = 1.0 / G%H_to_m + G%Angstrom = G%Angstrom_z*G%m_to_H else G%kg_m2_to_H = 1.0 / G%H_to_kg_m2 G%m_to_H = G%Rho0 * G%kg_m2_to_H diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index e117cd4c2a..4d333b664a 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -114,12 +114,12 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) !$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 + do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -G%bathyT(i,j); enddo ; enddo if (G%Boussinesq) then !$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) + eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)*G%H_to_m enddo ; enddo ; enddo if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height @@ -127,10 +127,10 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size) !$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)) + dilate(i) = (eta_bt(i,j)+G%bathyT(i,j)*G%m_to_H) / (eta(i,j,1)+G%bathyT(i,j)) enddo do k=1,nz ; do i=isv,iev - eta(i,j,K) = dilate(i) * (eta(i,j,K) + G%bathyT(i,j)) - G%bathyT(i,j) + eta(i,j,K) = dilate(i) * (eta(i,j,K) + G%bathyT(i,j)*G%m_to_H) - G%bathyT(i,j) enddo ; enddo enddo endif @@ -234,7 +234,7 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size) else !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie - eta(i,j) = eta(i,j) + h(i,j,k) + eta(i,j) = eta(i,j) + h(i,j,k)*G%H_to_m enddo ; enddo ; enddo endif else diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 5cda60ab51..18bfb93398 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -638,7 +638,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, CS, tracer_CSp) do j=js,je ; do i=is,ie hbelow = 0.0 do k=nz,1,-1 - hbelow = hbelow + h(i,j,k) + hbelow = hbelow + h(i,j,k) * G%H_to_m hint = (H_0APE(K) + hbelow - G%bathyT(i,j)) hbot = H_0APE(K) - G%bathyT(i,j) hbot = (hbot + ABS(hbot)) * 0.5 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 450da3ebdf..6174f2647f 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -325,7 +325,19 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, PF, dirs, & "units of m to kg m-2 or vice versa, depending on whether \n"//& "BOUSSINESQ is defined. This does not apply if a restart \n"//& "file is read.", default=.false.) - if (convert) call convert_thickness(h, G, PF, tv) + if (convert .and. .not. G%Boussinesq) then + ! Convert h from m to kg m-2 then to thickness units (H) + call convert_thickness(h, G, PF, tv) + elseif (G%Boussinesq) then + ! Convert h from m to thickness units (H) + h = h*G%m_to_H + endif + + if (debug) then + call hchksum(h*G%H_to_m, "MOM_initialize_state: h ", G, haloshift=1) + if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G, haloshift=1) + if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G, haloshift=1) + endif ! Remove the mass that would be displaced by an ice shelf or inverse barometer. call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, & diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index be40cf8151..a82dc56717 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -570,7 +570,7 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, & km1 = max(1, k-1) kk = 3*(k-1) deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1) - N2_1d(k) = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / (0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff) + N2_1d(k) = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / ((0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff)*G%H_to_m) N_1d(k) = sqrt( max( N2_1d(k), 0.) ) enddo N2_1d(G%ke+1 ) = 0.0 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 9a81cf5131..1bfeb343c0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2412,7 +2412,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, & ! To accommodate vanishing upper layers, we need to allow for an instantaneous ! distribution of forcing over some finite vertical extent. The bulk mixed layer ! code handles this issue properly. - H_limit_fluxes = max(G%Angstrom, 1.E-30) + H_limit_fluxes = max(G%Angstrom_z, 1.E-30) ! The inverse scale, IforcingDepthScale, is a hack which ! should not be tickled in Eulerian mode. It stops all of the forcing from @@ -2445,7 +2445,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, & h2d(i,k) = h(i,j,k) T2d(i,k) = tv%T(i,j,k) do n=1,nsw - opacityBand(n,i,k) = G%H_to_m*optics%opacity_band(n,i,j,k) + opacityBand(n,i,k) = (1.0 / G%m_to_H)*optics%opacity_band(n,i,j,k) enddo enddo ; enddo @@ -2589,7 +2589,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, & ! For layers thin relative to 1/IforcingDepthScale, then distribute ! forcing into deeper layers. ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale. - fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale) + fractionOfForcing = min(1.0, h2d(i,k)*G%H_to_m*IforcingDepthScale) ! In the case with (-1)*netMassOut greater than 0.8*h, then we limit ! applied to the top cell, and distribute the fluxes downwards. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 44ea4cf499..c728f9f217 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -1654,7 +1654,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, G, CS, ! This is energy loss in addition to work done as mixing, apparently to Joule heating. TKE_remaining = exp(-Idecay*dh) * TKE_remaining - z = z + h(i,j,k) ! Distance between upper interface of layer and the bottom, in m. + z = z + h(i,j,k)*G%H_to_m ! Distance between upper interface of layer and the bottom, in m. D_minus_z = max(total_thickness - z, 0.) ! Thickness above layer, m. ! Diffusivity using law of the wall, limited by rotation, at height z, in m2/s. diff --git a/src/parameterizations/vertical/MOM_shortwave_abs.F90 b/src/parameterizations/vertical/MOM_shortwave_abs.F90 index 80af047ff0..327d515b7f 100644 --- a/src/parameterizations/vertical/MOM_shortwave_abs.F90 +++ b/src/parameterizations/vertical/MOM_shortwave_abs.F90 @@ -385,7 +385,7 @@ subroutine sumSWoverBands(G, h, opacity_band, nsw, j, dt, & if (h(i,k) > 0.0) then do n=1,nsw ; if (Pen_SW_bnd(n,i) > 0.0) then ! SW_trans is the SW that is transmitted THROUGH the layer - opt_depth = h(i,k) * opacity_band(n,i,k) + opt_depth = h(i,k)*G%H_to_m * opacity_band(n,i,k) exp_OD = exp(-opt_depth) SW_trans = exp_OD ! Heating at a rate of less than 10-4 W m-2 = 10-3 K m / Century,