Comment and whitespace fixes. mom-ocean#189

src/core/MOM_PressureForce_Montgomery.F90

@@ -62,6 +62,7 @@ module MOM_PressureForce_Mont
 use MOM_variables, only : thermo_var_ptrs
 use MOM_EOS, only : calculate_density, calculate_density_derivs
 use MOM_EOS, only : int_specific_vol_dp, query_compressible
+
 implicit none ; private
 
 #include <MOM_memory.h>
@@ -689,7 +690,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
   real, dimension(NIMEM_,NJMEM_,NKMEM_), optional, intent(in) :: rho_star
 !    This subroutine determines the partial derivative of the acceleration due 
 !  to pressure forces with the free surface height.
-! Arguments: e - Interface height, in m.
+! Arguments: e - Interface height, in H.
 !  (in)      tv - A structure containing pointers to any available
 !                 thermodynamic fields, including potential temperature and
 !                 salinity or mixed layer density. Absent fields have NULL ptrs.
@@ -703,7 +704,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
 !  (in)      CS - The control structure returned by a previous call to
 !                 PressureForce_init.
 !  (out)     pbce - the baroclinic pressure anomaly in each layer
-!                   due to free surface height anomalies, in m s-2.
+!                   due to free surface height anomalies, in m2 H-1 s-2.
 !  (in,opt)  rho_star - The layer densities (maybe compressibility compensated),
 !                       times g/rho_0, in m s-2.
 

src/core/MOM_PressureForce_analytic_FV.F90

@@ -490,7 +490,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
 !  (in)      p_atm - the pressure at the ice-ocean or atmosphere-ocean
 !                    interface in Pa.
 !  (out)     pbce - the baroclinic pressure anomaly in each layer
-!                   due to free surface height anomalies, in m s-2.
+!                   due to free surface height anomalies, in m2 H-1 s-2.
 !  (out)     eta - the free surface height used to calculate PFu and PFv, in m,
 !                  with any tidal contributions or compressibility compensation.
 
@@ -510,7 +510,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
     dpa_bk, &    ! The change in pressure anomaly between the top and bottom
                  ! of a layer, in Pa.
     intz_dpa_bk  ! The vertical integral in depth of the pressure anomaly less
-                 ! the pressure anomaly at the top of the layer, in m Pa.
+                 ! the pressure anomaly at the top of the layer, in H Pa.
   real, dimension(SZIB_BK_(G),SZJ_BK_(G)) :: & ! on block indices
     intx_pa_bk, & ! The zonal integral of the pressure anomaly along the interface
                   ! atop a layer, divided by the grid spacing, in Pa.

src/core/MOM_barotropic.F90

@@ -435,7 +435,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
 ! Arguments: U_in - The initial (3-D) zonal velocity, in m s-1.
 !  (in)      V_in - The initial (3-D) meridional velocity, in m s-1.
 !  (in)      eta_in - The initial barotropic free surface height anomaly or
-!                     column mass anomaly, in m or kg m-2.
+!                     column mass anomaly, in H.
 !  (in)      dt - The time increment to integrate over.
 !  (in)      bc_accel_u - The zonal baroclinic accelerations, in m s-2.
 !  (in)      bc_accel_v - The meridional baroclinic accelerations, in m s-2.
@@ -572,11 +572,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
                   ! spacing, in H m.
   real, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
     eta, &        ! The barotropic free surface height anomaly or column mass
-                  ! anomaly, in m or kg m-2.
-    eta_pred      ! A predictor value of eta, in m or kg m-2 like eta.
+                  ! anomaly, in H
+    eta_pred      ! A predictor value of eta, in H like eta.
   real, pointer, dimension(:,:) :: &
     eta_PF_BT     ! A pointer to the eta array (either eta or eta_pred) that
-                  ! determines the barotropic pressure force, in m or kg m-2.
+                  ! determines the barotropic pressure force, in H
   real, dimension(SZIW_(CS),SZJW_(CS)) :: &
     eta_sum, &    ! eta summed across the timesteps, in m or kg m-2.
     eta_wtd, &    ! A weighted estimate used to calculate eta_out, in m or kg m-2.
@@ -1764,6 +1764,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
                      (eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j)) * &
                     dgeo_de * CS%IdxCu(I,j)
       enddo ; enddo
+
       if (CS%dynamic_psurf) then
 !$OMP do
         do j=jsv,jev ; do I=isv-1,iev
@@ -2200,7 +2201,7 @@ subroutine set_dtbt(G, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
 !  (in)      CS - The control structure returned by a previous call to
 !                 barotropic_init.
 !  (in,opt)  eta - The barotropic free surface height anomaly or
-!                  column mass anomaly, in m or kg m-2.
+!                  column mass anomaly, in H
 !  (in,opt)  pbce - The baroclinic pressure anomaly in each layer
 !                   due to free surface height anomalies, in m2 H-1 s-2.
 !  (in,opt)  BT_cont - A structure with elements that describe the effective
@@ -3466,7 +3467,7 @@ subroutine find_face_areas(Datu, Datv, G, CS, MS, eta, halo, add_max)
 !                 barotropic_init.
 !  (in)      MS - A type that describes the memory sizes of the argument arrays.
 !  (in, opt) eta - The barotropic free surface height anomaly or
-!                  column mass anomaly, in m or kg m-2.
+!                  column mass anomaly, in H.
 !  (in, opt) halo - The halo size to use, default = 1.
 !  (in, opt) add_max - A value to add to the maximum depth (used to overestimate
 !                      the external wave speed) in m.

src/core/MOM_forcing_type.F90

@@ -660,8 +660,8 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, &
   type(thermo_var_ptrs),                 intent(inout) :: tv             ! thermodynamics type
   integer,                               intent(in)    :: j              ! j-row to work on 
   real, dimension(NIMEM_,NK_INTERFACE_), intent(inout) :: buoyancyFlux   ! buoyancy flux (m^2/s^3)
-  real, dimension(NIMEM_),               intent(inout) :: netHeatMinusSW ! surf Heat flux (K m/s)
-  real, dimension(NIMEM_),               intent(inout) :: netSalt        ! surf salt flux (ppt m/s)
+  real, dimension(NIMEM_),               intent(inout) :: netHeatMinusSW ! surf Heat flux (K H)
+  real, dimension(NIMEM_),               intent(inout) :: netSalt        ! surf salt flux (ppt H)
 
 
   ! Local variables
@@ -699,7 +699,7 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, &
   ! We aggregate the thermodynamic forcing for a time step into the following:
   ! netH       = water (H units) added/removed via surface fluxes
   ! netHeat    = heat (degC * H) via surface fluxes
-  ! netSalt    = salt ( g(salt)/m2 for non-Bouss and ppt*m/s for Bouss ) via surface fluxes
+  ! netSalt    = salt ( g(salt)/m2 for non-Bouss and ppt*m for Bouss ) via surface fluxes
   call extractFluxes1d(G, fluxes, optics, nsw, j, dt,                                 &
                 depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
                 h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW,                 &
@@ -748,8 +748,8 @@ subroutine calculateBuoyancyFlux2d(G, fluxes, optics, h, Temp, Salt, tv, &
   real, dimension(NIMEM_,NJMEM_,NKMEM_),       intent(in)    :: Salt           ! salinity (ppt)
   type(thermo_var_ptrs),                       intent(inout) :: tv             ! Thermodynamics type
   real, dimension(NIMEM_,NJMEM_,NK_INTERFACE_),intent(inout) :: buoyancyFlux   ! buoy flux (m^2/s^3)
-  real, dimension(NIMEM_,NJMEM_),optional,     intent(inout) :: netHeatMinusSW ! surf temp flux (K m/s)
-  real, dimension(NIMEM_,NJMEM_),optional,     intent(inout) :: netSalt        ! surf salt flux (ppt m/s)
+  real, dimension(NIMEM_,NJMEM_),optional,     intent(inout) :: netHeatMinusSW ! surf temp flux (K H)
+  real, dimension(NIMEM_,NJMEM_),optional,     intent(inout) :: netSalt        ! surf salt flux (ppt H)
 
   real, dimension( SZI_(G) ) :: netT ! net temperature flux (K m/s)
   real, dimension( SZI_(G) ) :: netS ! net saln flux (ppt m/s)

src/core/MOM_interface_heights.F90

@@ -199,10 +199,10 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size)
 !  (in)      tv      - structure pointing to various thermodynamic variables
 !  (in)      G_Earth - Earth gravitational acceleration (m/s2)
 !  (in)      G       - ocean grid structure
-!  (out)     eta     - free surface height relative to mean sea level (z=0) (m or kg/m2)
+!  (out)     eta     - free surface height relative to mean sea level (z=0) (m)
 !  (in,opt)  eta_bt  - optional barotropic variable that gives the "correct"
 !                      free surface height (Boussinesq) or total water column
-!                      mass per unit aread (non-Boussinesq), in m or kg m-2.
+!                      mass per unit aread (non-Boussinesq), in m
 !  (in,opt)  halo_size - width of halo points on which to calculate eta
 
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -112,8 +112,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
 !  limited to guarantee stability.
 
 ! Arguments: h - Layer thickness, in m.
-!  (in/out)  uhtr - Accumulated zonal mass fluxes in m3.
-!  (in/out)  vhtr - Accumulated meridional mass fluxes in m3.
+!  (in/out)  uhtr - Accumulated zonal mass fluxes in m2 H.
+!  (in/out)  vhtr - Accumulated meridional mass fluxes in m2 H.
 !  (in)      tv - A structure containing pointers to any available
 !                 thermodynamic fields. Absent fields have NULL ptrs.
 !  (in)      dt - Time increment in s.
@@ -131,7 +131,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
   real :: e(SZI_(G), SZJ_(G), SZK_(G)+1) ! The heights of the interfaces, relative
                                     ! to mean sea level, in m, pos. upward.
   real :: uhD(SZIB_(G), SZJ_(G), SZK_(G)) ! uhD & vhD are the diffusive u*h &
-  real :: vhD(SZI_(G), SZJB_(G), SZK_(G)) ! v*h fluxes, in m3 s-1.
+  real :: vhD(SZI_(G), SZJB_(G), SZK_(G)) ! v*h fluxes, in m2 H s-1.
 
   real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
     KH_u, &       ! The interface height diffusivities in u-columns, in m2 s-1.
@@ -266,14 +266,14 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
       enddo; enddo
     endif
   endif
-  
+
   if (Resoln_scaled) then
 !$OMP do
     do J=js-1,je ; do i=is,ie
       Khth_Loc(i,j) = Khth_Loc(i,j) * VarMix%Res_fn_v(i,J)
     enddo; enddo
   endif
-   
+
   if (CS%Khth_Max > 0) then
 !$OMP do
     do J=js-1,je ; do i=is,ie
@@ -285,6 +285,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
       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))
@@ -1017,7 +1018,7 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, CS
   type(thickness_diffuse_CS),            pointer       :: CS
   real, dimension(NIMEMB_,NJMEM_,NK_INTERFACE_), intent(inout) :: int_slope_u
   real, dimension(NIMEM_,NJMEMB_,NK_INTERFACE_), intent(inout) :: int_slope_v
-! Arguments: h - Layer thickness, in m.
+! Arguments: h - Layer thickness, in H.
 !  (in)      e - Interface heights relative to mean sea level, in m.
 !  (inout)   Kh_u - Thickness diffusivity on interfaces at u points, in m2 s-1.
 !  (inout)   Kh_v - Thickness diffusivity on interfaces at v points, in m2 s-1.

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -385,11 +385,11 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
   ! make_frazil is deliberately called at both the beginning and at
   ! the end of the diabatic processes.
   if (ASSOCIATED(tv%T) .AND. ASSOCIATED(tv%frazil)) then
-    if (ASSOCIATED(fluxes%p_surf_full)) then                                                                                                                                                       
-        call make_frazil(h,tv,G,CS,fluxes%p_surf_full)                                                                                                                                              
-    else                                                                                                                                                                                           
-        call make_frazil(h,tv,G,CS)                                                                                                                                                                 
-    endif   
+    if (ASSOCIATED(fluxes%p_surf_full)) then
+        call make_frazil(h,tv,G,CS,fluxes%p_surf_full)
+    else
+        call make_frazil(h,tv,G,CS)
+    endif
     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, tv%T, tv%S, G)
@@ -400,6 +400,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
       h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
     enddo ; enddo ; enddo
   endif
+
   if (CS%use_geothermal) then
     call cpu_clock_begin(id_clock_geothermal)
     call geothermal(h, tv, dt, eaml, ebml, G, CS%geothermal_CSp)
@@ -512,6 +513,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
     ! KPP needs the surface buoyancy flux but does not update state variables.
     ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.  ????
     ! Sets: CS%buoyancyFlux, CS%netHeatMinusSW, CS%netSalt
+
     call calculateBuoyancyFlux2d(G, fluxes, CS%optics, h, tv%T, tv%S, tv, &
                                  CS%buoyancyFlux, CS%netHeatMinusSW, CS%netSalt)
     ! The KPP scheme calculates the boundary layer diffusivities and non-local transport.
@@ -2463,7 +2465,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, &
     !                netMassOut < 0 means mass leaves ocean. 
     ! netHeat      = heat (degC * H) via surface fluxes, excluding the part 
     !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0. 
-    ! netSalt      = surface salt fluxes ( g(salt)/m2 for non-Bouss and ppt*m/s for Bouss )
+    ! netSalt      = surface salt fluxes ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
     ! Pen_SW_bnd   = components to penetrative shortwave radiation 
     call extractFluxes1d(G, fluxes, optics, nsw, j, dt,                        &
                   H_limit_fluxes, use_riverHeatContent, useCalvingHeatContent, &

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -736,6 +736,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, G, C
         Kd(i,j,k) = max( Kd(i,j,k) , &  ! Apply floor to Kd
            dissip * (CS%FluxRi_max / (G%Rho0 * (N2_lay(i,k) + Omega2))) )
       enddo ; enddo
+
       if (present(Kd_int)) then ; do K=2,nz ; do i=is,ie
       ! This calculates the dissipation ONLY from Kd calculated in this routine
       ! dissip has units of W/m3 (kg/m3 * m2/s * 1/s2 = J/s/m3)

src/parameterizations/vertical/MOM_shortwave_abs.F90

@@ -329,21 +329,21 @@ subroutine sumSWoverBands(G, h, opacity_band, nsw, j, dt, &
   real,                                  intent(in)    :: H_limit_fluxes
   logical,                               intent(in)    :: absorbAllSW
   real, dimension(:,:),                  intent(in)    :: iPen_SW_bnd
-  real, dimension(NIMEM_,NK_INTERFACE_), intent(inout) :: netPen ! Units of K m
+  real, dimension(NIMEM_,NK_INTERFACE_), intent(inout) :: netPen ! Units of K H
 
 ! Arguments:
 !  (in)      G             = ocean grid structure
 !  (in)      h             = layer thickness (m or kg/m^2)
 !                            units of h are referred to as H below.
 !  (in)      opacity_band  = opacity in each band of penetrating shortwave
-!                            radiation, in H-1. The indicies are band, i, k.
+!                            radiation, in m-1. The indicies are band, i, k.
 !  (in)      nsw           =  number of bands of penetrating shortwave radiation.
 !  (in)      j             = j-index to work on
 !  (in)      dt            = time step (seconds)
 !  (inout)   Pen_SW_bnd    = penetrating shortwave heating in each band that
 !                            hits the bottom and will be redistributed through
 !                            the water column (K H units) size nsw x NIMEM_.
-!  (out)     netPen        = attenuated flux at interfaces, summed over bands (K m units)
+!  (out)     netPen        = attenuated flux at interfaces, summed over bands (K H units)
 
   real :: h_heat(SZI_(G))     !  thickness of the water column that receives
                               !  remaining shortwave radiation, in H.

src/parameterizations/vertical/MOM_vert_friction.F90

@@ -569,7 +569,7 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, CS)
 
 ! Arguments: u - Zonal velocity, in m s-1.  Intent in.
 !  (in)      v - Meridional velocity, in m s-1.
-!  (in)      h - Layer thickness, in m.
+!  (in)      h - Layer thickness, in H.
 !  (in)      fluxes - A structure containing pointers to any possible
 !                     forcing fields.  Unused fields have NULL ptrs.
 !  (in)      visc - The vertical viscosity type, containing information about