Merge pull request #578 from NCAR/dev/ncar

This rolls dev/master forward. MOM6-examples will follow shortly...

Remember, to sync dev/master:
```
git checkout dev/master
git pull https://github.com/NOAA-GFDL/MOM6.git dev/master
```

config_src/coupled_driver/ocean_model_MOM.F90

@@ -89,7 +89,9 @@ module ocean_model_mod
 public ocean_model_restart
 public ice_ocn_bnd_type_chksum
 public ocean_public_type_chksum
-public    ocean_model_data_get
+public ocean_model_data_get
+public get_state_pointers
+
 interface ocean_model_data_get
   module procedure ocean_model_data1D_get
   module procedure ocean_model_data2D_get
@@ -1086,4 +1088,15 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)
 end subroutine ocean_public_type_chksum
 
+!> Returns pointers to objects within ocean_state_type
+subroutine get_state_pointers(OS, grid, surf)
+  type(ocean_state_type),          pointer :: OS !< Ocean state type
+  type(ocean_grid_type), optional, pointer :: grid !< Ocean grid
+  type(surface), optional, pointer         :: surf !< Ocean surface state
+
+  if (present(grid)) grid => OS%grid
+  if (present(surf)) surf=> OS%state
+
+end subroutine get_state_pointers
+
 end module ocean_model_mod

src/core/MOM_grid.F90

@@ -5,6 +5,7 @@ module MOM_grid
 
 use MOM_hor_index, only : hor_index_type, hor_index_init
 use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent
+use MOM_domains, only : get_global_shape
 use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL
 use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
 
@@ -13,7 +14,7 @@ module MOM_grid
 #include <MOM_memory.h>
 
 public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction
-public isPointInCell, hor_index_type
+public isPointInCell, hor_index_type, get_global_grid_size
 
 !> Ocean grid type. See mom_grid for details.
 type, public :: ocean_grid_type
@@ -443,6 +444,16 @@ subroutine set_first_direction(G, y_first)
   G%first_direction = y_first
 end subroutine set_first_direction
 
+!> Return global shape of horizontal grid
+subroutine get_global_grid_size(G, niglobal, njglobal)
+  type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type
+  integer,               intent(out)   :: niglobal !< i-index global size of grid
+  integer,               intent(out)   :: njglobal !< j-index global size of grid
+
+  call get_global_shape(G%domain, niglobal, njglobal)
+
+end subroutine get_global_grid_size
+
 !> Allocate memory used by the ocean_grid_type and related structures.
 subroutine allocate_metrics(G)
   type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type

src/framework/MOM_domains.F90

@@ -60,7 +60,7 @@ module MOM_domains
 public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners
 public :: create_group_pass, do_group_pass, group_pass_type
 public :: start_group_pass, complete_group_pass
-public :: compute_block_extent
+public :: compute_block_extent, get_global_shape
 
 interface pass_var
   module procedure pass_var_3d, pass_var_2d
@@ -1896,4 +1896,15 @@ subroutine get_domain_extent(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed, &
 
 end subroutine get_domain_extent
 
+!> Returns the global shape of h-point arrays
+subroutine get_global_shape(domain, niglobal, njglobal)
+  type(MOM_domain_type), intent(in)  :: domain !< MOM domain
+  integer,               intent(out) :: niglobal !< i-index global size of h-point arrays
+  integer,               intent(out) :: njglobal !< j-index global size of h-point arrays
+
+  niglobal = domain%niglobal
+  njglobal = domain%njglobal
+
+end subroutine get_global_shape
+
 end module MOM_domains

src/ice_shelf/MOM_ice_shelf.F90

@@ -1053,17 +1053,6 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
       endif
 
 
-      ! Add frazil formation diagnosed by the ocean model (J m-2) in the
-      ! form of surface layer evaporation (kg m-2 s-1). Update lprec in the
-      ! control structure for diagnostic purposes.
-
-      if (associated(state%frazil)) then
-        fraz = state%frazil(i,j) / CS%time_step / CS%Lat_fusion
-        if (associated(fluxes%evap)) fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
-        CS%lprec(i,j)=CS%lprec(i,j) - fraz
-        state%frazil(i,j) = 0.0
-      endif
-
       if (associated(fluxes%sens)) fluxes%sens(i,j) = -frac_area*CS%t_flux(i,j)*CS%flux_factor
       if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = frac_area * CS%salt_flux(i,j)*CS%flux_factor
       if (associated(fluxes%p_surf)) fluxes%p_surf(i,j) = frac_area * CS%g_Earth * CS%mass_shelf(i,j)

src/initialization/MOM_coord_initialization.F90

@@ -58,6 +58,7 @@ subroutine MOM_initialize_coord(GV, PF, write_geom, output_dir, tv, max_depth)
 ! Set-up the layer densities, GV%Rlay, and reduced gravities, GV%g_prime.
   call get_param(PF, mdl, "COORD_CONFIG", config, &
                  "This specifies how layers are to be defined: \n"//&
+                 " \t ALE or none - used to avoid defining layers in ALE mode \n"//&
                  " \t file - read coordinate information from the file \n"//&
                  " \t\t specified by (COORD_FILE).\n"//&
                  " \t BFB - Custom coords for buoyancy-forced basin case \n"//&
@@ -93,7 +94,8 @@ subroutine MOM_initialize_coord(GV, PF, write_geom, output_dir, tv, max_depth)
       call user_set_coord(GV%Rlay, GV%g_prime, GV, PF, eos)
     case ("BFB")
       call BFB_set_coord(GV%Rlay, GV%g_prime, GV, PF, eos)
-    case ("none")
+    case ("none", "ALE")
+      call set_coord_to_none(GV%Rlay, GV%g_prime, GV, PF)
     case default ; call MOM_error(FATAL,"MOM_initialize_coord: "// &
       "Unrecognized coordinate setup"//trim(config))
   end select
@@ -493,7 +495,38 @@ subroutine set_coord_linear(Rlay, g_prime, GV, param_file)
   call callTree_leave(trim(mdl)//'()')
 end subroutine set_coord_linear
 
-! -----------------------------------------------------------------------------
+!> Sets Rlay to Rho0 and g_prime to zero except for the free surface.
+!! This is for use only in ALE mode where Rlay should not be used and g_prime(1) alone
+!! might be used.
+subroutine set_coord_to_none(Rlay, g_prime, GV, param_file)
+  real, dimension(:),      intent(out) :: Rlay       !< The layers' target coordinate values
+                                                     !! (potential density).
+  real, dimension(:),      intent(out) :: g_prime    !< A structure indicating the open file to
+                                                     !! parse for model parameter values.
+  type(verticalGrid_type), intent(in)  :: GV         !< The ocean's vertical grid structure.
+  type(param_file_type),   intent(in)  :: param_file !< A structure to parse for run-time parameters
+  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
+  character(len=40)  :: mdl = "set_coord_to_none" ! This subroutine's name.
+  integer :: k, nz
+  nz = GV%ke
+
+  call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
+
+  call get_param(param_file, mdl, "GFS" , g_fs, &
+                 "The reduced gravity at the free surface.", units="m s-2", &
+                 default=GV%g_Earth)
+
+  g_prime(1) = g_fs
+  do k=2,nz ; g_prime(k) = 0. ; enddo
+  Rlay(1) = GV%Rho0
+  do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
+
+  call callTree_leave(trim(mdl)//'()')
+
+end subroutine set_coord_to_none
+
+!>   This subroutine writes out a file containing any available data related
+!! to the vertical grid used by the MOM ocean model.
 subroutine write_vertgrid_file(GV, param_file, directory)
   type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
   type(param_file_type), intent(in)    :: param_file !< A structure to parse for run-time parameters

src/initialization/MOM_state_initialization.F90

@@ -1795,10 +1795,14 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, just_read_params)
   logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                       !! only read parameters without changing h.
 
-  character(len=200) :: filename   ! The name of an input file containing temperature
-                                   ! and salinity in z-space; also used for  ice shelf area.
-  character(len=200) :: shelf_file ! The name of an input file used for  ice shelf area.
-  character(len=200) :: inputdir ! The directory where NetCDF input filesare.
+  character(len=200) :: filename   !< The name of an input file containing temperature
+                                   !! and salinity in z-space; also used for  ice shelf area.
+  character(len=200) :: tfilename  !< The name of an input file containing only temperature
+                                   !! in z-space.
+  character(len=200) :: sfilename  !< The name of an input file containing only salinity
+                                   !! in z-space.
+  character(len=200) :: shelf_file !< The name of an input file used for  ice shelf area.
+  character(len=200) :: inputdir   !! The directory where NetCDF input filesare.
   character(len=200) :: mesg, area_varname, ice_shelf_file
 
   type(EOS_type), pointer :: eos => NULL()
@@ -1892,15 +1896,24 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, just_read_params)
 
   call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE",filename, &
                  "The name of the z-space input file used to initialize \n"//&
-                 "the layer thicknesses, temperatures and salinities.", &
-                 default="temp_salt_z.nc", do_not_log=just_read)
+                 "temperatures (T) and salinities (S). If T and S are not \n" //&
+                 "in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE \n" //&
+                 "must be set.",default="temp_salt_z.nc",do_not_log=just_read)
+  call get_param(PF, mdl, "TEMP_Z_INIT_FILE",tfilename, &
+                 "The name of the z-space input file used to initialize \n"//&
+                 "temperatures, only.", default=trim(filename),do_not_log=just_read)
+  call get_param(PF, mdl, "SALT_Z_INIT_FILE",sfilename, &
+                 "The name of the z-space input file used to initialize \n"//&
+                 "temperatures, only.", default=trim(filename),do_not_log=just_read)
   filename = trim(inputdir)//trim(filename)
+  tfilename = trim(inputdir)//trim(tfilename)
+  sfilename = trim(inputdir)//trim(sfilename)
   call get_param(PF, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
                  "The name of the potential temperature variable in \n"//&
-                 "TEMP_SALT_Z_INIT_FILE.", default="ptemp", do_not_log=just_read)
+                 "TEMP_Z_INIT_FILE.", default="ptemp",do_not_log=just_read)
   call get_param(PF, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
                  "The name of the salinity variable in \n"//&
-                 "TEMP_SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read)
+                 "SALT_Z_INIT_FILE.", default="salt",do_not_log=just_read)
   call get_param(PF, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
                  "If True, then horizontally homogenize the interpolated \n"//&
                  "initial conditions.", default=.false., do_not_log=just_read)
@@ -1965,10 +1978,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, just_read_params)
 !   to the North/South Pole past the limits of the input data, they are extrapolated using the average
 !   value at the northernmost/southernmost latitude.
 
-  call horiz_interp_and_extrap_tracer(filename, potemp_var,1.0,1, &
+  call horiz_interp_and_extrap_tracer(tfilename, potemp_var,1.0,1, &
        G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, tripolar_n, homogenize)
 
-  call horiz_interp_and_extrap_tracer(filename, salin_var,1.0,1, &
+  call horiz_interp_and_extrap_tracer(sfilename, salin_var,1.0,1, &
        G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, tripolar_n, homogenize)
 
   kd = size(z_in,1)

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -580,7 +580,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
         v0(i,J) = CS%IDXDY2v(i,J)*(CS%DY2q(I,J)*sh_xy(I,J) - CS%DY2q(I-1,J)*sh_xy(I-1,J)) - &
                   CS%IDX2dyCv(i,J)*(CS%DX2h(i,j+1)*sh_xx(i,j+1) - CS%DX2h(i,j)*sh_xx(i,j))
       enddo ; enddo
-      if (apply_OBC .and. OBC%zero_biharmonic) then
+      if (apply_OBC) then; if (OBC%zero_biharmonic) then
         do n=1,OBC%number_of_segments
           I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
           if (OBC%segment(n)%is_N_or_S .and. (J >= Jsq-1) .and. (J <= Jeq+1)) then
@@ -593,7 +593,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
             enddo
           endif
         enddo
-      endif
+      endif; endif
     endif
 
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1