Merge branch 'dev/gfdl' into ice_dynamics

config_src/drivers/mct_cap/ocn_comp_mct.F90

@@ -722,7 +722,7 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn)
   call mct_gGrid_importRattr(dom_ocn,"lat",data,lsize)
 
   k = 0
-  L2_to_rad2 = grid%US%L_to_m**2 / grid%Rad_Earth**2
+  L2_to_rad2 = 1.0 / grid%Rad_Earth_L**2
   do j = grid%jsc, grid%jec
     do i = grid%isc, grid%iec
       k = k + 1 ! Increment position within gindex

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -1111,7 +1111,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
            k = k + 1 ! Increment position within gindex
            if (mask(k) /= 0) then
               mesh_areas(k) = dataPtr_mesh_areas(k)
-              model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
+              model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth_L**2
               mod2med_areacor(k) = model_areas(k) / mesh_areas(k)
               med2mod_areacor(k) = mesh_areas(k) / model_areas(k)
            end if

src/core/MOM.F90

@@ -188,10 +188,10 @@ module MOM
     vh, &           !< vh = v * h * dx at v grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
     vhtr            !< accumulated meridional thickness fluxes to advect tracers [H L2 ~> m3 or kg]
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ssh_rint
-                    !< A running time integral of the sea surface height [T m ~> s m].
+                    !< A running time integral of the sea surface height [T Z ~> s m].
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc
                     !< time-averaged (over a forcing time step) sea surface height
-                    !! with a correction for the inverse barometer [m]
+                    !! with a correction for the inverse barometer [Z ~> m]
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_av_bc
                     !< free surface height or column mass time averaged over the last
                     !! baroclinic dynamics time step [H ~> m or kg m-2]
@@ -515,7 +515,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
   logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
   real :: cycle_time    ! The length of the coupled time-stepping cycle [T ~> s].
   real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
-    ssh         ! sea surface height, which may be based on eta_av [m]
+    ssh         ! sea surface height, which may be based on eta_av [Z ~> m]
 
   real, dimension(:,:,:), pointer :: &
     u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1]
@@ -868,7 +868,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
       ! Determining the time-average sea surface height is part of the algorithm.
       ! This may be eta_av if Boussinesq, or need to be diagnosed if not.
       CS%time_in_cycle = CS%time_in_cycle + dt
-      call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0, dZref=G%Z_ref)
+      call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, dZref=G%Z_ref)
       do j=js,je ; do i=is,ie
         CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j)
       enddo ; enddo
@@ -2867,11 +2867,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     endif
   endif
 
-  if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then
+  if (query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then
+    if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z /= US%m_to_Z_restart) ) then
+      Z_rescale = US%m_to_Z / US%m_to_Z_restart
+      do j=js,je ; do i=is,ie
+        CS%ave_ssh_ibc(i,j) = Z_rescale * CS%ave_ssh_ibc(i,j)
+      enddo ; enddo
+    endif
+  else
     if (CS%split) then
-      call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0, dZref=G%Z_ref)
+      call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, dZref=G%Z_ref)
     else
-      call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0, dZref=G%Z_ref)
+      call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, dZref=G%Z_ref)
     endif
   endif
   if (CS%split) deallocate(eta)
@@ -2977,7 +2984,7 @@ subroutine register_diags(Time, G, GV, US, IDs, diag)
       'Layer Thickness after the dynamics update', thickness_units, conversion=GV%H_to_MKS, &
       v_extensive=.true.)
   IDs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
-      Time, 'Instantaneous Sea Surface Height', 'm')
+      Time, 'Instantaneous Sea Surface Height', 'm', conversion=US%Z_to_m)
 end subroutine register_diags
 
 !> Set up CPU clock IDs for timing various subroutines.
@@ -3097,14 +3104,14 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
   type(ocean_grid_type),             intent(in)    :: G   !< ocean grid structure
   type(verticalGrid_type),           intent(in)    :: GV  !< ocean vertical grid structure
   type(unit_scale_type),             intent(in)    :: US  !< A dimensional unit scaling type
-  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: ssh !< time mean surface height [m]
+  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: ssh !< time mean surface height [Z ~> m]
   real, dimension(:,:),              pointer       :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa]
   logical,                           intent(in)    :: use_EOS !< If true, calculate the density for
                                                        !! the SSH correction using the equation of state.
 
   real :: Rho_conv(SZI_(G))  ! The density used to convert surface pressure to
                       ! a corrected effective SSH [R ~> kg m-3].
-  real :: IgR0        ! The SSH conversion factor from R L2 T-2 to m [m T2 R-1 L-2 ~> m Pa-1].
+  real :: IgR0        ! The SSH conversion factor from R L2 T-2 to Z [Z T2 R-1 L-2 ~> m Pa-1].
   logical :: calc_rho
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: i, j, is, ie, js, je
@@ -3119,12 +3126,13 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, &
                                tv%eqn_of_state, EOSdom)
         do i=is,ie
-          IgR0 = US%Z_to_m / (Rho_conv(i) * GV%g_Earth)
+          IgR0 = 1.0 / (Rho_conv(i) * GV%g_Earth)
           ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
         enddo
       else
+        IgR0 = 1.0 / (GV%Rho0 * GV%g_Earth)
         do i=is,ie
-          ssh(i,j) = ssh(i,j) + p_atm(i,j) * (US%Z_to_m / (GV%Rho0 * GV%g_Earth))
+          ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
         enddo
       endif
     enddo
@@ -3209,7 +3217,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
   sfc_state%S_is_absS = CS%tv%S_is_absS
 
   do j=js,je ; do i=is,ie
-    sfc_state%sea_lev(i,j) = US%m_to_Z*CS%ave_ssh_ibc(i,j)
+    sfc_state%sea_lev(i,j) = CS%ave_ssh_ibc(i,j)
   enddo ; enddo
 
   if (allocated(sfc_state%frazil) .and. associated(CS%tv%frazil)) then ; do j=js,je ; do i=is,ie

src/diagnostics/MOM_diagnostics.F90

@@ -209,18 +209,23 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
 
 
   ! tmp array for surface properties
-  real :: surface_field(SZI_(G),SZJ_(G))
+  real :: surface_field(SZI_(G),SZJ_(G)) ! The surface temperature or salinity [degC] or [ppt]
   real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa]
-  real :: wt, wt_p
-
+  real :: wt, wt_p ! The fractional weights of two successive values when interpolating from
+                   ! a list [nondim], scaled so that wt + wt_p = 1.
   real :: f2_h     ! Squared Coriolis parameter at to h-points [T-2 ~> s-2]
   real :: mag_beta ! Magnitude of the gradient of f [T-1 L-1 ~> s-1 m-1]
   real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2]
 
   integer :: k_list
 
-  real, dimension(SZK_(GV)) :: temp_layer_ave, salt_layer_ave
-  real :: thetaoga, soga, masso, tosga, sosga
+  real, dimension(SZK_(GV)) :: temp_layer_ave ! The average temperature in a layer [degC]
+  real, dimension(SZK_(GV)) :: salt_layer_ave ! The average salinity in a layer [degC]
+  real :: thetaoga  ! The volume mean potential temperature [degC]
+  real :: soga      ! The volume mean ocean salinity [ppt]
+  real :: masso     ! The total mass of the ocean [kg]
+  real :: tosga     ! The area mean sea surface temperature [degC]
+  real :: sosga     ! The area mean sea surface salinity [ppt]
 
   is  = G%isc  ; ie   = G%iec  ; js  = G%jsc  ; je  = G%jec
   Isq = G%IscB ; Ieq  = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -437,7 +442,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     do j=js,je ; do i=is,ie
       surface_field(i,j) = tv%T(i,j,1)
     enddo ; enddo
-    tosga = global_area_mean(surface_field, G)
+    tosga = global_area_mean(tv%T(:,:,1), G)
     call post_data(CS%id_tosga, tosga, CS%diag)
   endif
 
@@ -1240,7 +1245,8 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh)
   type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
   type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
   real, dimension(SZI_(G),SZJ_(G)), &
-                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
+                            intent(in) :: ssh !< Time mean surface height without corrections
+                                              !! for ice displacement [Z ~> m]
 
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: speed  ! The surface speed [L T-1 ~> m s-1]
@@ -1280,23 +1286,25 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
   real,                     intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s].
   type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
   type(thermo_var_ptrs),    intent(in) :: tv  !< A structure pointing to various thermodynamic variables
-  real, dimension(SZI_(G),SZJ_(G)), &
-                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
+  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< Time mean surface height without corrections
+                                              !! for ice displacement [Z ~> m]
   real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
-                                                  !! for ice displacement and the inverse barometer [m]
+                                              !! for ice displacement and the inverse barometer [Z ~> m]
 
   real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
   real, dimension(SZI_(G),SZJ_(G)) :: &
     zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m]
   real :: I_time_int    ! The inverse of the time interval [T-1 ~> s-1].
-  real :: zos_area_mean, volo, ssh_ga
+  real :: zos_area_mean ! Global area mean sea surface height [m]
+  real :: volo          ! Total volume of the ocean [m3]
+  real :: ssh_ga        ! Global ocean area weighted mean sea seaface height [m]
   integer :: i, j, is, ie, js, je
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
   ! area mean SSH
   if (IDs%id_ssh_ga > 0) then
-    ssh_ga = global_area_mean(ssh, G)
+    ssh_ga = global_area_mean(ssh, G, scale=US%Z_to_m)
     call post_data(IDs%id_ssh_ga, ssh_ga, diag)
   endif
 
@@ -1306,7 +1314,7 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
   if (IDs%id_zos > 0 .or. IDs%id_zossq > 0) then
     zos(:,:) = 0.0
     do j=js,je ; do i=is,ie
-      zos(i,j) = ssh_ibc(i,j)
+      zos(i,j) = US%Z_to_m*ssh_ibc(i,j)
     enddo ; enddo
     zos_area_mean = global_area_mean(zos, G)
     do j=js,je ; do i=is,ie
@@ -1324,9 +1332,9 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
   ! post total volume of the liquid ocean
   if (IDs%id_volo > 0) then
     do j=js,je ; do i=is,ie
-      work_2d(i,j) = G%mask2dT(i,j)*(ssh(i,j) + US%Z_to_m*G%bathyT(i,j))
+      work_2d(i,j) = G%mask2dT(i,j) * (ssh(i,j) + G%bathyT(i,j))
     enddo ; enddo
-    volo = global_area_integral(work_2d, G)
+    volo = global_area_integral(work_2d, G, scale=US%Z_to_m)
     call post_data(IDs%id_volo, volo, diag)
   endif
 
@@ -1841,7 +1849,7 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv)
       standard_name='square_of_sea_surface_height_above_geoid',             &
       long_name='Square of sea surface height above geoid', units='m2')
   IDs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, Time, &
-      'Sea Surface Height', 'm')
+      'Sea Surface Height', 'm', conversion=US%Z_to_m)
   IDs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', Time, diag,&
       long_name='Area averaged sea surface height', units='m',            &
       standard_name='area_averaged_sea_surface_height')