Merge branch 'dev/gfdl' into new_save_restart

.testing/tc3/MOM_input

@@ -283,10 +283,10 @@ HMIX_FIXED = 20.0               !   [m]
 KV = 1.0E-04                    !   [m2 s-1]
                                 ! The background kinematic viscosity in the interior.
                                 ! The molecular value, ~1e-6 m2 s-1, may be used.
-KVML = 0.01                     !   [m2 s-1] default = 1.0E-04
-                                ! The kinematic viscosity in the mixed layer.  A typical
-                                ! value is ~1e-2 m2 s-1. KVML is not used if
-                                ! BULKMIXEDLAYER is true.  The default is set by KV.
+KV_ML_INVZ2 = 0.01              !   [m2 s-1] default = 0.0
+                                ! An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, with
+                                ! the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the distance
+                                ! from the surface, to allow for finite wind stresses to be transmitted through.
 HBBL = 10.0                     !   [m]
                                 ! The thickness of a bottom boundary layer with a
                                 ! viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or

config_src/infra/FMS1/MOM_coms_infra.F90

@@ -25,7 +25,7 @@ module MOM_coms_infra
 !> Communicate an array, string or scalar from one PE to others
 interface broadcast
   module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D
-  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D
+  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D
 end interface broadcast
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is
@@ -260,6 +260,28 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking)
 
 end subroutine broadcast_real2D
 
+
+!> Communicate a 3-D array of reals from one PE to others
+subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking)
+  real, dimension(:,:,:), intent(inout) :: dat       !< The data to communicate and destination
+  integer,              intent(in)      :: length    !< The total number of data elements
+  integer,    optional, intent(in)      :: from_PE   !< The source PE, by default the root PE
+  integer,    optional, intent(in)      :: PElist(:) !< The list of participating PEs, by default the
+                                                     !! active PE set as previously set via Set_PElist.
+  logical,    optional, intent(in)      :: blocking  !< If true, barriers are added around the call
+
+  integer :: src_PE   ! The processor that is sending the data
+  logical :: do_block ! If true add synchronizing barriers
+
+  do_block = .false. ; if (present(blocking)) do_block = blocking
+  if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif
+
+  if (do_block) call mpp_sync(PElist)
+  call mpp_broadcast(dat, length, src_PE, PElist)
+  if (do_block) call mpp_sync_self(PElist)
+
+end subroutine broadcast_real3D
+
 ! field_chksum wrappers
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is

config_src/infra/FMS1/MOM_io_infra.F90

@@ -57,7 +57,7 @@ module MOM_io_infra
 !> Read a data field from a file
 interface read_field
   module procedure read_field_4d
-  module procedure read_field_3d
+  module procedure read_field_3d, read_field_3d_region
   module procedure read_field_2d, read_field_2d_region
   module procedure read_field_1d, read_field_1d_int
   module procedure read_field_0d, read_field_0d_int
@@ -696,6 +696,45 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, &
   endif ; endif
 end subroutine read_field_3d
 
+!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
+!! global 3-D data field named "fieldname" from file "filename".
+subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
+                                no_domain, scale)
+  character(len=*),       intent(in)    :: filename  !< The name of the file to read
+  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
+  real, dimension(:,:,:),   intent(inout) :: data    !< The 3-dimensional array into which the data
+                                                     !! should be read
+  integer, dimension(:),  intent(in)    :: start     !< The starting index to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  integer, dimension(:),  intent(in)    :: nread     !< The number of points to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  type(MOM_domain_type), &
+                optional, intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
+  logical,      optional, intent(in)    :: no_domain !< If present and true, this variable does not
+                                                     !! use domain decomposion.
+  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
+                                                     !! by before it is returned.
+
+  if (present(MOM_Domain)) then
+    call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, &
+                   no_domain=no_domain)
+  else
+    call read_data(filename, fieldname, data, start, nread, no_domain=no_domain)
+  endif
+
+  if (present(scale)) then ; if (scale /= 1.0) then
+    if (present(MOM_Domain)) then
+      call rescale_comp_data(MOM_Domain, data, scale)
+    else
+      ! Dangerously rescale the whole array
+      data(:,:,:) = scale*data(:,:,:)
+    endif
+  endif ; endif
+end subroutine read_field_3d_region
+
+
 !> This routine uses the fms_io subroutine read_data to read a distributed
 !! 4-D data field named "fieldname" from file "filename".  Valid values for
 !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.

config_src/infra/FMS2/MOM_coms_infra.F90

@@ -25,7 +25,7 @@ module MOM_coms_infra
 !> Communicate an array, string or scalar from one PE to others
 interface broadcast
   module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D
-  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D
+  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D
 end interface broadcast
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is
@@ -260,6 +260,27 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking)
 
 end subroutine broadcast_real2D
 
+!> Communicate a 3-D array of reals from one PE to others
+subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking)
+  real, dimension(:,:,:), intent(inout) :: dat       !< The data to communicate and destination
+  integer,              intent(in)    :: length    !< The total number of data elements
+  integer,    optional, intent(in)    :: from_PE   !< The source PE, by default the root PE
+  integer,    optional, intent(in)    :: PElist(:) !< The list of participating PEs, by default the
+                                                   !! active PE set as previously set via Set_PElist.
+  logical,    optional, intent(in)    :: blocking  !< If true, barriers are added around the call
+
+  integer :: src_PE   ! The processor that is sending the data
+  logical :: do_block ! If true add synchronizing barriers
+
+  do_block = .false. ; if (present(blocking)) do_block = blocking
+  if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif
+
+  if (do_block) call mpp_sync(PElist)
+  call mpp_broadcast(dat, length, src_PE, PElist)
+  if (do_block) call mpp_sync_self(PElist)
+
+end subroutine broadcast_real3D
+
 ! field_chksum wrappers
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is

config_src/infra/FMS2/MOM_io_infra.F90

@@ -66,7 +66,7 @@ module MOM_io_infra
 !> Read a data field from a file
 interface read_field
   module procedure read_field_4d
-  module procedure read_field_3d
+  module procedure read_field_3d, read_field_3d_region
   module procedure read_field_2d, read_field_2d_region
   module procedure read_field_1d, read_field_1d_int
   module procedure read_field_0d, read_field_0d_int
@@ -1030,6 +1030,74 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, &
 
 end subroutine read_field_3d
 
+!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
+!! global 3-D data field named "fieldname" from file "filename".
+subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
+                                no_domain, scale)
+  character(len=*),       intent(in)    :: filename  !< The name of the file to read
+  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
+  real, dimension(:,:,:),   intent(inout) :: data      !< The 3-dimensional array into which the data
+                                                     !! should be read
+  integer, dimension(:),  intent(in)    :: start     !< The starting index to read in each of 3
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th value is always 1.
+  integer, dimension(:),  intent(in)    :: nread     !< The number of points to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  type(MOM_domain_type), &
+                optional, intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
+  logical,      optional, intent(in)    :: no_domain !< If present and true, this variable does not
+                                                     !! use domain decomposion.
+  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
+                                                     !! by before it is returned.
+
+  ! Local variables
+  type(FmsNetcdfFile_t)       :: fileObj ! A handle to a non-domain-decomposed file
+  type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
+  character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
+  logical :: success               ! True if the file was successfully opened
+
+  if (present(MOM_Domain)) then
+    ! Open the FMS2 file-set.
+    success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
+    if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
+
+    ! Find the matching case-insensitive variable name in the file and prepare to read it.
+    call prepare_to_read_var(fileobj_DD, fieldname, "read_field_2d_region: ", &
+                             filename, var_to_read)
+
+    ! Read the data.
+    call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
+
+    ! Close the file-set.
+    if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
+  else
+    ! Open the FMS2 file-set.
+    success = fms2_open_file(fileObj, trim(filename), "read")
+    if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
+
+    ! Find the matching case-insensitive variable name in the file, and determine whether it
+    ! has a time dimension.
+    call find_varname_in_file(fileObj, fieldname, "read_field_2d_region: ", filename, var_to_read)
+
+    ! Read the data.
+    call fms2_read_data(fileobj, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
+
+    ! Close the file-set.
+    if (check_if_open(fileobj)) call fms2_close_file(fileobj)
+  endif
+
+  if (present(scale)) then ; if (scale /= 1.0) then
+    if (present(MOM_Domain)) then
+      call rescale_comp_data(MOM_Domain, data, scale)
+    else
+      ! Dangerously rescale the whole array
+      data(:,:,:) = scale*data(:,:,:)
+    endif
+  endif ; endif
+
+end subroutine read_field_3d_region
+
 !> This routine uses the fms_io subroutine read_data to read a distributed
 !! 4-D data field named "fieldname" from file "filename".  Valid values for
 !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.

src/core/MOM_PressureForce_FV.F90

@@ -519,7 +519,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS
 
   h_neglect = GV%H_subroundoff
-  dz_neglect = GV%H_subroundoff * GV%H_to_Z
+  dz_neglect = GV%dZ_subroundoff
   I_Rho0 = 1.0 / GV%Rho0
   G_Rho0 = GV%g_Earth / GV%Rho0
   rho_ref = CS%Rho0
@@ -827,12 +827,11 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
 
   mdl = "MOM_PressureForce_FV"
   call log_version(param_file, mdl, version, "")
-  call get_param(param_file, mdl, "RHO_0", CS%Rho0, &
-                 "The mean ocean density used with BOUSSINESQ true to "//&
-                 "calculate accelerations and the mass for conservation "//&
-                 "properties, or with BOUSSINSEQ false to convert some "//&
-                 "parameters from vertical units of m to kg m-2.", &
-                 units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
+  call get_param(param_file, mdl, "RHO_PGF_REF", CS%Rho0, &
+                 "The reference density that is subtracted off when calculating pressure "//&
+                 "gradient forces.  Its inverse is subtracted off of specific volumes when "//&
+                 "in non-Boussinesq mode.  The default is RHO_0.", &
+                 units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R)
   call get_param(param_file, mdl, "TIDES", CS%tides, &
                  "If true, apply tidal momentum forcing.", default=.false.)
   call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &

src/diagnostics/MOM_sum_output.F90

@@ -8,7 +8,7 @@ module MOM_sum_output
 use MOM_coms,          only : sum_across_PEs, PE_here, root_PE, num_PEs, max_across_PEs, field_chksum
 use MOM_coms,          only : reproducing_sum, reproducing_sum_EFP, EFP_to_real, real_to_EFP
 use MOM_coms,          only : EFP_type, operator(+), operator(-), assignment(=), EFP_sum_across_PEs
-use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, MOM_mesg
+use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe
 use MOM_file_parser,   only : get_param, log_param, log_version, param_file_type
 use MOM_forcing_type,  only : forcing
 use MOM_grid,          only : ocean_grid_type
@@ -1077,7 +1077,7 @@ subroutine depth_list_setup(G, GV, US, DL, CS)
         valid_DL_read = .true. ! Otherwise there would have been a fatal error.
       endif
     else
-      if (is_root_pe()) call MOM_error(WARNING, "depth_list_setup: "// &
+      if (is_root_pe()) call MOM_error(NOTE, "depth_list_setup: "// &
         trim(CS%depth_list_file)//" does not exist.  Creating a new file.")
       valid_DL_read = .false.
     endif

src/framework/MOM_horizontal_regridding.F90

@@ -309,6 +309,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
   real, dimension(:,:),  allocatable   :: tr_in      !< A 2-d array for holding input data on its
                                                      !! native horizontal grid, with units that change
                                                      !! as the input data is interpreted [a] then [A ~> a]
+  real, dimension(:,:,:),  allocatable   :: tr_in_full  !< A 3-d array for holding input data on the
+                                                     !! model horizontal grid, with units that change
+                                                     !! as the input data is interpreted [a] then [A ~> a]
   real, dimension(:,:),  allocatable   :: tr_inp     !< Native horizontal grid data extended to the poles
                                                      !! with units that change as the input data is
                                                      !! interpreted [a] then [A ~> a]
@@ -448,6 +451,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   if (is_ongrid) then
     allocate(tr_in(is:ie,js:je), source=0.0)
+    allocate(tr_in_full(is:ie,js:je,kd), source=0.0)
     allocate(mask_in(is:ie,js:je), source=0.0)
   else
     call horizontal_interp_init()
@@ -470,14 +474,19 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   ! Loop through each data level and interpolate to model grid.
   ! After interpolating, fill in points which will be needed to define the layers.
+
+  if (is_ongrid) then
+    start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = 1
+    count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = kd ; start(4) = 1 ; count(4) = 1
+    call MOM_read_data(trim(filename), trim(varnam), tr_in_full, start, count, G%Domain)
+  endif
+
   do k=1,kd
     mask_in(:,:)  = 0.0
     tr_out(:,:) = 0.0
 
     if (is_ongrid) then
-      start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = k
-      count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = 1 ; start(4) = 1 ; count(4) = 1
-      call MOM_read_data(trim(filename), trim(varnam), tr_in, start, count, G%Domain)
+      tr_in(is:ie,js:je) = tr_in_full(is:ie,js:je,k)
       do j=js,je
         do i=is,ie
           if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
@@ -594,7 +603,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   enddo ! kd
 
-  deallocate(lon_in, lat_in)
+  if (allocated(lat_inp)) deallocate(lat_inp)
+  deallocate(tr_in)
+  if (allocated(tr_inp)) deallocate(tr_inp)
+  if (allocated(tr_in_full)) deallocate(tr_in_full)
 
 end subroutine horiz_interp_and_extrap_tracer_record