Merge branch 'dev/gfdl' into add_diagnostic

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -27,7 +27,8 @@ module MOM_surface_forcing_gfdl
 use MOM_grid,             only : ocean_grid_type
 use MOM_interpolate,      only : init_external_field, time_interp_external
 use MOM_interpolate,      only : time_interp_external_init
-use MOM_io,               only : slasher, write_version_number, MOM_read_data, stdout
+use MOM_io,               only : slasher, write_version_number, MOM_read_data
+use MOM_io,               only : stdout_if_root
 use MOM_restart,          only : register_restart_field, restart_init, MOM_restart_CS
 use MOM_restart,          only : restart_init_end, save_restart, restore_state
 use MOM_string_functions, only : uppercase
@@ -1628,8 +1629,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   logical :: root    ! True only on the root PE
   integer :: outunit ! The output unit to write to
 
-  outunit = stdout
   root = is_root_pe()
+  outunit = stdout_if_root()
 
   if (root) write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
   chks = field_chksum( iobt%u_flux         ) ; if (root) write(outunit,100) 'iobt%u_flux         ', chks

config_src/drivers/FMS_cap/ocean_model_MOM.F90

@@ -35,7 +35,7 @@ module ocean_model_mod
 use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags
 use MOM_get_input, only : Get_MOM_Input, directories
 use MOM_grid, only : ocean_grid_type
-use MOM_io, only : write_version_number, stdout
+use MOM_io, only : write_version_number, stdout_if_root
 use MOM_marine_ice, only : iceberg_forces, iceberg_fluxes, marine_ice_init, marine_ice_CS
 use MOM_restart, only : MOM_restart_CS, save_restart
 use MOM_string_functions, only : uppercase
@@ -1107,8 +1107,8 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
   logical :: root    ! True only on the root PE
   integer :: outunit ! The output unit to write to
 
-  outunit = stdout
   root = is_root_pe()
+  outunit = stdout_if_root()
 
   if (root) write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep
   chks = field_chksum(ocn%t_surf ) ; if (root) write(outunit,100) 'ocean%t_surf   ', chks

config_src/infra/FMS1/MOM_io_infra.F90

@@ -17,6 +17,7 @@ module MOM_io_infra
 use mpp_io_mod,           only : mpp_get_fields, fieldtype
 use mpp_io_mod,           only : mpp_get_info, mpp_get_times
 use mpp_io_mod,           only : mpp_io_init
+use mpp_mod,              only : stdout_if_root=>stdout
 ! These are encoding constants.
 use mpp_io_mod,           only : APPEND_FILE=>MPP_APPEND, WRITEONLY_FILE=>MPP_WRONLY
 use mpp_io_mod,           only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY
@@ -33,6 +34,7 @@ module MOM_io_infra
 public :: MOM_read_data, MOM_read_vector, write_metadata, write_field
 public :: field_exists, get_field_atts, get_field_size, get_axis_data, read_field_chksum
 public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version
+public :: stdout_if_root
 ! These types are inherited from underlying infrastructure code, to act as containers for
 ! information about fields and axes, respectively, and are opaque to this module.
 public :: fieldtype, axistype

config_src/infra/FMS2/MOM_io_infra.F90

@@ -29,6 +29,7 @@ module MOM_io_infra
 use mpp_io_mod,           only : mpp_get_fields, mpp_fieldtype=>fieldtype
 use mpp_io_mod,           only : mpp_get_info, mpp_get_times
 use mpp_io_mod,           only : mpp_io_init
+use mpp_mod,              only : stdout_if_root=>stdout
 ! These are encoding constants.
 use mpp_io_mod,           only : APPEND_FILE=>MPP_APPEND, WRITEONLY_FILE=>MPP_WRONLY
 use mpp_io_mod,           only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY
@@ -44,6 +45,7 @@ module MOM_io_infra
 public :: MOM_read_data, MOM_read_vector, write_metadata, write_field
 public :: field_exists, get_field_atts, get_field_size, get_axis_data, read_field_chksum
 public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version
+public :: stdout_if_root
 ! These types act as containers for information about files, fields and axes, respectively,
 ! and may also wrap opaque types from the underlying infrastructure.
 public :: file_type, fieldtype, axistype

src/framework/MOM_diag_mediator.F90

@@ -1051,6 +1051,7 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num
   else
     axes%v_cell_method = ''
   endif
+
   if (present(nz)) axes%nz = nz
   if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number
   if (present(is_h_point)) axes%is_h_point = is_h_point
@@ -1971,38 +1972,45 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
   type(diag_ctrl), pointer :: diag_cs => NULL()
   type(axes_grp), pointer :: remap_axes => null()
   type(axes_grp), pointer :: axes => null()
+  type(axes_grp), pointer :: axes_d2 => null()
   integer :: dm_id, i, dl
   character(len=256) :: msg, cm_string
   character(len=256) :: new_module_name
   character(len=480) :: module_list, var_list
   integer :: num_modnm, num_varnm
   logical :: active
 
-  axes => axes_in
-  MOM_missing_value = axes%diag_cs%missing_value
-  if (present(missing_value)) MOM_missing_value = missing_value
-
-  diag_cs => axes%diag_cs
-  dm_id = -1
+  diag_cs => axes_in%diag_cs
 
+  ! Check if the axes match a standard grid axis.
+  ! If not, allocate the new axis and copy the contents.
   if (axes_in%id == diag_cs%axesTL%id) then
     axes => diag_cs%axesTL
   elseif (axes_in%id == diag_cs%axesBL%id) then
     axes => diag_cs%axesBL
-  elseif (axes_in%id == diag_cs%axesCuL%id ) then
+  elseif (axes_in%id == diag_cs%axesCuL%id) then
     axes => diag_cs%axesCuL
   elseif (axes_in%id == diag_cs%axesCvL%id) then
     axes => diag_cs%axesCvL
   elseif (axes_in%id == diag_cs%axesTi%id) then
     axes => diag_cs%axesTi
   elseif (axes_in%id == diag_cs%axesBi%id) then
     axes => diag_cs%axesBi
-  elseif (axes_in%id == diag_cs%axesCui%id ) then
+  elseif (axes_in%id == diag_cs%axesCui%id) then
     axes => diag_cs%axesCui
   elseif (axes_in%id == diag_cs%axesCvi%id) then
     axes => diag_cs%axesCvi
+  else
+    allocate(axes)
+    axes = axes_in
   endif
 
+  MOM_missing_value = axes%diag_cs%missing_value
+  if (present(missing_value)) MOM_missing_value = missing_value
+
+  diag_cs => axes%diag_cs
+  dm_id = -1
+
   module_list = "{"//trim(module_name)
   num_modnm = 1
 
@@ -2090,40 +2098,40 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
     new_module_name = trim(module_name)//'_d2'
 
     if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then
-      axes => null()
+      axes_d2 => null()
       if (axes_in%id == diag_cs%axesTL%id) then
-        axes => diag_cs%dsamp(dl)%axesTL
+        axes_d2 => diag_cs%dsamp(dl)%axesTL
       elseif (axes_in%id == diag_cs%axesBL%id) then
-        axes => diag_cs%dsamp(dl)%axesBL
+        axes_d2 => diag_cs%dsamp(dl)%axesBL
       elseif (axes_in%id == diag_cs%axesCuL%id ) then
-        axes => diag_cs%dsamp(dl)%axesCuL
+        axes_d2 => diag_cs%dsamp(dl)%axesCuL
       elseif (axes_in%id == diag_cs%axesCvL%id) then
-        axes => diag_cs%dsamp(dl)%axesCvL
+        axes_d2 => diag_cs%dsamp(dl)%axesCvL
       elseif (axes_in%id == diag_cs%axesTi%id) then
-        axes => diag_cs%dsamp(dl)%axesTi
+        axes_d2 => diag_cs%dsamp(dl)%axesTi
       elseif (axes_in%id == diag_cs%axesBi%id) then
-        axes => diag_cs%dsamp(dl)%axesBi
+        axes_d2 => diag_cs%dsamp(dl)%axesBi
       elseif (axes_in%id == diag_cs%axesCui%id ) then
-        axes => diag_cs%dsamp(dl)%axesCui
+        axes_d2 => diag_cs%dsamp(dl)%axesCui
       elseif (axes_in%id == diag_cs%axesCvi%id) then
-        axes => diag_cs%dsamp(dl)%axesCvi
+        axes_d2 => diag_cs%dsamp(dl)%axesCvi
       elseif (axes_in%id == diag_cs%axesT1%id) then
-        axes => diag_cs%dsamp(dl)%axesT1
+        axes_d2 => diag_cs%dsamp(dl)%axesT1
       elseif (axes_in%id == diag_cs%axesB1%id) then
-        axes => diag_cs%dsamp(dl)%axesB1
+        axes_d2 => diag_cs%dsamp(dl)%axesB1
       elseif (axes_in%id == diag_cs%axesCu1%id ) then
-        axes => diag_cs%dsamp(dl)%axesCu1
+        axes_d2 => diag_cs%dsamp(dl)%axesCu1
       elseif (axes_in%id == diag_cs%axesCv1%id) then
-        axes => diag_cs%dsamp(dl)%axesCv1
+        axes_d2 => diag_cs%dsamp(dl)%axesCv1
       else
         !Niki: Should we worry about these, e.g., diag_to_Z_CS?
         call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " &
               //trim(new_module_name)//"-"//trim(field_name))
       endif
     endif
     ! Register the native diagnostic
-    if (associated(axes)) then
-       active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, &
+    if (associated(axes_d2)) then
+       active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, &
                 init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, &
                 range=range, mask_variant=mask_variant, standard_name=standard_name, &
                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
@@ -2196,7 +2204,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
                              v_extensive=v_extensive)
     module_list = trim(module_list)//"}"
     if (num_modnm <= 1) module_list = module_name
-    if (num_varnm <= 1) var_list = ""
+    if (num_varnm <= 1) var_list = ''
 
     call log_available_diag(dm_id>0, module_list, field_name, cm_string, msg, diag_CS, &
                             long_name, units, standard_name, variants=var_list)
@@ -2216,7 +2224,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name,
   integer,          intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group
   character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model"
   character(len=*), intent(in) :: field_name !< Name of the diagnostic field
-  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes
+  type(axes_grp),   intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes
                                              !! for this field
   type(time_type),  intent(in) :: init_time !< Time at which a field is first available?
   character(len=*), optional, intent(in) :: long_name !< Long name of a field.

src/framework/MOM_horizontal_regridding.F90

@@ -259,7 +259,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
 
   character(len=*),      intent(in)    :: filename   !< Path to file containing tracer to be
                                                      !! interpolated.
-  character(len=*),      intent(in)    :: varnam     !< Name of tracer in filee.
+  character(len=*),      intent(in)    :: varnam     !< Name of tracer in file.
   real,                  intent(in)    :: conversion !< Conversion factor for tracer.
   integer,               intent(in)    :: recnum     !< Record number of tracer to be read.
   type(ocean_grid_type), intent(inout) :: G          !< Grid object
@@ -348,9 +348,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
                                  " in file "//trim(filename)//" in hinterp_extrap")
 
   rcode = NF90_INQUIRE_VARIABLE(ncid, varid, ndims=ndims, dimids=dims)
-  if (rcode /= 0) call MOM_error(FATAL,'error inquiring dimensions hinterp_extrap')
-  if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(varnam)//" in file "// &
-              trim(filename)//" has too few dimensions.")
+  if (rcode /= 0) call MOM_error(FATAL, "Error inquiring about the dimensions of "//trim(varnam)//&
+                                 " in file "//trim(filename)//" in hinterp_extrap")
+  if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(varnam)//" in file "//trim(filename)// &
+                                " has too few dimensions to be read as a 3-d array.")
 
   rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id)
   if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 data for "// &
@@ -373,8 +374,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
 
   missing_value=0.0
   rcode = NF90_GET_ATT(ncid, varid, "_FillValue", missing_value)
-  if (rcode /= 0) call MOM_error(FATAL,"error finding missing value for "//&
-       trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
+  if (rcode /= 0) call MOM_error(FATAL,"error finding missing value for "//trim(varnam)//&
+                                 " in file "// trim(filename)//" in hinterp_extrap")
 
   rcode = NF90_GET_ATT(ncid, varid, "add_offset", add_offset)
   if (rcode /= 0) add_offset = 0.0
@@ -465,7 +466,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
       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
       rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count)
-      if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//&
+      if (rcode /= 0) call MOM_error(FATAL,"horiz_interp_and_extrap_tracer_record: "//&
            "error reading level "//trim(laynum)//" of variable "//&
            trim(varnam)//" in file "// trim(filename))
 
@@ -484,7 +485,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
       if (is_root_pe()) then
         start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd
         rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count)
-        if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//&
+        if (rcode /= 0) call MOM_error(FATAL,"horiz_interp_and_extrap_tracer_record: "//&
              "error reading level "//trim(laynum)//" of variable "//&
              trim(varnam)//" in file "// trim(filename))
 

src/framework/MOM_io.F90

@@ -19,6 +19,7 @@ module MOM_io
 use MOM_io_infra,         only : get_file_times, axistype, get_axis_data, get_filename_suffix
 use MOM_io_infra,         only : write_field, write_metadata, write_version
 use MOM_io_infra,         only : MOM_namelist_file, check_namelist_error, io_infra_init, io_infra_end
+use MOM_io_infra,         only : stdout_if_root
 use MOM_io_infra,         only : APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE
 use MOM_io_infra,         only : READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE
 use MOM_io_infra,         only : CENTER, CORNER, NORTH_FACE, EAST_FACE
@@ -47,6 +48,7 @@ module MOM_io
 public :: MOM_read_data, MOM_read_vector, read_field_chksum
 public :: slasher, write_field, write_version_number
 public :: io_infra_init, io_infra_end
+public :: stdout_if_root
 ! This is used to set up information descibing non-domain-decomposed axes.
 public :: axis_info, set_axis_info, delete_axis_info
 ! This is used to set up global file attributes
@@ -647,15 +649,16 @@ subroutine get_var_sizes(filename, varname, ndims, sizes, match_case, caller, al
     do n=2,nval ;  sizes(n-1) = size_msg(n) ; enddo
     deallocate(size_msg)
 
-    if (present(dim_names)) then
+    if (present(dim_names) .and. (ndims > 0)) then
       nval = min(ndims, size(dim_names))
       call broadcast(dim_names(1:nval), len(dim_names(1)), blocking=.true.)
     endif
   endif
 
 end subroutine get_var_sizes
 
-!> read_var_sizes returns the number and size of dimensions associate with a variable in a file.
+!> read_var_sizes returns the number and size of dimensions associated with a variable in a file.
+!! If the variable is not in the file the returned sizes are all 0 and ndims is -1.
 !! Every processor for which this is called does the reading.
 subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, dim_names, ncid_in)
   character(len=*),      intent(in)  :: filename   !< Name of the file to read, used here in messages
@@ -675,7 +678,7 @@ subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, d
   character(len=256) :: hdr, dimname
   integer, allocatable :: dimids(:)
   integer :: varid, ncid, n, status
-  logical :: success
+  logical :: success, found
   hdr = "get_var_size: " ; if (present(caller)) hdr = trim(hdr)//": "
   sizes(:) = 0 ; ndims = -1
 
@@ -687,8 +690,8 @@ subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, d
   endif
 
   ! Get the dimension sizes of the variable varname.
-  call get_varid(varname, ncid, filename, varid, match_case=match_case)
-  if (varid < 0) return
+  call get_varid(varname, ncid, filename, varid, match_case=match_case, found=found)
+  if (.not.found) return
 
   status = NF90_inquire_variable(ncid, varid, ndims=ndims)
   if (status /= NF90_NOERR) then