Merge branch 'dev/gfdl' into opacity_rescale

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -1168,8 +1168,8 @@ subroutine apply_force_adjustments(G, US, CS, Time, forces)
   tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0
   ! Either reads data or leaves contents unchanged
   overrode_x = .false. ; overrode_y = .false.
-  call data_override(G%Domain, 'taux_adj', tempx_at_h, Time, override=overrode_x, scale=Pa_conversion)
-  call data_override(G%Domain, 'tauy_adj', tempy_at_h, Time, override=overrode_y, scale=Pa_conversion)
+  call data_override(G%Domain, 'taux_adj', tempx_at_h(isc:iec,jsc:jec), Time, override=overrode_x, scale=Pa_conversion)
+  call data_override(G%Domain, 'tauy_adj', tempy_at_h(isc:iec,jsc:jec), Time, override=overrode_y, scale=Pa_conversion)
 
   if (overrode_x .or. overrode_y) then
     if (.not. (overrode_x .and. overrode_y)) call MOM_error(FATAL,"apply_flux_adjustments: "//&

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

config_src/drivers/unit_drivers/MOM_sum_driver.F90

@@ -29,6 +29,7 @@ program MOM_main
   use MOM_io, only : MOM_io_init, file_exists, open_file, close_file
   use MOM_io, only : check_nml_error, io_infra_init, io_infra_end
   use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE
+  use MOM_unit_scaling, only : unit_scale_type, unit_no_scaling_init, unit_scaling_end
 
   implicit none
 
@@ -39,6 +40,8 @@ program MOM_main
   type(hor_index_type)   :: HI        ! A hor_index_type for array extents
   type(param_file_type)  :: param_file ! The structure indicating the file(s)
                                 ! containing all run-time parameters.
+  type(unit_scale_type), pointer :: US => NULL() !< A structure containing various unit
+                                ! conversion factors, but in this case all are 1.
   real    :: max_depth          ! The maximum ocean depth [m]
   integer :: verbosity
   integer :: num_sums
@@ -104,7 +107,8 @@ program MOM_main
   allocate(depth_tot_fastR(num_sums)) ; depth_tot_fastR(:) = 0.0
 
 ! Set up the parameters of the physical grid
-  call set_grid_metrics(grid, param_file)
+  call unit_no_scaling_init(US)
+  call set_grid_metrics(grid, param_file, US)
 
 ! Set up the bottom depth, grid%bathyT either analytically or from file
   call get_param(param_file, "MOM", "MAXIMUM_DEPTH", max_depth, &
@@ -162,6 +166,7 @@ program MOM_main
   enddo
 
   call destroy_dyn_horgrid(grid)
+  call unit_scaling_end(US)
   call io_infra_end ; call MOM_infra_end
 
 contains

config_src/infra/FMS1/MOM_cpu_clock_infra.F90

@@ -85,8 +85,13 @@ integer function cpu_clock_id(name, sync, grain)
   integer :: clock_flags
 
   clock_flags = clock_flag_default
-  if (present(sync)) &
-    clock_flags = ibset(clock_flags, 0)
+  if (present(sync)) then
+    if (sync) then
+      clock_flags = ibset(clock_flags, 0)
+    else
+      clock_flags = ibclr(clock_flags, 0)
+    endif
+  endif
 
   cpu_clock_id = mpp_clock_id(name, flags=clock_flags, grain=grain)
 end function cpu_clock_id

config_src/infra/FMS2/MOM_cpu_clock_infra.F90

@@ -85,8 +85,13 @@ integer function cpu_clock_id(name, sync, grain)
   integer :: clock_flags
 
   clock_flags = clock_flag_default
-  if (present(sync)) &
-    clock_flags = ibset(clock_flags, 0)
+  if (present(sync)) then
+    if (sync) then
+      clock_flags = ibset(clock_flags, 0)
+    else
+      clock_flags = ibclr(clock_flags, 0)
+    endif
+  endif
 
   cpu_clock_id = mpp_clock_id(name, flags=clock_flags, grain=grain)
 end function cpu_clock_id

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/core/MOM_checksum_packages.F90

@@ -92,23 +92,21 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric)
                                    intent(in) :: v  !< Meridional velocity [L T-1 ~> m s-1] or [m s-1]..
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  &
                                    intent(in) :: h  !< Layer thicknesses [H ~> m or kg m-2].
-  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type, which is
+  type(unit_scale_type),            intent(in) :: US !< A dimensional unit scaling type, which is
                                                     !! used to rescale u and v if present.
   integer,               optional, intent(in) :: haloshift !< The width of halos to check (default 0).
   logical,               optional, intent(in) :: symmetric !< If true, do checksums on the fully
                                                     !! symmetric computational domain.
-  real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> 1] or [1]
+
   integer :: hs
   logical :: sym
 
-  L_T_to_m_s = 1.0 ; if (present(US)) L_T_to_m_s = US%L_T_to_m_s
-
   ! Note that for the chksum calls to be useful for reproducing across PE
   ! counts, there must be no redundant points, so all variables use is..ie
   ! and js...je as their extent.
   hs = 1 ; if (present(haloshift)) hs = haloshift
   sym = .false. ; if (present(symmetric)) sym = symmetric
-  call uvchksum(mesg//" u", u, v, G%HI, haloshift=hs, symmetric=sym, scale=L_T_to_m_s)
+  call uvchksum(mesg//" u", u, v, G%HI, haloshift=hs, symmetric=sym, scale=US%L_T_to_m_s)
   call hchksum(h, mesg//" h",G%HI, haloshift=hs, scale=GV%H_to_m)
 end subroutine MOM_state_chksum_3arg
 

src/core/MOM_continuity_PPM.F90

@@ -223,7 +223,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are
                                                  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
   real,                    intent(in)    :: dt   !< Time increment [T ~> s].
   type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
-  type(continuity_PPM_CS), intent(in)       :: CS   !< This module's control structure.
+  type(continuity_PPM_CS), intent(in)    :: CS   !< This module's control structure.
   type(loop_bounds_type),  intent(in)    :: LB   !< Loop bounds structure.
   type(ocean_OBC_type),    pointer       :: OBC  !< Open boundaries control structure.
   real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), &
@@ -512,10 +512,10 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are
   if  (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then
     if (present(u_cor)) then
       call zonal_face_thickness(u_cor, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, &
-                                CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU(:,j,k), visc_rem_u)
+                                CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
     else
       call zonal_face_thickness(u, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, &
-                                CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU(:,j,k), visc_rem_u)
+                                CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU, visc_rem_u)
     endif
   endif ; endif
 
@@ -672,9 +672,11 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
     else ; h_u(I,j,k) = h_avg ; endif
   enddo ; enddo ; enddo
   if (present(visc_rem_u)) then
+    !### The expression setting h_u should also be multiplied by por_face_areaU in this case,
+    !    and in the two OBC cases below with visc_rem_u.
     !$OMP parallel do default(shared)
     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)
+      h_u(I,j,k) = h_u(I,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
     enddo ; enddo ; enddo
   endif
 
@@ -687,7 +689,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
         if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
           if (present(visc_rem_u)) then ; do k=1,nz
             do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
-              h_u(I,j,k) = h(i,j,k) * visc_rem_u(I,j,k)
+              h_u(I,j,k) = h(i,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
             enddo
           enddo ; else ; do k=1,nz
             do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
@@ -697,7 +699,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL,
         else
           if (present(visc_rem_u)) then ; do k=1,nz
             do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
-              h_u(I,j,k) = h(i+1,j,k) * visc_rem_u(I,j,k)
+              h_u(I,j,k) = h(i+1,j,k) * visc_rem_u(I,j,k) !### * por_face_areaU(I,j,k)
             enddo
           enddo ; else ; do k=1,nz
             do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed
@@ -1495,9 +1497,11 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
   enddo ; enddo ; enddo
 
   if (present(visc_rem_v)) then
+    !### This expression setting h_v should also be multiplied by por_face_areaU in this case,
+    !    and in the two OBC cases below with visc_rem_u.
     !$OMP parallel do default(shared)
     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)
+      h_v(i,J,k) = h_v(i,J,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
     enddo ; enddo ; enddo
   endif
 
@@ -1510,7 +1514,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
         if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
           if (present(visc_rem_v)) then ; do k=1,nz
             do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
-              h_v(i,J,k) = h(i,j,k) * visc_rem_v(i,J,k)
+              h_v(i,J,k) = h(i,j,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
             enddo
           enddo ; else ; do k=1,nz
             do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
@@ -1520,7 +1524,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL,
         else
           if (present(visc_rem_v)) then ; do k=1,nz
             do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied
-              h_v(i,J,k) = h(i,j+1,k) * visc_rem_v(i,J,k)
+              h_v(i,J,k) = h(i,j+1,k) * visc_rem_v(i,J,k) !### * por_face_areaV(i,J,k)
             enddo
           enddo ; else ; do k=1,nz
             do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied