diff --git a/.github/workflows/hola_tierra.yml b/.github/workflows/hola_tierra.yml deleted file mode 100644 index 82f9b7d87e..0000000000 --- a/.github/workflows/hola_tierra.yml +++ /dev/null @@ -1,27 +0,0 @@ -# This is a basic workflow to help you get started with Actions - -name: CI - -# Controls when the action will run. Triggers the workflow on push or pull request -# events but only for the dev/gfdl branch -on: - push: - branches: [ user/jml/add_fms2io_to_MOM_restart ] - pull_request: - branches: [ dev/gfdl ] - -# A workflow run is made up of one or more jobs that can run sequentially or in parallel -jobs: - # This workflow contains a single job called "build" - build: - # The type of runner that the job will run on - runs-on: ubuntu-latest - - # Steps represent a sequence of tasks that will be executed as part of the job - steps: - # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - - uses: actions/checkout@v2 - - # Runs a single command using the runners shell - - name: Run a one-line script - run: echo Hola, tierra! diff --git a/.testing/Makefile b/.testing/Makefile index ab978fdadc..05fb630a31 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -20,7 +20,7 @@ MKMF := $(abspath $(DEPS)/mkmf/bin/mkmf) # FMS framework FMS_URL ?= https://github.com/NOAA-GFDL/FMS.git -FMS_COMMIT ?= 2019.01.03 +FMS_COMMIT ?= 2020.03-alpha1 FMS := $(DEPS)/fms #--- diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index 7075fb7c10..4a730d6e6d 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -1224,7 +1224,8 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & if (.not.associated(CS)) return if (.not.associated(CS%restart_CSp)) return - call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface; omit this argument to use the old interface + call save_restart(directory, Time, G, CS%restart_CSp, time_stamped, use_fms2=.true.) end subroutine forcing_save_restart @@ -1589,8 +1590,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. if (.not.new_sim) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface; omit this argument to use the old interface call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & - G, CS%restart_CSp) + G, CS%restart_CSp, use_fms2=.true.) endif endif diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 082099158c..ff365a9e78 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -684,8 +684,9 @@ subroutine ocean_model_restart(OS, timestamp) "restart files can only be created after the buoyancy forcing is applied.") if (BTEST(OS%Restart_control,1)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, .true., GV=OS%GV) + OS%restart_CSp, .true., GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, .true.) if (OS%use_ice_shelf) then @@ -693,8 +694,9 @@ subroutine ocean_model_restart(OS, timestamp) endif endif if (BTEST(OS%Restart_control,0)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV) + OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then @@ -746,8 +748,8 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) if (present(directory)) then ; restart_dir = directory else ; restart_dir = OS%dirs%restart_output_dir ; endif - - call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index f8a4a19532..3c75cb12eb 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -690,8 +690,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) "restart files can only be created after the buoyancy forcing is applied.") if (present(restartname)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV, filename=restartname) + OS%restart_CSp, GV=OS%GV, filename=restartname, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) ! Is this needed? if (OS%use_ice_shelf) then @@ -700,8 +701,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) endif else if (BTEST(OS%Restart_control,1)) then + ! NOTE:use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, .true., GV=OS%GV) + OS%restart_CSp, .true., GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, .true.) if (OS%use_ice_shelf) then @@ -709,8 +711,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) endif endif if (BTEST(OS%Restart_control,0)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV) + OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then @@ -766,7 +769,7 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) if (present(directory)) then ; restart_dir = directory else ; restart_dir = OS%dirs%restart_output_dir ; endif - call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV) + call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index a42a8c3015..88b7f01654 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -1001,7 +1001,8 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & if (.not.associated(CS)) return if (.not.associated(CS%restart_CSp)) return - call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(.true., directory, Time, G, CS%restart_CSp, time_stamped, use_fms2=.true.) end subroutine forcing_save_restart @@ -1325,8 +1326,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. if (.not.new_sim) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & - G, CS%restart_CSp) + G, CS%restart_CSp, use_fms2=.true.) endif endif diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index b1ce9a60c0..9f1912d79f 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -534,7 +534,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') trim(runid), year, month, day, seconds call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & - glb%ocn_state%restart_CSp, .false., filename=restartname, GV=glb%ocn_state%GV) + glb%ocn_state%restart_CSp, .false., filename=restartname, GV=glb%ocn_state%GV, use_fms2=.true.) ! write name of restart file in the rpointer file nu = shr_file_getUnit() diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 9946aec4f9..a8765bdc08 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -686,8 +686,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) "restart files can only be created after the buoyancy forcing is applied.") if (present(restartname)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV, filename=restartname) + OS%restart_CSp, GV=OS%GV, filename=restartname, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) ! Is this needed? if (OS%use_ice_shelf) then @@ -696,8 +697,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) endif else if (BTEST(OS%Restart_control,1)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, .true., GV=OS%GV) + OS%restart_CSp, .true., GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, .true.) if (OS%use_ice_shelf) then @@ -705,8 +707,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname) endif endif if (BTEST(OS%Restart_control,0)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(OS%dirs%restart_output_dir, OS%Time, OS%grid, & - OS%restart_CSp, GV=OS%GV) + OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then @@ -760,8 +763,8 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) if (present(directory)) then ; restart_dir = directory else ; restart_dir = OS%dirs%restart_output_dir ; endif - - call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(restart_dir, Time, OS%grid, OS%restart_CSp, GV=OS%GV, use_fms2=.true.) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index 3d49c66ce6..a565da3d93 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -1000,7 +1000,8 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & if (.not.associated(CS)) return if (.not.associated(CS%restart_CSp)) return - call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(directory, Time, G, CS%restart_CSp, time_stamped, use_fms2=.true.) end subroutine forcing_save_restart @@ -1330,8 +1331,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. if (.not.new_sim) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & - G, CS%restart_CSp) + G, CS%restart_CSp, use_fms2=.true.) endif endif diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index f180cd9717..c6fbe0e4e6 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -583,16 +583,18 @@ program MOM_main if ((permit_incr_restart) .and. (fluxes%fluxes_used) .and. & (Time + (Time_step_ocean/2) > restart_time)) then if (BTEST(Restart_control,1)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(dirs%restart_output_dir, Time, grid, & - restart_CSp, .true., GV=GV) + restart_CSp, .true., GV=GV, use_fms2=.true.) call forcing_save_restart(surface_forcing_CSp, grid, Time, & dirs%restart_output_dir, .true.) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & dirs%restart_output_dir, .true.) endif if (BTEST(Restart_control,0)) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call save_restart(dirs%restart_output_dir, Time, grid, & - restart_CSp, GV=GV) + restart_CSp, GV=GV, use_fms2=.true.) call forcing_save_restart(surface_forcing_CSp, grid, Time, & dirs%restart_output_dir) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & @@ -616,8 +618,8 @@ program MOM_main "End of MOM_main reached with unused buoyancy fluxes. "//& "For conservation, the ocean restart files can only be "//& "created after the buoyancy forcing is applied.") - - call save_restart(dirs%restart_output_dir, Time, grid, restart_CSp, GV=GV) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(dirs%restart_output_dir, Time, grid, restart_CSp, GV=GV, use_fms2=.true.) if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & dirs%restart_output_dir) ! Write ocean solo restart file. diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 0a56abb681..5b10ea46e4 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -1524,8 +1524,8 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & if (.not.associated(CS)) return if (.not.associated(CS%restart_CSp)) return - - call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call save_restart(directory, Time, G, CS%restart_CSp, time_stamped, use_fms2=.true.) end subroutine forcing_save_restart @@ -1925,8 +1925,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. if (.not.new_sim) then + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & - G, CS%restart_CSp) + G, CS%restart_CSp, use_fms2=.true.) endif endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a9b9c7fec4..2c48796f57 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2782,9 +2782,9 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0) call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, & "Interface heights", "meter", z_grid='i') - + ! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface call save_restart(dirs%output_directory, Time, CS%G_in, & - restart_CSp_tmp, filename=CS%IC_file, GV=GV) + restart_CSp_tmp, filename=CS%IC_file, GV=GV, write_ic=.true.) deallocate(z_interface) deallocate(restart_CSp_tmp) endif diff --git a/src/framework/MOM_axis.F90 b/src/framework/MOM_axis.F90 new file mode 100644 index 0000000000..48f70bec70 --- /dev/null +++ b/src/framework/MOM_axis.F90 @@ -0,0 +1,625 @@ +!> This module contains routines that define and register axes to files +module MOM_axis + +! This file is part of MOM6. See LICENSE.md for the license. +use MOM_domains, only : MOM_domain_type +use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING +use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_string_functions, only : lowercase +use MOM_verticalGrid, only : verticalGrid_type +use fms2_io_mod, only : is_dimension_registered, register_axis, is_dimension_unlimited +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited +use fms2_io_mod, only : get_variable_size, get_variable_num_dimensions, check_if_open +use fms2_io_mod, only : fms2_open_file=>open_file, fms2_close_file=>close_file +use fms2_io_mod, only : get_variable_dimension_names, read_data, get_unlimited_dimension_name +use fms2_io_mod, only : get_dimension_size +use mpp_domains_mod, only : domain2d, CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST +use mpp_domains_mod, only : mpp_get_compute_domain +use netcdf +implicit none ; private + +public MOM_register_diagnostic_axis, get_var_dimension_metadata, get_time_units +public MOM_get_diagnostic_axis_data, MOM_register_variable_axes, get_time_index +public convert_checksum_to_string +!> A type for making arrays of pointers to real 1-d arrays +type p1d + real, dimension(:), pointer :: p => NULL() !< A pointer to a 1d array +end type p1d + +!> A structure with information about a single axis variable +type axis_atts + character(len=64) :: name !< Names of the axis + character(len=48) :: units !< Physical dimensions of the axis + character(len=240) :: longname !< Long name of the axis + character(len=8) :: positive !< Positive-definite direction: up, down, east, west, north, south + integer :: horgrid_position !< Horizontal grid position + logical :: is_domain_decomposed !< if .true. the axis data are domain-decomposed + !! and need to be indexed by the compute domain + !! before passing to write_data +end type axis_atts + +!> Type for describing an axis variable (e.g., lath, lonh, Time) +type, public :: axis_data_type + !> An array of descriptions of the registered axes + type(axis_atts), pointer :: axis(:) => NULL() !< structure with axis attributes + type(p1d), pointer :: data(:) => NULL() !< pointer to the axis data +end type axis_data_type + +!> interface for registering axes associated with a variable to a netCDF file object +interface MOM_register_variable_axes + module procedure MOM_register_variable_axes_subdomain + module procedure MOM_register_variable_axes_full +end interface MOM_register_variable_axes + +contains + +!> register a MOM diagnostic axis to a domain-decomposed file +subroutine MOM_register_diagnostic_axis(fileObj, axisName, axisLength) + type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file + character(len=*), intent(in) :: axisName !< name of the axis to register to file + integer, intent(in), optional :: axisLength !< length of axis/dimension ;only needed for Layer, Interface, Time, + !! Period + select case (trim(lowercase(axisName))) + case ('latq'); call register_axis(fileObj,'latq','y', domain_position=NORTH_FACE) + case ('lath'); call register_axis(fileObj,'lath','y', domain_position=CENTER) + case ('lonq'); call register_axis(fileObj,'lonq','x', domain_position=EAST_FACE) + case ('lonh'); call register_axis(fileObj,'lonh','x', domain_position=CENTER) + case default + if (.not. present(axisLength)) call MOM_error(FATAL,"MOM_io:register_diagnostic_axis: "//& + "An axis_length argument is required to register the axis "//trim(axisName)) + call register_axis(fileObj, trim(axisName), axisLength) + end select +end subroutine MOM_register_diagnostic_axis + + +!> Get the horizontal grid, vertical grid, and/or time dimension names and lengths +!! for a single variable from the hor_grid, t_grid, and z_grid values returned by a prior call to query_vardesc +subroutine get_var_dimension_metadata(hor_grid, z_grid, t_grid_in, & + dim_names, dim_lengths, num_dims, G, dG, GV) + + character(len=*), intent(in) :: hor_grid !< horizontal grid + character(len=*), intent(in) :: z_grid !< vertical grid + character(len=*), intent(in) :: t_grid_in !< time grid + character(len=*), dimension(:), intent(inout) :: dim_names !< array of dimension names + integer, dimension(:), intent(inout) :: dim_lengths !< array of dimension sizes + integer, intent(inout) :: num_dims !< number of axes to register in the restart file + type(ocean_grid_type), optional, intent(in) :: G !< The ocean's grid structure + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure + + ! local + logical :: use_lath + logical :: use_lonh + logical :: use_latq + logical :: use_lonq + character(len=8) :: t_grid + character(len=8) :: t_grid_read + integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB + !integer :: npes + real, pointer, dimension(:) :: gridLatT => NULL(), & ! The latitude or longitude of T or B points for + gridLatB => NULL(), & ! the purpose of labeling the output axes. + gridLonT => NULL(), & + gridLonB => NULL() + type(MOM_domain_type), pointer :: domain => NULL() ! Domain used to get the pe count + + use_lath = .false. + use_lonh = .false. + use_latq = .false. + use_lonq = .false. + + ! set the ocean grid coordinates + + if (present(G)) then + gridLatT => G%gridLatT ; gridLatB => G%gridLatB + gridLonT => G%gridLonT ; gridLonB => G%gridLonB + isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg + IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB + + call get_horizontal_grid_logic(hor_grid, use_lath, use_lonh, use_latq, use_lonq) + elseif (present(dG)) then + gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB + gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB + isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg + IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB + + call get_horizontal_grid_logic(hor_grid, use_lath, use_lonh, use_latq, use_lonq) + endif + + ! add longitude name to dimension name array + if (use_lonh) then + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("lonh")) = "lonh" + dim_lengths(num_dims) = size(gridLonT(isg:ieg)) + elseif (use_lonq) then + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("lonq")) = "lonq" + dim_lengths(num_dims) = size(gridLonB(IsgB:IegB)) + endif + ! add latitude name to dimension name array + if (use_lath) then + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("lath")) = "lath" + dim_lengths(num_dims) = size(gridLatT(jsg:jeg)) + elseif (use_latq) then + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("latq")) = "latq" + dim_lengths(num_dims) = size(gridLatB(JsgB:JegB)) + endif + + if (present(GV)) then + ! vertical grid + select case (trim(z_grid)) + case ('L') + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("Layer")) = "Layer" + dim_lengths(num_dims) = GV%ke + case ('i') + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("Interface")) = "Interface" + dim_lengths(num_dims) = GV%ke+1 + case ('1') ! Do nothing. + case default + call MOM_error(FATAL, "MOM_io: get_var_dimension_features: "//& + " has an unrecognized z_grid argument"//trim(z_grid)) + end select + endif + ! time + t_grid = adjustl(t_grid_in) + select case (t_grid(1:1)) + case ('s', 'a', 'm') + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("Time")) = "Time" + dim_lengths(num_dims) = unlimited + case ('p') + if (len_trim(t_grid(2:8)) <= 0) then + call MOM_error(FATAL,"MOM_io:get_var_dimension_features: "//& + "No periodic axis length was specified in "//trim(t_grid)) + endif + num_dims = num_dims+1 + dim_names(num_dims) = "" + dim_names(num_dims)(1:len_trim("Period")) = "Period" + dim_lengths(num_dims) = unlimited + case ('1') ! Do nothing. + case default + call MOM_error(WARNING, "MOM_io: get_var_dimension_metadata: "//& + "Unrecognized t_grid "//trim(t_grid)) + end select +end subroutine get_var_dimension_metadata + + +!> Populate the axis_data structure with axis data and attributes for diagnostic and restart files +subroutine MOM_get_diagnostic_axis_data(axis_data_CS, axis_name, axis_number, G, dG, GV, time_val, time_units) + + type(axis_data_type), intent(inout) :: axis_data_CS !< structure containing the axis data and metadata + character(len=*), intent(in) :: axis_name !< name of the axis + integer, intent(in) :: axis_number !< positional value (wrt to file) of the axis to register + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the file uses any + !! horizontal grid axes. + type(verticalGrid_type), target, optional, intent(in) :: GV !< ocean vertical grid structure + real,dimension(:), target, optional, intent(in) :: time_val !< time value + character(len=*), optional,intent(in) :: time_units!< units for non-periodic time axis + ! local + character(len=40) :: x_axis_units='', y_axis_units='' + integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB + real, pointer, dimension(:) :: gridLatT => NULL(), & ! The latitude or longitude of T or B points for + gridLatB => NULL(), & ! the purpose of labeling the output axes. + gridLonT => NULL(), & + gridLonB => NULL() + + ! initialize axis_data_CS elements + axis_data_CS%axis(axis_number)%name = '' + axis_data_CS%axis(axis_number)%longname = '' + axis_data_CS%axis(axis_number)%units = '' + axis_data_CS%axis(axis_number)%horgrid_position = 0 + axis_data_CS%axis(axis_number)%is_domain_decomposed = .false. + axis_data_CS%axis(axis_number)%positive = '' + axis_data_CS%data(axis_number)%p => NULL() + + ! set the ocean grid coordinates and metadata + if (present(G)) then + gridLatT => G%gridLatT ; gridLatB => G%gridLatB + gridLonT => G%gridLonT ; gridLonB => G%gridLonB + x_axis_units = G%x_axis_units ; y_axis_units = G%y_axis_units + isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg + IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB + elseif (present(dG)) then + gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB + gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB + x_axis_units = dG%x_axis_units ; y_axis_units = dG%y_axis_units + isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg + IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB + endif + + select case(trim(lowercase(axis_name))) + case('lath') + if (associated(gridLatT)) & + axis_data_CS%data(axis_number)%p=>gridLatT(jsg:jeg) + + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Latitude' + axis_data_CS%axis(axis_number)%units = y_axis_units + axis_data_CS%axis(axis_number)%horgrid_position = CENTER + axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. + case('lonh') + if (associated(gridLonT)) & + axis_data_CS%data(axis_number)%p=>gridLonT(isg:ieg) + + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%horgrid_position = CENTER + axis_data_CS%axis(axis_number)%longname = 'Longitude' + axis_data_CS%axis(axis_number)%units = x_axis_units + axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. + case('latq') + if (associated(gridLatB)) & + axis_data_CS%data(axis_number)%p=>gridLatB(JsgB:JegB) + + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Latitude' + axis_data_CS%axis(axis_number)%units = y_axis_units + axis_data_CS%axis(axis_number)%horgrid_position = NORTH_FACE + axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. + case('lonq') + if (associated(gridLonB)) & + axis_data_CS%data(axis_number)%p=>gridLonB(IsgB:IegB) + + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Longitude' + axis_data_CS%axis(axis_number)%units = x_axis_units + axis_data_CS%axis(axis_number)%horgrid_position = EAST_FACE + axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. + case('layer') + if (present(GV)) then + axis_data_CS%data(axis_number)%p=>GV%sLayer(1:GV%ke) + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Layer pseudo-depth, -z*' + axis_data_CS%axis(axis_number)%units = GV%zAxisUnits + axis_data_CS%axis(axis_number)%positive = 'up' + endif + case('interface') + if (present(GV)) then + axis_data_CS%data(axis_number)%p=>GV%sInterface(1:GV%ke+1) + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Interface pseudo-depth, -z*' + axis_data_CS%axis(axis_number)%units = GV%zAxisUnits + axis_data_CS%axis(axis_number)%positive = 'up' + endif + case('time') + if (.not.(present(time_val))) & + call MOM_error(FATAL, "MOM_io::get_diagnostic_axis_data: requires time_val"//& + " and time_units arguments for "//trim(axis_name)) + + axis_data_CS%data(axis_number)%p=>time_val + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Time' + + if (present(time_units)) then + axis_data_CS%axis(axis_number)%units = time_units + else + axis_data_CS%axis(axis_number)%units = 'days' + endif + case('period') + if (.not.(present(time_val))) & + call MOM_error(FATAL, "MOM_axis::get_diagnostic_axis_data: requires a time_val argument "// & + "for "//trim(axis_name)) + axis_data_CS%data(axis_number)%p=>time_val + axis_data_CS%axis(axis_number)%name = trim(axis_name) + axis_data_CS%axis(axis_number)%longname = 'Periods for cyclical variables' + case default + call MOM_error(WARNING, "MOM_axis::get_diagnostic_axis_data:"//trim(axis_name)//" is an unrecognized axis") + end select + +end subroutine MOM_get_diagnostic_axis_data + + +!> set the logical variables that determine which diagnositic axes to use +subroutine get_horizontal_grid_logic(grid_string_id, use_lath, use_lonh, use_latq, use_lonq) + character(len=*), intent(in) :: grid_string_id !< horizontal grid string + logical, intent(out) :: use_lath !< if .true., y-axis is oriented in CENTER position + logical, intent(out) :: use_lonh !< if .true., x-axis is oriented in CENTER position + logical, intent(out) :: use_latq !< if .true., y-axis is oriented in NORTH_FACE position + logical, intent(out) :: use_lonq !< if .true., x-axis is oriented in EAST_FACE position + + use_lath = .false. + use_lonh = .false. + use_latq = .false. + use_lonq = .false. + select case (trim(grid_string_id)) + case ('h') ; use_lath = .true. ; use_lonh = .true. ! x=CENTER, y=CENTER + case ('q') ; use_latq = .true. ; use_lonq = .true. ! x=EAST_FACE, y=NORTH_FACE + case ('u') ; use_lath = .true. ; use_lonq = .true. ! x=EAST_FACE, y=CENTER + case ('v') ; use_latq = .true. ; use_lonh = .true. ! x=CENTER, y=NORTH_FACE + case ('T') ; use_lath = .true. ; use_lonh = .true. ! x=CENTER, y=CENTER + case ('Bu') ; use_latq = .true. ; use_lonq = .true. ! x=EAST_FACE, y=NORTH_FACE + case ('Cu') ; use_lath = .true. ; use_lonq = .true. ! x=EAST_FACE, y=CENTER + case ('Cv') ; use_latq = .true. ; use_lonh = .true. ! x=CENTER, y=NORTH_FACE + case ('1') ; ! x=0, y=0 + case default + call MOM_error(FATAL, "MOM_axis:get_var_dimension_features "//& + "Unrecognized hor_grid argument "//trim(grid_string_id)) + end select +end subroutine get_horizontal_grid_logic + +!> Define the time units for the input time value +function get_time_units(time_value) result(time_units_out) + real, intent(in) :: time_value !< numerical time value in seconds + !! i.e., before dividing by 86400. + ! local + character(len=10) :: time_units ! time units + character(len=10) :: time_units_out ! time units trimmed + time_units = '' + time_units_out = '' + if (time_value < 0.0) then + time_units = "days" ! The default value. + elseif (mod(time_value,86400.0)==0.0) then + time_units = "days" + elseif ((time_value >= 0.99) .and. (time_value < 1.01)) then + time_units = "seconds" + elseif ((time_value >= 3599.0) .and. (time_value < 3601.0)) then + time_units = "hours" + elseif ((time_value >= 86399.0) .and. (time_value < 86401.0)) then + time_units = "days" + elseif ((time_value >= 3.0e7) .and. (time_value < 3.2e7)) then + time_units = "years" + else + write(time_units,'(es8.2," s")') time_value + endif + time_units_out = trim(time_units) +end function get_time_units + +!> function to get the index of a time_value from a netCDF file +function get_time_index(filename, time_to_find) result (time_index) + character(len=*) :: filename ! name of the file to read in + real, intent(in) :: time_to_find ! time value to search for in file + ! local + type(fmsNetcdfFile_t) :: fileobj ! netCDF file object returned by open_file + real, allocatable, dimension(:) :: file_times ! array of time values read from file + integer :: dim_unlim_size, i, time_index + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + logical :: file_open_success + + time_index = 1 + dim_unlim_size = 0 + dim_unlim_name = "" + file_open_success = .false. + + if (.not. check_if_open(fileobj)) & + !call MOM_error(FATAL, "get_time_index_nodd: netcdf file object must be open.") + file_open_success=fms2_open_file(fileobj, trim(filename), "read", is_restart=.false.) + + call get_unlimited_dimension_name(fileobj, dim_unlim_name) + call get_dimension_size(fileObj, trim(dim_unlim_name), dim_unlim_size) + ! time index will be one more than the unlimited dimension size if the time_to_find is not in the file + if (dim_unlim_size .gt. 0) then + time_index = dim_unlim_size+1 + allocate(file_times(dim_unlim_size)) + call read_data(fileobj,trim(dim_unlim_name), file_times) + + do i=1,dim_unlim_size + if (ABS(file_times(i)-time_to_find) .gt. TINY(time_to_find)) then + continue + else + time_index = i + exit + endif + enddo + deallocate(file_times) + endif + if (check_if_open(fileobj)) call fms2_close_file(fileobj) +end function get_time_index + +!> register axes associated with a variable from a domain-decomposed netCDF file that are mapped to +!! a sub-domain (e.g., a supergrid). +!> \note The user must specify units for variables with longitude/x-axis and/or latitude/y-axis axes to obtain +!! the correct domain decomposition for the data buffer. +subroutine MOM_register_variable_axes_subdomain(fileObj, variableName, io_domain, xPosition, yPosition) + type(FmsNetcdfFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file + character(len=*), intent(in) :: variableName !< name of the variable + type(domain2d), intent(in) :: io_domain !< type that contains the mpp io domain + integer, intent(in), optional :: xPosition !< domain position of the x-axis + integer, intent(in), optional :: yPosition !< domain position of the y-axi + ! local + character(len=40) :: units ! units corresponding to a specific variable dimension + character(len=40), allocatable, dimension(:) :: dim_names ! variable dimension names + integer :: i, isg, ieg, isc, iec, jsg, jeg, jsc, jec, xlen, ylen + integer :: ndims ! number of dimensions + integer :: xPos, yPos, pos ! domain positions for x and y axes. Default is CENTER + integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes + + if (.not. check_if_open(fileObj)) call MOM_error(FATAL,"MOM_axis:register_variable_axes_subdomain: The fileObj "// & + " has not been opened. Call fms2_open_file(fileObj,...) "// & + "before passing the fileObj argument to this function.") + xPos=CENTER + yPos=CENTER + if (present(xPosition)) xPos=xPosition + if (present(yPosition)) yPos=yPosition + ! get variable dimension names and lengths + ndims = get_variable_num_dimensions(fileObj, trim(variableName)) + allocate(dimSizes(ndims)) + allocate(dim_names(ndims)) + call get_variable_size(fileObj, trim(variableName), dimSizes, broadcast=.true.) + call get_variable_dimension_names(fileObj, trim(variableName), dim_names) + ! determine the position to pass to the mpp domain calls + if (xPos .eq. EAST_FACE) then + if (yPos .eq. NORTH_FACE) then + pos = CORNER + else + pos = EAST_FACE + endif + elseif (yPos .eq. NORTH_FACE) then + pos = NORTH_FACE + endif + ! Get the lengths of the global indicies + call mpp_get_compute_domain(io_domain, xsize=xlen, ysize=ylen, position=pos) + ! register the axes + !>\note: This is not a comprehensive check for all possible supported horizontal axes associated with variables + !! read from netCDF files. Developers should add/remove cases as needed. + do i=1,ndims + !if (.not.(is_dimension_registered(fileObj, trim(dim_names(i))))) then + select case(trim(lowercase(dim_names(i)))) + case ("grid_x_t") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case ("nx") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("nxp") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("longitude") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("long") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("lon") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("lonh") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("lonq") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case("xh") + call register_axis(fileObj, trim(dim_names(i)), xlen) + case ("grid_y_t") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case ("ny") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("nyp") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("latitude") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("lat") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("lath") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("latq") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case("yh") + call register_axis(fileObj, trim(dim_names(i)), ylen) + case default ! assumes that the axis is not domain-decomposed + if (.not. is_dimension_unlimited(fileObj, trim(dim_names(i)))) & + call MOM_error(WARNING,"MOM_register_variable_axes_subdomain: the axis "//trim(dim_names(i))//& + "is not included in the valid x and y dimension cases. If the code hangs, check the whether "//& + "an x or y axis is being registered as a non-domain-decomposed variable, "//& + "and add it to the accepted cases if necessary.") + call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) + end select + ! endif + enddo + + if (allocated(dimSizes)) deallocate(dimSizes) + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_register_variable_axes_subdomain + +!> register axes associated with a variable from a domain-decomposed netCDF file +!> @note The user must specify units for variables with longitude/x-axis and/or latitude/y-axis axes +!! to obtain the correct domain decomposition for the data buffer. +subroutine MOM_register_variable_axes_full(fileObj, variableName, xPosition, yPosition) + type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file + character(len=*), intent(in) :: variableName !< name of the variable + integer, intent(in), optional :: xPosition !< domain position of the x-axis + integer, intent(in), optional :: yPosition !< domain position of the y-axis + ! local + character(len=40) :: units ! units corresponding to a specific variable dimension + character(len=40), allocatable, dimension(:) :: dim_names ! variable dimension names + integer :: i + integer :: ndims ! number of dimensions + integer :: xPos, yPos ! domain positions for x and y axes. Default is CENTER + integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes + + if (.not. check_if_open(fileObj)) call MOM_error(FATAL,"MOM_axis:register_variable_axes: The fileObj has "// & + "not been opened. Call fms2_open_file(fileObj,...) before "// & + "passing the fileObj argument to this function.") + xPos=CENTER + yPos=CENTER + if (present(xPosition)) xPos=xPosition + if (present(yPosition)) yPos=yPosition + ! get variable dimension names and lengths + ndims = get_variable_num_dimensions(fileObj, trim(variableName)) + allocate(dimSizes(ndims)) + allocate(dim_names(ndims)) + call get_variable_size(fileObj, trim(variableName), dimSizes) + call get_variable_dimension_names(fileObj, trim(variableName), dim_names) + ! register the axes + !>@note: This is not a comprehensive check for all possible supported horizontal axes associated with variables + !! read from netCDF files. Developers should add/remove cases as needed. + do i=1,ndims + if (.not.(is_dimension_registered(fileobj, trim(dim_names(i))))) then + select case(trim(lowercase(dim_names(i)))) + case ("grid_x_t") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case ("nx") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("nxp") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("longitude") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("long") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("lon") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("lonh") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("lonq") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("xh") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case("i") + call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) + case ("grid_y_t") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case ("ny") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("nyp") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("latitude") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("lat") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("lath") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("latq") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("yh") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case("j") + call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) + case default ! assumes that the axis is not domain-decomposed + if (.not. is_dimension_unlimited(fileObj, trim(dim_names(i)))) & + call MOM_error(WARNING,"MOM_register_variable_axes_full: the axis "//trim(dim_names(i))//" is not "//& + "included in the valid x and y dimension cases. If the code hangs, check the whether "//& + "an x or y axis is being registered as a non-domain-decomposed variable, "//& + "and add it to the accepted cases if necessary.") + call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) + end select + endif + enddo + + deallocate(dimSizes) + deallocate(dim_names) +end subroutine MOM_register_variable_axes_full + + +!> convert the variable checksum integer(s) to a single string +!! If there is more than 1 checksum, commas are inserted between +!! each checksum value in the output string +function convert_checksum_to_string(checksum_int) result (checksum_string) + integer(kind=8), intent(in) :: checksum_int !< checksum integer values +! local + character(len=64) :: checksum_string + integer :: i + + checksum_string = '' + + write (checksum_string,'(Z16)') checksum_int ! Z16 is the hexadecimal format code + +end function convert_checksum_to_string + + +end module MOM_axis diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index c516c96e86..6768e47dfa 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -4,13 +4,17 @@ module MOM_io ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_axis, only : MOM_get_diagnostic_axis_data, MOM_register_diagnostic_axis +use MOM_axis, only : axis_data_type, get_time_index, get_var_dimension_metadata +use MOM_axis, only : get_time_units, convert_checksum_to_string use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING use MOM_domains, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE use MOM_domains, only : get_simple_array_i_ind, get_simple_array_j_ind use MOM_file_parser, only : log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_dyn_horgrid, only : dyn_horgrid_type -use MOM_string_functions, only : lowercase, slasher +use MOM_string_functions, only : lowercase, slasher, append_substring +use MOM_time_manager, only : time_type, time_type_to_real use MOM_verticalGrid, only : verticalGrid_type use ensemble_manager_mod, only : get_ensemble_id @@ -18,8 +22,12 @@ module MOM_io use fms_io_mod, only : file_exist, field_size, read_data use fms_io_mod, only : field_exists => field_exist, io_infra_end=>fms_io_exit use fms_io_mod, only : get_filename_appendix => get_filename_appendix +use mpp_mod, only : mpp_pe, mpp_max, mpp_npes use mpp_domains_mod, only : domain1d, domain2d, mpp_get_domain_components use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST +use mpp_domains_mod, only : mpp_get_data_domain, mpp_get_compute_domain, mpp_get_global_domain +use mpp_domains_mod, only : mpp_get_domain_npes, mpp_define_io_domain, mpp_get_io_domain +use mpp_domains_mod, only : mpp_get_io_domain_layout use mpp_io_mod, only : open_file => mpp_open, close_file => mpp_close use mpp_io_mod, only : mpp_write_meta, write_field => mpp_write, mpp_get_info use mpp_io_mod, only : mpp_get_atts, mpp_get_axes, get_axis_data=>mpp_get_axis_data, axistype @@ -33,6 +41,28 @@ module MOM_io use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times use mpp_io_mod, only : io_infra_init=>mpp_io_init +! fms2_io +use fms2_io_mod, only : check_if_open, get_dimension_names,get_dimension_size +use fms2_io_mod, only : get_compute_domain_dimension_indices, get_global_attribute +use fms2_io_mod, only : get_global_io_domain_indices, get_num_dimensions, get_num_variables +use fms2_io_mod, only : get_unlimited_dimension_name, get_variable_dimension_names, get_variable_names +use fms2_io_mod, only : get_variable_num_dimensions, get_variable_size, get_variable_units +use fms2_io_mod, only : get_variable_unlimited_dimension_index, global_att_exists, is_dimension_unlimited +use fms2_io_mod, only : is_dimension_registered, register_restart_field, register_axis +use fms2_io_mod, only : register_field, register_variable_attribute, fms2_open_file => open_file +use fms2_io_mod, only : fms2_close_file => close_file, write_data, attribute_exists => variable_att_exists +use fms2_io_mod, only : dimension_exists, variable_exists, fms2_io_file_exists => file_exists +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited + +!use, intrinsic :: iso_fortran_env + +!NOTE: uncomment when ready to replace mpp_read calls +!use MOM_read_data_fms2, only : MOM_read_data_4d_DD, MOM_read_data_3d_DD, MOM_read_data_2d_DD +!use MOM_read_data_fms2, only : MOM_read_data_1d_DD, MOM_read_data_scalar +!use MOM_read_data_fms2, only : MOM_read_data_4d_noDD, MOM_read_data_3d_noDD, MOM_read_data_2d_noDD +!use MOM_read_data_fms2, only : MOM_read_data_1d_noDD, MOM_read_vector_2d_fms2, MOM_read_vector_3d_fms2 + +! use MOM_write_field_fms2, only : write_field !NOTE: uncomment when ready to replace mpp_write calls use netcdf implicit none ; private @@ -66,22 +96,30 @@ module MOM_io !> Indicate whether a file exists, perhaps with domain decomposition interface file_exists + module procedure FMS2_file_exists module procedure FMS_file_exists module procedure MOM_file_exists end interface -!> Read a data field from a file +!> Read a pair of data fields representing the two components of a vector from a file +interface MOM_read_vector + module procedure MOM_read_vector_3d + module procedure MOM_read_vector_2d +end interface + +!> interface to read data from a netcdf file interface MOM_read_data module procedure MOM_read_data_4d module procedure MOM_read_data_3d module procedure MOM_read_data_2d module procedure MOM_read_data_1d -end interface +end interface MOM_read_data -!> Read a pair of data fields representing the two components of a vector from a file -interface MOM_read_vector - module procedure MOM_read_vector_3d - module procedure MOM_read_vector_2d +!> Open a netcdf file in write or overwrite mode using the fms-io or fms2-io netcdf interfaces +interface create_file + module procedure create_file_old + module procedure create_file_fms2_filename + module procedure create_file_fms2_fileobj end interface contains @@ -89,7 +127,7 @@ module MOM_io !> Routine creates a new NetCDF file. It also sets up !! structures that describe this file and variables that will !! later be written to this file. Type for describing a variable, typically a tracer -subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV, checksums) +subroutine create_file_old(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV, checksums) integer, intent(out) :: unit !< unit id of an open file or -1 on a !! nonwriting PE with single file output character(len=*), intent(in) :: filename !< full path to the file to create @@ -342,7 +380,463 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit if (use_int) call write_field(unit, axis_int) if (use_periodic) call write_field(unit, axis_periodic) -end subroutine create_file +end subroutine create_file_old + + +!> This routine opens a netcdf file in "write" or "overwrite" mode, registers the global diagnostic axes, and writes +!! the axis data and metadata to the file +subroutine create_file_fms2_filename(filename, vars, numVariables, use_fms2, register_time, G, DG, GV, checksums, & + is_restart) + character(len=*), intent(in) :: filename !< full path to the netcdf file + type(vardesc), dimension(:), intent(in) :: vars !< structures describing the output + integer, intent(in) :: numVariables !< number of variables to write to the file + logical, intent(in) :: use_fms2 !< flag indicating whether to use this routine + logical, optional, intent(in) :: register_time !< if .true., register a time dimension to the file + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums(:,:) !< checksums of the variables + logical, optional, intent(in) :: is_restart !< indicates whether file is a restart file + + ! local + type(FmsNetcdfFile_t) :: fileObjNoDD ! non-domain-decomposed netcdf file object returned by open_file + type(FmsNetcdfDomainFile_t) :: fileObjDD ! domain-decomposed netcdf file object returned by open_file + type(axis_data_type) :: axis_data_CS ! structure for coordinate variable metadata + type(MOM_domain_type), pointer :: Domain => NULL() + logical :: file_open_successDD, file_open_successNoDD ! true if netcdf file is opened + logical :: one_file, domain_set ! indicates whether the file will be domain-decomposed or not + logical :: reg_time ! register the time if .true. + logical :: is_restart_file + character(len=10) :: nc_mode + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=1024) :: filename_temp + character(len=48), allocatable, dimension(:,:) :: dim_names ! variable dimension names + integer :: i, is, ie, j, substring_index, total_axes + integer :: num_dims ! number of dimensions + integer :: thread ! indicates whether threading is used + integer, dimension(4) :: dim_lengths ! variable dimension lengths + integer, allocatable :: pelist(:) ! list of pes associated with the file + real :: time + + ! determine whether the file will be domain-decomposed or not + domain_set=.false. + if (present(G)) then + domain_set = .true. ; Domain => G%Domain + elseif (present(dG)) then + domain_set = .true. ; Domain => dG%Domain + endif + + is_restart_file = .false. + if (present(is_restart)) is_restart_file = is_restart + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index < 1) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + + nc_mode = "" + if (file_exists(trim(filename_temp), .true.)) then + nc_mode = "overwrite" + else + nc_mode = "write" + endif + + reg_time = .false. + if (present(register_time)) reg_time = register_time + + ! open the file + file_open_successNoDD=.false. + file_open_successDD=.false. + + if (domain_set) then + ! define the io domain if on one pe and the io domain is not set + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + + if (.not. check_if_open(fileObjDD)) & + file_open_successDD=fms2_open_file(fileObjDD, filename_temp, trim(nc_mode), Domain%mpp_domain, & + is_restart=is_restart_file) + else + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + + if (.not. check_if_open(fileObjNoDD)) & + file_open_successNoDD=fms2_open_file(fileObjNoDD, filename_temp, trim(nc_mode), & + is_restart=is_restart_file, pelist=pelist) + endif + ! allocate the output data variable dimension attributes + allocate(dim_names(numVariables,4)) + dim_names(:,:) = "" + ! allocate the axis data and attribute types for the file + !> \note The user should increase the sizes of the axis and data attributes to accommodate more axes if necessary. + allocate(axis_data_CS%axis(7)) + allocate(axis_data_CS%data(7)) + ! axis registration procedure for the domain-decomposed case + if (file_open_successDD) then + do i=1,numVariables + num_dims=0 + dim_lengths(:) = 0 + + !> \note The time dimension is registered separately at the end of the procedure if reg_time = .true. + !! so the t_grid argument in get_var_dimension_metadata is set to '1' (do nothing) + if (present(G)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, dG=dG) + endif + + if(present(GV)) & + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, GV=GV) + !> \note num_dims will be 0 for scalar values + if (num_dims .le. 0) cycle + + do j=1,num_dims + ! register the variable axes to the file if they are not already registered + if (dim_lengths(j) .gt. 0) then + if (.not.(dimension_exists(fileObjDD, dim_names(i,j)))) then + + if (present(G)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G) + endif + elseif (present(dG)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG) + endif + elseif (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, GV=GV) + endif + call MOM_register_diagnostic_axis(fileObjDD, trim(dim_names(i,j)), dim_lengths(j)) + endif + ! register the axis attributes and write the axis data to the file + if (.not.(variable_exists(fileObjDD, trim(axis_data_CS%axis(j)%name)))) then + if (associated(axis_data_CS%data(j)%p)) then + + call register_field(fileObjDD, trim(axis_data_CS%axis(j)%name), & + "double", dimensions=(/trim(axis_data_CS%axis(j)%name)/)) + + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'long_name', axis_data_CS%axis(j)%longname) + + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'units', trim(axis_data_CS%axis(j)%units)) + + if (len_trim(axis_data_CS%axis(j)%positive)>1) & + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'positive', trim(axis_data_CS%axis(j)%positive)) + + if (axis_data_CS%axis(j)%is_domain_decomposed) then + call get_global_io_domain_indices(fileObjDD, trim(axis_data_CS%axis(j)%name), is, ie) + call write_data(fileObjDD, trim(axis_data_CS%axis(j)%name), axis_data_CS%data(j)%p(is:ie)) + else + call write_data(fileObjDD, trim(axis_data_CS%axis(j)%name), axis_data_CS%data(j)%p) + endif + endif + endif + endif + enddo + enddo + + if (reg_time) then + if (.not.(dimension_exists(fileObjDD,"Time"))) & + call register_axis(fileObjDD, "Time", unlimited) + endif + + if (check_if_open(fileObjDD)) call fms2_close_file(fileObjDD) + ! axis registration and write procedure for the non-domain-decomposed case + elseif (file_open_successNoDD) then + do i=1,numVariables + num_dims=0 + dim_lengths(:) = 0 + + !> \note The time dimension is registered separately at the end of the procedure if reg_time = .true. + !! so the t_grid argument in get_var_dimension_metadata is set to '1' (do nothing) + if (present(G)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, dG=dG) + endif + + if(present(GV)) & + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, GV=GV) + !> \note num_dims will be 0 for scalar variables + if (num_dims .le. 0) cycle + + do j=1,num_dims + ! register the variable axes to the file if they are not already registered + if (dim_lengths(j) .gt. 0) then + if (.not.(dimension_exists(fileObjNoDD, dim_names(i,j)))) then + if (present(G)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G) + endif + elseif (present(dG)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG) + endif + elseif (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, GV=GV) + endif + call register_axis(fileObjNoDD, trim(dim_names(i,j)), dim_lengths(j)) + endif + ! register the axis attributes and write the axis data to the file + if (.not.(variable_exists(fileObjNoDD, trim(axis_data_CS%axis(j)%name)))) then + if (associated(axis_data_CS%data(j)%p)) then + call register_field(fileObjNoDD, trim(axis_data_CS%axis(j)%name), & + "double", dimensions=(/trim(axis_data_CS%axis(j)%name)/)) + + call register_variable_attribute(fileObjNoDD, trim(axis_data_CS%axis(j)%name), & + 'long_name', axis_data_CS%axis(j)%longname) + + call register_variable_attribute(fileObjNoDD, trim(axis_data_CS%axis(j)%name), & + 'units', trim(axis_data_CS%axis(j)%units)) + + if (len_trim(axis_data_CS%axis(j)%positive)>1) & + call register_variable_attribute(fileObjNoDD, trim(axis_data_CS%axis(j)%name), & + 'positive', trim(axis_data_CS%axis(j)%positive)) + + if (lowercase(trim(axis_data_CS%axis(j)%name)) .ne. 'time') then + call write_data(fileObjNoDD, trim(axis_data_CS%axis(j)%name), axis_data_CS%data(j)%p) + endif + endif + endif + endif + enddo + enddo + + if (reg_time) then + if (.not.(dimension_exists(fileObjNoDD,"Time"))) & + call register_axis(fileObjNoDD, "Time" , unlimited) + endif + + if (check_if_open(fileObjNoDD)) call fms2_close_file(fileObjNoDD) + endif + + deallocate(dim_names) + deallocate(axis_data_CS%axis) + deallocate(axis_data_CS%data) + if (allocated(pelist)) deallocate(pelist) + nullify(Domain) + +end subroutine create_file_fms2_filename + +!> This routine opens a netcdf file in "write" or "overwrite" mode, registers the global diagnostic axes, and writes +!! the axis data and metadata to the file. It returns the netcdf file object for additional writing. +subroutine create_file_fms2_fileobj(filename, fileObjDD, vars, numVariables, register_time, G, DG, GV, & + checksums, is_restart) + character(len=*), intent(in) :: filename !< full path to the netcdf file + type(FmsNetcdfDomainFile_t), intent(inout) :: fileObjDD !< domain-decomposed netcdf file object + !! returned by open_file + type(vardesc), dimension(:), intent(in) :: vars !< structures describing the output + integer, intent(in) :: numVariables !< number of variables to write to the file + logical, optional, intent(in) :: register_time !< if .true., register a time dimension to the file + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums(:,:) !< checksums of the variables + logical, optional, intent(in) :: is_restart !< indicates whether file is a restart file + + ! local + type(axis_data_type) :: axis_data_CS ! structure for coordinate variable metadata + type(MOM_domain_type), pointer :: Domain => NULL() + logical :: file_open_successDD ! true if netcdf file is opened + logical :: one_file, domain_set ! indicates whether the file will be domain-decomposed or not + logical :: reg_time ! register the time if .true. + logical :: is_restart_file ! .true. if the file is a restart file + character(len=10) :: nc_mode + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=1024) :: filename_temp + character(len=48), allocatable, dimension(:,:) :: dim_names ! variable dimension names + integer :: i, is, ie, j, substring_index, total_axes + integer :: num_dims ! number of dimensions + integer :: thread ! indicates whether threading is used + integer, dimension(4) :: dim_lengths ! variable dimension lengths + integer, allocatable :: pelist(:) ! list of pes associated with the file + real :: time + + ! determine whether the file will be domain-decomposed or not + domain_set=.false. + if (present(G)) then + domain_set = .true. ; Domain => G%Domain + elseif (present(dG)) then + domain_set = .true. ; Domain => dG%Domain + endif + + is_restart_file = .false. + if (present(is_restart)) is_restart_file = is_restart + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index < 1) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + + nc_mode = "" + if (file_exists(trim(filename_temp), .true.)) then + nc_mode = "overwrite" + else + nc_mode = "write" + endif + + reg_time = .false. + if (present(register_time)) reg_time = register_time + ! open the file + file_open_successDD=.false. + ! define the io domain if on one pe and the io domain is not set + if (domain_set) then + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + else + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + if (.not. check_if_open(fileObjDD)) & + !write(output_unit, '(A)'), "Create_file: Opening file ", trim(filename_temp) + file_open_successDD=fms2_open_file(fileObjDD, filename_temp, trim(nc_mode), Domain%mpp_domain, & + is_restart=is_restart_file) + ! allocate the output data variable dimension attributes + allocate(dim_names(numVariables,4)) + dim_names(:,:) = "" + ! allocate the axis data and attribute types for the file + !> \note The user should increase the sizes of the axis and data attributes to accommodate more axes if necessary. + allocate(axis_data_CS%axis(7)) + allocate(axis_data_CS%data(7)) + ! axis registration procedure for the domain-decomposed case + if (file_open_successDD) then + do i=1,numVariables + num_dims=0 + dim_lengths(:) = 0 + !> \note The time dimension is registered separately at the end of the procedure if reg_time = .true. + !! so the t_grid argument in get_var_dimension_metadata is set to '1' (do nothing) + if (present(G)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, dG=dG) + endif + + if(present(GV)) & + call get_var_dimension_metadata(vars(i)%hor_grid, vars(i)%z_grid, '1', dim_names(i,:), & + dim_lengths, num_dims, GV=GV) + !> \note num_dims will be 0 for scalar values + if (num_dims .le. 0) cycle + + do j=1,num_dims + ! register the variable axes to the file if they are not already registered + if (dim_lengths(j) .gt. 0) then + if (.not.(dimension_exists(fileObjDD, dim_names(i,j)))) then + + if (present(G)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, G=G) + endif + elseif (present(dG)) then + if (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG, GV=GV) + else + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, dG=dG) + endif + elseif (present(GV)) then + call MOM_get_diagnostic_axis_data(axis_data_CS, dim_names(i,j), j, GV=GV) + endif + call MOM_register_diagnostic_axis(fileObjDD, trim(dim_names(i,j)), dim_lengths(j)) + endif + ! register the axis attributes and write the axis data to the file + if (.not.(variable_exists(fileObjDD, trim(axis_data_CS%axis(j)%name)))) then + if (associated(axis_data_CS%data(j)%p)) then + + call register_field(fileObjDD, trim(axis_data_CS%axis(j)%name), & + "double", dimensions=(/trim(axis_data_CS%axis(j)%name)/)) + + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'long_name', axis_data_CS%axis(j)%longname) + + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'units', trim(axis_data_CS%axis(j)%units)) + + if (len_trim(axis_data_CS%axis(j)%positive)>1) & + call register_variable_attribute(fileObjDD, trim(axis_data_CS%axis(j)%name), & + 'positive', trim(axis_data_CS%axis(j)%positive)) + + if (axis_data_CS%axis(j)%is_domain_decomposed) then + call get_global_io_domain_indices(fileObjDD, trim(axis_data_CS%axis(j)%name), is, ie) + call write_data(fileObjDD, trim(axis_data_CS%axis(j)%name), axis_data_CS%data(j)%p(is:ie)) + else + call write_data(fileObjDD, trim(axis_data_CS%axis(j)%name), axis_data_CS%data(j)%p) + endif + endif + endif + endif + enddo + enddo + + if (reg_time) then + if (.not.(dimension_exists(fileObjDD,"Time"))) & + call register_axis(fileObjDD, "Time", unlimited) + endif + else + call MOM_error(FATAL, "MOM_io::create_file_fms2_filobj: unable to open file "//trim(filename)) + endif + + deallocate(dim_names) + deallocate(axis_data_CS%axis) + deallocate(axis_data_CS%data) + if (allocated(pelist)) deallocate(pelist) + nullify(Domain) + +end subroutine create_file_fms2_fileobj !> This routine opens an existing NetCDF file for output. If it @@ -844,6 +1338,19 @@ function FMS_file_exists(filename, domain, no_domain) end function FMS_file_exists +!> Returns true if the named file exists +function FMS2_file_exists(filename, use_fms2) + character(len=*), intent(in) :: filename !< The name of the file being inquired about + logical, intent(in) :: use_fms2 !< flag indicating to use the fms2-io interface +! This function uses the fms2_io function file_exists to determine whether +! a named file (or its decomposed variant) exists. + + logical :: FMS2_file_exists + + FMS2_file_exists = fms2_io_file_exists(filename) + +end function FMS2_file_exists + !> This function uses the fms_io function read_data to read 1-D !! data field named "fieldname" from file "filename". subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale) diff --git a/src/framework/MOM_read_data_fms2.F90 b/src/framework/MOM_read_data_fms2.F90 new file mode 100644 index 0000000000..d15d5a3085 --- /dev/null +++ b/src/framework/MOM_read_data_fms2.F90 @@ -0,0 +1,1540 @@ +!> This module contains routines that wrap the fms2 read_data calls +module MOM_read_data_fms2 + +! This file is part of MOM6. See LICENSE.md for the license. +use MOM_axis, only : MOM_register_variable_axes +use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING +use MOM_domains, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE +use MOM_domains, only : get_simple_array_i_ind, get_simple_array_j_ind +use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_string_functions, only : lowercase +use MOM_verticalGrid, only : verticalGrid_type +use fms2_io_mod, only : read_data, attribute_exists => variable_att_exists, variable_exists +use fms2_io_mod, only : register_field, register_variable_attribute, fms2_open_file => open_file +use fms2_io_mod, only : fms2_close_file => close_file, write_data, get_variable_dimension_names +use fms2_io_mod, only : check_if_open, get_dimension_names, get_dimension_size +use fms2_io_mod, only : is_dimension_registered, register_axis, get_variable_size +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited, get_variable_names +use fms2_io_mod, only : get_variable_num_dimensions, get_variable_units, is_dimension_unlimited +use fms2_io_mod, only : get_num_variables +use mpp_domains_mod, only : domain2d +use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST +use mpp_domains_mod, only : mpp_get_domain_npes, mpp_define_io_domain, mpp_get_io_domain +use mpp_domains_mod, only : mpp_get_data_domain, mpp_get_compute_domain, mpp_get_global_domain + +implicit none ; private + +public MOM_read_data_scalar, MOM_read_vector_2d_fms2, MOM_read_vector_3d_fms2 +public MOM_read_data_4d_noDD, MOM_read_data_3d_noDD, MOM_read_data_2d_noDD, MOM_read_data_1d_noDD +public MOM_read_data_4d_DD, MOM_read_data_3d_DD, MOM_read_data_2d_DD, MOM_read_data_1d_DD + +! CAUTION: The following variables are saved by default, and are only necessary for consecutive calls to +! MOM_read_data with the same file name. The user should ensure that fms2_close_file on +! the fileobj_read structures are called at every requisite time step at after the last +! variable is written to the file by omitting the optional leave_file_open argument, or setting it to .false. + +!> netCDF domain-decomposed file object returned by call to +!! open_file in MOM_read_data_DD calls +type(FmsNetcdfDomainFile_t), private :: fileobj_read_dd + +!> netCDF domain-decomposed file object returned by call to +!! open_file in MOM_read_data_noDD calls +type(FmsNetcdfFile_t), private :: fileobj_read + +!> Type with variable metadata for a netCDF file opened to read domain-decomposed data +type file_variable_meta_DD + integer :: nvars = 0!< number of variables in a netCDF file opened to read domain-decomposed data + character(len=96), allocatable, dimension(:) :: var_names !< array for names of variables in a netCDF + !! file opened to read domain-decomposed data +end type file_variable_meta_DD + +!> Type with variable metadata for a netCDF file opened to read non-domain-decomposed data +type file_variable_meta_noDD + integer :: nvars = 0 !< number of variables in a netCDF file opened to read non-domain-decomposed data + character(len=96), allocatable, dimension(:) :: var_names !< array for names of variables in a netCDF + !! file opened to read non-domain-decomposed data +end type file_variable_meta_noDD +!> type to hold metadata for variables in a domain-decomposed file +type (file_variable_meta_DD), private :: file_var_meta_DD + +!> type to hold metadata for variables in a non-domain-decomposed file +type (file_variable_meta_noDD), private :: file_var_meta_noDD + +!> index of the time_level value that is written to netCDF file bythe write_field routines. +integer, private :: write_field_time_index + +!> interface to apply a scale factor to an array after reading in a field +interface scale_data + module procedure scale_data_4d + module procedure scale_data_3d + module procedure scale_data_2d + module procedure scale_data_1d +end interface + +contains + +!> This routine calls the fms_io read_data subroutine to read 1-D domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_1d_DD(filename, fieldname, data, domain, use_fms2, start_index, edge_lengths, & + timelevel, scale, x_position, y_position, leave_file_open) + 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 1-dimensional data array to pass to read_data + type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + logical, intent(in) :: use_fms2 !< flag to distinguish interface from the old write_field interface + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in + !! default is the variable size + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by + integer, optional, intent(in) :: x_position !< domain position of x-dimension; CENTER (default) or EAST_FACE + integer, optional, intent(in) :: y_position !< domain position of y-dimension; CENTER (default) or NORTH_FACE + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, num_var_dims, dim_unlim_size + integer, dimension(1) :: start, nread ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + integer :: xpos, ypos ! x and y domain positions + + xpos = CENTER + ypos = CENTER + if (present(x_position)) xpos = x_position + if (present(y_position)) ypos = y_position + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + ! define the io domain for 1-pe jobs because it is required to read domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", domain%mpp_domain, is_restart=.false.) + file_var_meta_DD%nvars = get_num_variables(fileobj_read_dd) + if (file_var_meta_DD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_DD%var_names))) allocate(file_var_meta_DD%var_names(file_var_meta_DD%nvars)) + call get_variable_names(fileobj_read_dd, file_var_meta_DD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_DD%nvars + if (lowercase(trim(file_var_meta_DD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_DD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_1d_DD: "//& + trim(fieldname)//" not found in"//trim(filename)) + ! register the variable axes + call MOM_register_variable_axes(fileobj_read_dd, trim(variable_to_read), xPosition=xpos, yPosition=ypos) + ! set the start and nread values that will be passed as the read_data start_index and edge_lengths arguments + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read_dd, trim(variable_to_read), dim_names) + + start(1)=1 + if (present(timelevel)) then + if (is_dimension_unlimited(fileobj_read_dd, dim_names(1))) start(1) = timelevel + elseif (present(start_index)) then + start(1) = start_index(1) + endif + + if (present(edge_lengths)) then + nread(1) = edge_lengths(1) + else + call get_dimension_size(fileobj_read_dd, trim(dim_names(1)), nread(1)) + endif + ! read the data + dim_unlim_size = 0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read_dd, dim_names(i))) then + call get_dimension_size(fileobj_read_dd, dim_names(i), dim_unlim_size) + exit + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call MOM_error(WARNING, "MOM_read_data_fms2::MOM_read_data_1d_DD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + if (allocated(file_var_meta_DD%var_names)) deallocate(file_var_meta_DD%var_names) + file_var_meta_DD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_read_data_1d_DD + +!> This routine calls the fms_io read_data subroutine to read 2-D domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_2d_DD(filename, fieldname, data, domain, use_fms2, start_index, edge_lengths, & + timelevel, scale, x_position, y_position, leave_file_open) + 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 2-dimensional data array to pass to read_data + type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + logical, intent(in) :: use_fms2 !< flag to distinguish interface from the old write_field interface + integer, dimension(2), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + integer, intent(in), optional :: x_position !< domain position of x-dimension; CENTER (default) or EAST_FACE + integer, intent(in), optional :: y_position !< domain position of y-dimension; CENTER (default) or NORTH_FACE + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims, first(2), last(2) + integer :: start(2), nread(2) ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + integer :: xpos, ypos, pos ! x and y domain positions + integer :: isc, iec, jsc, jec, isg, ieg, jsg, jeg + type(domain2D), pointer :: io_domain => NULL() + + xpos = CENTER + ypos = CENTER + if (present(x_position)) xpos = x_position + if (present(y_position)) ypos = y_position + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + ! define the io domain for 1-pe jobs because it is required to read domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", domain%mpp_domain, is_restart=.false.) + file_var_meta_DD%nvars = get_num_variables(fileobj_read_dd) + if (file_var_meta_DD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_DD%var_names))) allocate(file_var_meta_DD%var_names(file_var_meta_DD%nvars)) + call get_variable_names(fileobj_read_dd, file_var_meta_DD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_DD%nvars + if (lowercase(trim(file_var_meta_DD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_DD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_2d_DD: "//& + trim(fieldname)//" not found in "//trim(filename)) + ! register the variable axes + call MOM_register_variable_axes(fileobj_read_dd, trim(variable_to_read), xPosition=xpos, yPosition=ypos) + + pos = CENTER + if (present(x_position)) then + if (present(y_position)) then + pos = CORNER + else + pos = xpos + endif + elseif (present(y_position)) then + pos = ypos + endif + ! set the start and nread values that will be passed as the read_data corner and edge_lengths argument + num_var_dims = get_variable_num_dimensions(fileobj_read_dd, trim(variable_to_read)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read_dd, trim(variable_to_read), dim_names) + ! get the IO domain + !io_domain => mpp_get_io_domain(domain%mpp_domain) + ! Get the global indicies + !call mpp_get_global_domain(io_domain, xbegin=isg, xend=ieg, ybegin=jsg, yend=jeg, position=pos) + ! Get the compute indicies + !call mpp_get_compute_domain(io_domain, xbegin=isc, xend=iec, ybegin=jsc, yend=jec, position=pos) + !last(1) = iec - isg + 1 ! get array indices for the axis data + !last(2) = jec - jsg + 1 + !first(1) = isc - isg + 1 + !first(2) = jsc - jsg + 1 + + start(:) = 1 + if (present(start_index)) then + start = start_index + !else + ! start(:) = first(:) + endif + + if (present(edge_lengths)) then + nread = edge_lengths + else + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read_dd, dim_names(i))) then + call get_dimension_size(fileobj_read_dd, dim_names(i), dim_unlim_size) + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + if (allocated(file_var_meta_DD%var_names)) deallocate(file_var_meta_DD%var_names) + file_var_meta_DD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) + if (associated(io_domain)) nullify(io_domain) +end subroutine MOM_read_data_2d_DD + +!> This routine calls the fms_io read_data subroutine to read 3-D domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_3d_DD(filename, fieldname, data, domain, use_fms2, start_index, edge_lengths, & + timelevel, scale, x_position, y_position, leave_file_open) + 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 data array to pass to read_data + type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + logical, intent(in) :: use_fms2 !< flag indicating whether to use this routine + integer, dimension(3), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(3), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + integer, intent(in), optional :: x_position !< domain position of x-dimension; CENTER (default) or EAST_FACE + integer, intent(in), optional :: y_position !< domain position of y-dimension; CENTER (default) or NORTH_FACE + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! if .true., the variable was found in the netCDF file + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims + integer, dimension(3) :: start, nread, first, last ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + integer :: xpos, ypos, pos ! x and y domain positions + integer :: isc, iec, jsc, jec, isg, ieg, jsg, jeg + type(domain2D), pointer :: io_domain => NULL() + + xpos = CENTER + ypos = CENTER + if (present(x_position)) xpos = x_position + if (present(y_position)) ypos = y_position + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + ! define the io domain for 1-pe jobs because it is required to read domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", domain%mpp_domain, is_restart=.false.) + file_var_meta_DD%nvars = get_num_variables(fileobj_read_dd) + if (file_var_meta_DD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_DD%var_names))) allocate(file_var_meta_DD%var_names(file_var_meta_DD%nvars)) + call get_variable_names(fileobj_read_dd, file_var_meta_DD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_DD%nvars + if (lowercase(trim(file_var_meta_DD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_DD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_3d_DD: "//& + trim(fieldname)//" not found in"//trim(filename)) + ! register the variable axes + call MOM_register_variable_axes(fileobj_read_dd, trim(variable_to_read), xPosition=xpos, yPosition=ypos) + pos = CENTER + if (present(x_position)) then + if (present(y_position)) then + pos = CORNER + else + pos = xpos + endif + elseif (present(y_position)) then + pos = ypos + endif + ! set the start and nread values that will be passed as the read_data corner and edge_lengths argument + num_var_dims = get_variable_num_dimensions(fileobj_read_dd, trim(variable_to_read)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read_dd, trim(variable_to_read), dim_names) + ! get the IO domain + io_domain => mpp_get_io_domain(domain%mpp_domain) + ! Get the global indicies + ! call mpp_get_global_domain(io_domain, xbegin=isg, xend=ieg, ybegin=jsg, yend=jeg, position=pos) + ! call mpp_get_compute_domain(io_domain, xbegin=isc, xend=iec, ybegin=jsc, yend=jec, position=pos) + !last(1) = iec - isg + 1 ! get array indices for the axis data + !last(2) = jec - jsg + 1 + !first(1) = isc - isg + 1 + !first(2) = jsc - jsg + 1 + + start(:) = 1 + if (present(start_index)) then + start = start_index + !else + ! start(1:2) = first(1:2) + endif + + if (present(edge_lengths)) then + nread = edge_lengths + else + !nread(1) = last(1) - first(1) + 1 + !nread(2) = last(2) - first(2) + 1 + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read_dd, dim_names(i))) then + call get_dimension_size(fileobj_read_dd, dim_names(i), dim_unlim_size) + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call MOM_error(WARNING, "MOM_read_data_fms2::MOM_read_data_3d_DD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + if (allocated(file_var_meta_DD%var_names)) deallocate(file_var_meta_DD%var_names) + file_var_meta_DD%nvars = 0 + endif + + if (allocated(dim_names)) deallocate(dim_names) + if (associated(io_domain)) nullify(io_domain) +end subroutine MOM_read_data_3d_DD + +!> This routine calls the fms_io read_data subroutine to read 4-D domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_4d_DD(filename, fieldname, data, domain, use_fms2, start_index, edge_lengths, & + timelevel, scale, x_position, y_position, leave_file_open) + 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 1-dimensional data array to pass to read_data + type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + logical, intent(in) :: use_fms2 !< flag indicating whether to use this routine + integer, dimension(4), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + integer, intent(in), optional :: x_position !< domain position of x-dimension; CENTER (default) or EAST_FACE + integer, intent(in), optional :: y_position !< domain position of y-dimension; CENTER (default) or NORTH_FACE + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims + integer, dimension(4) :: start, nread, first, last ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + integer :: xpos, ypos, pos ! x and y domain positions + integer :: isc, iec, jsc, jec, isg, ieg, jsg, jeg + type(domain2D), pointer :: io_domain => NULL() + + xpos = CENTER + ypos = CENTER + if (present(x_position)) xpos = x_position + if (present(y_position)) ypos = y_position + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + ! define the io domain for 1-pe jobs because it is required to read domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", domain%mpp_domain, is_restart=.false.) + file_var_meta_DD%nvars = get_num_variables(fileobj_read_dd) + if (file_var_meta_DD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_DD%var_names))) allocate(file_var_meta_DD%var_names(file_var_meta_DD%nvars)) + call get_variable_names(fileobj_read_dd, file_var_meta_DD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_DD%nvars + if (lowercase(trim(file_var_meta_DD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_DD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) & + call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_4d_DD: "//trim(fieldname)//" not found in"//& + trim(filename)) + ! register the variable axes + call MOM_register_variable_axes(fileobj_read_dd, trim(variable_to_read), xPosition=xpos, yPosition=ypos) + pos = CENTER + if (present(x_position)) then + if (present(y_position)) then + pos = CORNER + else + pos = xpos + endif + elseif (present(y_position)) then + pos = ypos + endif + ! set the start and nread values that will be passed as the read_data corner and edge_lengths argument + num_var_dims = get_variable_num_dimensions(fileobj_read_dd, trim(variable_to_read)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read_dd, trim(variable_to_read), dim_names) + ! get the IO domain + !io_domain => mpp_get_io_domain(domain%mpp_domain) + ! Get the global indicies + !call mpp_get_global_domain(domain%mpp_domain, xbegin=isg, xend=ieg, ybegin=jsg, yend=jeg, position=pos) + ! Get the compute indicies + ! call mpp_get_compute_domain(domain%mpp_domain, xbegin=isc, xend=iec, ybegin=jsc, yend=jec, position=pos) + !last(1) = iec - isg + 1 ! get array indices for the axis data + !first(1) = isc - isg + 1 + + start(:) = 1 + if (present(start_index)) then + start(:) = start_index(:) + !else + !start(1:2) = first(1:2) + endif + + if (present(edge_lengths)) then + nread = edge_lengths + else + !nread(1) = last(1) - first(1) + 1 + !nread(2) = last(2) - first(2) + 1 + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1, num_var_dims + if (is_dimension_unlimited(fileobj_read_dd, dim_names(i))) then + call get_dimension_size(fileobj_read_dd, dim_names(i), dim_unlim_size) + endif + if (i .eq. 4) then + nread(i) = 1 + start(i) = timelevel + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call MOM_error(WARNING, "MOM_read_data_fms2::MOM_read_data_4d_DD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read_dd, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + if (allocated(file_var_meta_DD%var_names)) deallocate(file_var_meta_DD%var_names) + file_var_meta_DD%nvars = 0 + endif + if (associated(io_domain)) nullify(io_domain) + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_read_data_4d_DD + +!!> This routine calls the fms_io read_data subroutine to read a scalar (0-D) field named "fieldname" +!! from file "filename". +subroutine MOM_read_data_scalar(filename, fieldname, data, use_fms2, leave_file_open) + 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, intent(inout) :: data !< data buffer to pass to read_data + logical, intent(in) :: use_fms2 !< flag distinguishing interface from old MOM_read_data + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + integer :: i + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read))) then + file_open_success = fms2_open_file(fileobj_read, filename, "read", is_restart=.false.) + file_var_meta_noDD%nvars = get_num_variables(fileobj_read) + if (file_var_meta_noDD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_noDD%var_names))) & + allocate(file_var_meta_noDD%var_names(file_var_meta_noDD%nvars)) + call get_variable_names(fileobj_read, file_var_meta_noDD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_noDD%nvars + if (lowercase(trim(file_var_meta_noDD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_noDD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_scalar: "//trim(fieldname)// & + " not found in"//trim(filename)) + ! read the data + call read_data(fileobj_read, trim(fieldname), data) + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read)) call fms2_close_file(fileobj_read) + if (allocated(file_var_meta_noDD%var_names)) deallocate(file_var_meta_noDD%var_names) + file_var_meta_noDD%nvars = 0 + endif +end subroutine MOM_read_data_scalar + +!> This routine calls the fms_io read_data subroutine to read 1-D non-domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_1d_noDD(filename, fieldname, data, use_fms2, start_index, & + edge_lengths, timelevel, scale, leave_file_open) + 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 1-dimensional data array to pass to read_data + logical, intent(in) :: use_fms2 !< flag to distinguish interface from old MOM_read_data interface + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is + !! the variable size + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + integer :: i, num_var_dims, dim_unlim_size + integer, dimension(1) :: start, nread ! indices for first data value and number of values to read + character(len=40), allocatable:: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read))) then + file_open_success = fms2_open_file(fileobj_read, filename, "read", is_restart=.false.) + file_var_meta_noDD%nvars = get_num_variables(fileobj_read) + if (file_var_meta_noDD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_noDD%var_names))) & + allocate(file_var_meta_noDD%var_names(file_var_meta_noDD%nvars)) + call get_variable_names(fileobj_read, file_var_meta_noDD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_noDD%nvars + if (lowercase(trim(file_var_meta_noDD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_noDD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_io:MOM_read_data_1d_noDD: "//trim(fieldname)//& + " not found in "//trim(filename)) + + num_var_dims = get_variable_num_dimensions(fileobj_read, trim(fieldname)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read, trim(variable_to_read), dim_names) + + ! set the start and nread values that will be passed as the read_data start_index and edge_lengths arguments + start(1)=1 + if (present(timelevel)) then + if (is_dimension_unlimited(fileobj_read, dim_names(1))) start(1) = timelevel + elseif (present(start_index)) then + start(1) = start_index(1) + endif + + if (present(edge_lengths)) then + nread(1) = edge_lengths(1) + else + nread = shape(data) + endif + ! read the data + dim_unlim_size = 0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read, dim_names(i))) then + call get_dimension_size(fileobj_read, dim_names(i), dim_unlim_size) + exit + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call MOM_error(WARNING, "MOM_read_data_fms2::MOM_read_data_1d_noDD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read)) call fms2_close_file(fileobj_read) + if (allocated(file_var_meta_noDD%var_names)) deallocate(file_var_meta_noDD%var_names) + file_var_meta_noDD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_read_data_1d_noDD + +!> This routine calls the fms_io read_data subroutine to read 2-D non-domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_2d_noDD(filename, fieldname, data, use_fms2, start_index, & + edge_lengths, timelevel, scale, leave_file_open) + 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 2-dimensional data array to pass to read_data + logical, intent(in) :: use_fms2 !< flag to distinguish interface from old MOM_read_data interface + integer, dimension(2), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims + integer, dimension(2) :: start, nread ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read))) then + file_open_success = fms2_open_file(fileobj_read, filename, "read", is_restart=.false.) + file_var_meta_noDD%nvars = get_num_variables(fileobj_read) + if (file_var_meta_noDD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_noDD%var_names))) & + allocate(file_var_meta_noDD%var_names(file_var_meta_noDD%nvars)) + call get_variable_names(fileobj_read, file_var_meta_noDD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_noDD%nvars + if (lowercase(trim(file_var_meta_noDD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_noDD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) & + call MOM_error(FATAL, "MOM_io:MOM_read_data_2d_noDD: "//trim(fieldname)//& + " not found in "//trim(filename)) + ! set the start and nread values that will be passed as the read_data corner and edge_lengths arguments + start(:) = 1 + if (present(start_index)) start = start_index + + if (present(edge_lengths)) then + nread = edge_lengths + else + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + num_var_dims = get_variable_num_dimensions(fileobj_read, trim(fieldname)) + allocate(dim_names(num_var_dims)) + call get_variable_dimension_names(fileobj_read, trim(variable_to_read), dim_names) + dim_names(:) = "" + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read, dim_names(i))) then + call get_dimension_size(fileobj_read, dim_names(i), dim_unlim_size) + endif + enddo + if (dim_unlim_size .LE. 0) then + call MOM_error(WARNING, "MOM_io::MOM_read_data_2d_noDD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + endif + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read)) call fms2_close_file(fileobj_read) + if (allocated(file_var_meta_noDD%var_names)) deallocate(file_var_meta_noDD%var_names) + file_var_meta_noDD%nvars = 0 + endif + if(allocated(dim_names)) deallocate(dim_names) + +end subroutine MOM_read_data_2d_noDD + +!> This routine calls the fms_io read_data subroutine to read 3-D non-domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_3d_noDD(filename, fieldname, data, use_fms2, start_index, & + edge_lengths, timelevel, scale, leave_file_open) + 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 data array to pass to read_data + logical, intent(in) :: use_fms2 !< flag to distinguish interface from old MOM_read_data interface + integer, dimension(3), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(3), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims + integer, dimension(3) :: start, nread ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read))) then + file_open_success = fms2_open_file(fileobj_read, filename, "read", is_restart=.false.) + file_var_meta_noDD%nvars = get_num_variables(fileobj_read) + if (file_var_meta_noDD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_noDD%var_names))) & + allocate(file_var_meta_noDD%var_names(file_var_meta_noDD%nvars)) + call get_variable_names(fileobj_read, file_var_meta_noDD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_noDD%nvars + if (lowercase(trim(file_var_meta_noDD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_noDD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_io:MOM_read_data_3d_noDD: "//trim(fieldname)//& + " not found in "//trim(filename)) + ! get the variable dimensions + num_var_dims = get_variable_num_dimensions(fileobj_read, trim(fieldname)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read, trim(variable_to_read), dim_names) + ! set the start and nread values that will be passed as the read_data corner and edge_lengths arguments + start(:) = 1 + if (present(start_index)) start = start_index + + if (present(edge_lengths)) then + nread = edge_lengths + else + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read, dim_names(i))) then + call get_dimension_size(fileobj_read, dim_names(i), dim_unlim_size) + endif + enddo + if (dim_unlim_size .LE. 0) then + call MOM_error(WARNING, "MOM_read_data_fms2::MOM_read_data_3d_noDD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + endif + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read)) call fms2_close_file(fileobj_read) + if (allocated(file_var_meta_noDD%var_names)) deallocate(file_var_meta_noDD%var_names) + file_var_meta_noDD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_read_data_3d_noDD + +!> This routine calls the fms_io read_data subroutine to read 4-D non-domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +subroutine MOM_read_data_4d_noDD(filename, fieldname, data, use_fms2, start_index, & + edge_lengths, timelevel, scale, leave_file_open) + 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 4-dimensional data array to pass to read_data + logical, intent(in) :: use_fms2 !< flag to distinguish interface from old MOM_read_data interface + integer, dimension(4), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after read_data is called; default is .true. + integer :: i, dim_unlim_size, num_var_dims + integer, dimension(4) :: start, nread ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + ! open the file + if (.not.(check_if_open(fileobj_read))) then + file_open_success = fms2_open_file(fileobj_read, filename, "read", is_restart=.false.) + file_var_meta_noDD%nvars = get_num_variables(fileobj_read) + if (file_var_meta_noDD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_noDD%var_names))) & + allocate(file_var_meta_noDD%var_names(file_var_meta_noDD%nvars)) + call get_variable_names(fileobj_read, file_var_meta_noDD%var_names) + endif + ! search for the variable in the file + variable_to_read = "" + variable_found = .false. + do i=1,file_var_meta_noDD%nvars + if (lowercase(trim(file_var_meta_noDD%var_names(i))) .eq. lowercase(trim(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_noDD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(FATAL, "MOM_read_data_fms2:MOM_read_data_4d_noDD: "//& + trim(fieldname)//" not found in "//trim(filename)) + ! get the variable dimensions + num_var_dims = get_variable_num_dimensions(fileobj_read, trim(fieldname)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read, trim(variable_to_read), dim_names) + ! set the start and nread values that will be passed as the read_data corner and edge_lengths arguments + start(:) = 1 + if (present(start_index)) start = start_index + + if (present(edge_lengths)) then + nread = edge_lengths + else + nread = shape(data) + endif + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1, num_var_dims + if (is_dimension_unlimited(fileobj_read, dim_names(i))) then + call get_dimension_size(fileobj_read, dim_names(i), dim_unlim_size) + endif + if (i .eq. 4) then + nread(i) = 1 + start(i) = timelevel + endif + enddo + if (dim_unlim_size .LE. 0) then + call MOM_error(WARNING, "MOM_io::MOM_read_data_4d_noDD: time level specified, but the variable "//& + trim(fieldName)// " does not have an unlimited dimension.") + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + endif + else + call read_data(fileobj_read, trim(variable_to_read), data, corner=start, edge_lengths=nread) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read)) call fms2_close_file(fileobj_read) + if (allocated(file_var_meta_noDD%var_names)) deallocate(file_var_meta_noDD%var_names) + file_var_meta_noDD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) +end subroutine MOM_read_data_4d_noDD + +!> This routine calls the fms_io read_data subroutine to read 2-D domain-decomposed data field named "fieldname" +!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. +!!The supergrid variable axis lengths are determined from compute domain lengths, and +!! the domain indices are computed from the difference between the global and compute domain indices +subroutine MOM_read_data_2d_supergrid(filename, fieldname, data, domain, is_supergrid, start_index, edge_lengths, & + timelevel, scale, x_position, y_position, leave_file_open) + 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 2-dimensional data array to pass to read_data + type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + logical, intent(in) :: is_supergrid !< flag indicating whether to use supergrid + integer, dimension(2), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 + integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in. + !! Default values are the variable dimension sizes + integer, optional, intent(in) :: timelevel !< time level to read + real, optional, intent(in):: scale !< A scaling factor that the field is multiplied by + integer, intent(in), optional :: x_position !< domain position of x-dimension; CENTER (default) or EAST_FACE + integer, intent(in), optional :: y_position !< domain position of y-dimension; CENTER (default) or NORTH_FACE + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: variable_found ! .true. if lowercase(fieldname) matches one of the lowercase file variable names + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + integer :: i, dim_unlim_size, npes, num_var_dims, first(2), last(2) + integer :: start(2), nread(2) ! indices for first data value and number of values to read + character(len=40), allocatable :: dim_names(:) ! variable dimension names + character(len=96) :: variable_to_read ! variable to read from the netcdf file + integer :: xpos, ypos, pos ! x and y domain positions + integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, isg, ieg, jsg, jeg + integer :: xsize_c, ysize_c, xsize_d, ysize_d + real, allocatable :: array(:,:) ! dummy array to pass to read data + type(domain2D), pointer :: io_domain => NULL() + + xpos = CENTER + ypos = CENTER + if (present(x_position)) xpos = x_position + if (present(y_position)) ypos = y_position + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + npes=-1; npes = mpp_get_domain_npes(domain%mpp_domain) + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + ! define the io domain for 1-pe jobs because it is required to read domain-decomposed files + if (npes .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", domain%mpp_domain, is_restart=.false.) + file_var_meta_DD%nvars = get_num_variables(fileobj_read_dd) + if (file_var_meta_DD%nvars .lt. 1) call MOM_error(FATAL, "nvars is less than 1 for file "// & + trim(filename)) + if (.not.(allocated(file_var_meta_DD%var_names))) & + allocate(file_var_meta_DD%var_names(file_var_meta_DD%nvars)) + call get_variable_names(fileobj_read_dd, file_var_meta_DD%var_names) + endif + ! search for the variable in the file + variable_to_read = trim(fieldname) + variable_found = .false. + do i=1,file_var_meta_DD%nvars + if (trim(lowercase(file_var_meta_DD%var_names(i))) .eq. trim(lowercase(fieldname))) then + variable_found = .true. + variable_to_read = trim(file_var_meta_DD%var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) call MOM_error(WARNING, "MOM_read_data_fms2:MOM_read_data_2d_supergrid: "//& + trim(fieldname)//" not found in "//trim(filename)) + + pos = CENTER + if (xpos .eq. NORTH_FACE) then + if (ypos .eq. EAST_FACE) then + pos = CORNER + else + pos = xpos + endif + elseif (ypos .eq. EAST_FACE) then + pos = ypos + endif + ! set the start and nread values that will be passed as the read_data corner and edge_lengths argument + num_var_dims = get_variable_num_dimensions(fileobj_read_dd, trim(variable_to_read)) + allocate(dim_names(num_var_dims)) + dim_names(:) = "" + call get_variable_dimension_names(fileobj_read_dd, trim(variable_to_read), dim_names) + ! get the IO domain + io_domain => mpp_get_io_domain(domain%mpp_domain) + ! register the variable axes + !call MOM_register_variable_axes(fileobj_read, trim(variable_to_read), io_domain, xPosition=xpos, yPosition=ypos) + call mpp_get_data_domain(domain%mpp_domain,isd,ied,jsd,jed,xsize=xsize_d,ysize=ysize_d,position=pos) + call mpp_get_global_domain(domain%mpp_domain,isg,ieg,jsg,jeg,position=pos) + call mpp_get_compute_domain(domain%mpp_domain,isc,iec,jsc,jec,position=pos) + ! get the start indices + start(:) = 1 + if (present(start_index)) then + start = start_index + else!if((size(data,1) .eq. xsize_d) .and. (size(data,2) .eq. ysize_d)) then ! on_data_domain + if (npes .gt. 1) then + start(1) = isc - isg + 1 + start(2) = jsc - jsg + 1 + else + if (iec-isc+1 .ne. ieg-isg+1) start(1) = isc - isg + 1 + if (jec-jsc+1 .ne. jeg-jsg+1) start(2) = jsc - jsg + 1 + endif + endif + ! get the values for the edge_lengths (nread) + nread = shape(data) + if (present(edge_lengths)) then + nread = edge_lengths + else!if((size(data,1) .eq. xsize_d) .and. (size(data,2) .eq. ysize_d)) then ! on_data_domain + if (npes .gt. 1) then + nread(1) = iec - isc + 1 + nread(2) = jec - jsc + 1 + else + if (iec-isc+1 .ne. ieg-isg+1) nread(1) = iec - isc + 1 + if (jec-jsc+1 .ne. jeg-jsg+1) nread(2) = jec - jsc + 1 + endif + endif + ! allocate the dummy array + if (.not. allocated(array)) allocate(array(size(data,1),size(data,2))) + ! read the data + dim_unlim_size=0 + if (present(timelevel)) then + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj_read_dd, dim_names(i))) then + call get_dimension_size(fileobj_read_dd, dim_names(i), dim_unlim_size) + endif + enddo + if (dim_unlim_size .gt. 0) then + call read_data(fileobj_read_dd, trim(variable_to_read), array, corner=start, edge_lengths=nread, & + unlim_dim_level=timelevel) + else + call read_data(fileobj_read_dd, trim(variable_to_read), array, corner=start, edge_lengths=nread) + endif + else + call read_data(fileobj_read_dd, trim(variable_to_read), array, corner=start, edge_lengths=nread) + endif + if((size(array,1) .eq. xsize_d) .and. (size(array,2) .eq. ysize_d)) then ! on_data_domain + data(isc-isd+1:iec-isd+1,jsc-jsd+1:jec-jsd+1) = array(isc-isd+1:iec-isd+1,jsc-jsd+1:jec-jsd+1) + else + data = array + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data, scale) + endif ; endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + if (allocated(file_var_meta_DD%var_names)) deallocate(file_var_meta_DD%var_names) + file_var_meta_noDD%nvars = 0 + endif + if (allocated(dim_names)) deallocate(dim_names) + if (associated(io_domain)) nullify(io_domain) + if (allocated(array)) deallocate(array) +end subroutine MOM_read_data_2d_supergrid + + +!> This routine uses the fms2_io read_data interface to read a pair of distributed +!! 2-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for +!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. +subroutine MOM_read_vector_2d_fms2(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & + use_fms2, timelevel, stagger, scale, leave_file_open) + character(len=*), intent(in) :: filename !< name of the netcdf file to read + character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file + character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file + real, dimension(:,:), intent(inout) :: u_data !< The 2 dimensional array into which the + !! u-component of the data should be read + real, dimension(:,:), intent(inout) :: v_data !< The 2 dimensional array into which the + !! v-component of the data should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + logical, intent(in) :: use_fms2 !< flag indicating whether to use this routine + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied + !! by before they are returned. + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + integer :: is, ie, js, je, i, ndims, dim_unlim_index + integer :: u_pos, v_pos + integer, allocatable :: dim_sizes_u(:), dim_sizes_v(:) + character(len=32), allocatable :: dim_names_u(:), dim_names_v(:), units_u(:), units_v(:) + character(len=1) :: x_or_y ! orientation of cartesian coordinate axis + logical :: is_valid + logical :: file_open_success ! .true. if open file is successful + logical :: close_the_file ! indicates whether to close the file after MOM_read_vector is called; default is .true. + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) & + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + if (.not. file_open_success) call MOM_error(FATAL, "MOM_read_vector_2d_fms2: netcdf file "//& + trim(filename)//" not opened.") + + u_pos = EAST_FACE ; v_pos = NORTH_FACE + if (present(stagger)) then + if (stagger == CGRID_NE .or. stagger == BGRID_NE ) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE + elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif + endif + + ndims = get_variable_num_dimensions(fileobj_read_dd, u_fieldname) + allocate(dim_sizes_u(ndims)) + allocate(dim_sizes_v(ndims)) + allocate(dim_names_u(ndims)) + allocate(dim_names_v(ndims)) + allocate(units_u(ndims)) + allocate(units_v(ndims)) + + units_u(:) = "" + units_v(:) = "" + dim_names_u(:) = "" + dim_names_v(:) = "" + dim_sizes_u(:) = 0 + dim_sizes_v(:) = 0 + + call get_variable_size(fileobj_read_dd, u_fieldname, dim_sizes_u) + call get_variable_size(fileobj_read_dd, v_fieldname, dim_sizes_v) + call get_variable_dimension_names(fileobj_read_dd, u_fieldname, dim_names_u) + call get_variable_dimension_names(fileobj_read_dd, v_fieldname, dim_names_v) + + do i=1,ndims + ! register the u axes + if (.not.(is_dimension_registered(fileobj_read_dd, dim_names_u(i)))) then + call get_variable_units(fileobj_read_dd, dim_names_u(i), units_u(i)) + call validate_lat_lon_units(units_u(i), x_or_y, is_valid) + if (is_valid) then + call register_axis(fileobj_read_dd, dim_names_u(i), x_or_y, domain_position=u_pos) + else + call register_axis(fileobj_read_dd, dim_names_u(i), dim_sizes_u(i)) + endif + endif + ! Register the v axes if they differ from the u axes + if (trim(lowercase(dim_names_v(i))) .ne. trim(lowercase(dim_names_u(i)))) then + if (.not.(is_dimension_registered(fileobj_read_dd, dim_names_v(i)))) then + call get_variable_units(fileobj_read_dd, dim_names_v(i), units_v(i)) + call validate_lat_lon_units(units_v(i), x_or_y, is_valid) + if (is_valid) then + call register_axis(fileobj_read_dd, dim_names_v(i), x_or_y, domain_position=v_pos) + else + call register_axis(fileobj_read_dd, dim_names_v(i), dim_sizes_v(i)) + endif + endif + endif + enddo + ! read the data + dim_unlim_index = 0 + if (present(timelevel)) then + do i=1,ndims + if (is_dimension_unlimited(fileobj_read_dd, dim_names_u(i))) then + dim_unlim_index = i + exit + endif + enddo + if (dim_unlim_index .gt. 0) then + call read_data(fileobj_read_dd, u_fieldname,u_data, unlim_dim_level=timelevel) + call read_data(fileobj_read_dd, v_fieldname, v_data, unlim_dim_level=timelevel) + else + call read_data(fileobj_read_dd, u_fieldname, u_data) + call read_data(fileobj_read_dd, v_fieldname, v_data) + endif + else + call read_data(fileobj_read_dd, u_fieldname, u_data) + call read_data(fileobj_read_dd, v_fieldname, v_data) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) + u_data(is:ie,js:je) = scale*u_data(is:ie,js:je) + call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) + v_data(is:ie,js:je) = scale*v_data(is:ie,js:je) + endif ; endif + if (allocated(dim_names_u)) deallocate(dim_names_u) + if (allocated(dim_names_v)) deallocate(dim_names_v) + if (allocated(dim_sizes_u)) deallocate(dim_sizes_u) + if (allocated(dim_sizes_v)) deallocate(dim_sizes_v) + if (allocated(units_u)) deallocate(units_u) + if (allocated(units_v)) deallocate(units_v) +end subroutine MOM_read_vector_2d_fms2 + +!> This routine uses the fms2_io read_data interface to read a pair of distributed +!! 3-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for +!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. +subroutine MOM_read_vector_3d_fms2(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & + use_fms2, timelevel, stagger, scale, leave_file_open) + character(len=*), intent(in) :: filename !< name of the netcdf file to read + character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file + character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file + real, dimension(:,:,:), intent(inout) :: u_data !< The 3 dimensional array into which the + !! u-component of the data should be read + real, dimension(:,:,:), intent(inout) :: v_data !< The 3 dimensional array into which the + !! v-component of the data should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + logical, intent(in) :: use_fms2 !< flag indicating whether to call this routine + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied + !! by before they are returned. + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + ! local + integer :: is, ie, js, je, i, dim_unlim, ndims + integer :: u_pos, v_pos + integer, allocatable :: dim_sizes_u(:), dim_sizes_v(:) + character(len=32), allocatable :: dim_names_u(:), dim_names_v(:), units_u(:), units_v(:) + character(len=1) :: x_or_y + logical :: is_valid + logical :: file_open_success ! .true. if open file is successful + logical :: close_the_file ! indicates whether to close the file after MOM_read_vector is called; default is .true. + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + ! open the file + if (.not.(check_if_open(fileobj_read_dd))) then + file_open_success = fms2_open_file(fileobj_read_dd, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + if (.not. file_open_success) & + call MOM_error(FATAL, "MOM_read_vector_3d_fms2: netcdf file "//trim(filename)//" not opened.") + endif + + u_pos = EAST_FACE ; v_pos = NORTH_FACE + if (present(stagger)) then + if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE + elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER + elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif + endif + + ndims = get_variable_num_dimensions(fileobj_read_dd, u_fieldname) + allocate(dim_sizes_u(ndims)) + allocate(dim_sizes_v(ndims)) + allocate(dim_names_u(ndims)) + allocate(dim_names_v(ndims)) + allocate(units_u(ndims)) + allocate(units_v(ndims)) + + units_u(:) = "" + units_v(:) = "" + dim_names_u(:) = "" + dim_names_v(:) = "" + + call get_variable_size(fileobj_read_dd, u_fieldname, dim_sizes_u, broadcast=.true.) + call get_variable_size(fileobj_read_dd, v_fieldname, dim_sizes_v, broadcast=.true.) + call get_variable_dimension_names(fileobj_read_dd, u_fieldname, dim_names_u, broadcast=.true.) + call get_variable_dimension_names(fileobj_read_dd, v_fieldname, dim_names_v, broadcast=.true.) + + do i=1,ndims + ! register the u axes + if (.not.(is_dimension_registered(fileobj_read_dd, dim_names_u(i)))) then + call get_variable_units(fileobj_read_dd, dim_names_u(i), units_u(i)) + call validate_lat_lon_units(units_u(i), x_or_y, is_valid) + if (is_valid) then + call register_axis(fileobj_read_dd, dim_names_u(i), x_or_y, domain_position=u_pos) + else + call register_axis(fileobj_read_dd, dim_names_u(i), dim_sizes_u(i)) + endif + endif + ! Register the v axes if they differ from the u axes + if (trim(lowercase(dim_names_v(i))) .ne. trim(lowercase(dim_names_u(i)))) then + if (.not.(is_dimension_registered(fileobj_read_dd, dim_names_v(i)))) then + call get_variable_units(fileobj_read_dd, dim_names_v(i), units_v(i)) + call validate_lat_lon_units(units_v(i), x_or_y, is_valid) + if (is_valid) then + call register_axis(fileobj_read_dd, dim_names_v(i), x_or_y, domain_position=v_pos) + else + call register_axis(fileobj_read_dd, dim_names_v(i), dim_sizes_v(i)) + endif + endif + endif + enddo + ! read the data + dim_unlim = 0 + if (present(timelevel)) then + do i=1,ndims + if (is_dimension_unlimited(fileobj_read_dd, dim_names_u(i))) then + dim_unlim = i + exit + endif + enddo + if (dim_unlim .gt. 0) then + call read_data(fileobj_read_dd, u_fieldname, u_data, unlim_dim_level=timelevel) + call read_data(fileobj_read_dd, v_fieldname, v_data, unlim_dim_level=timelevel) + else + call read_data(fileobj_read_dd, u_fieldname, u_data, edge_lengths=dim_sizes_u) + call read_data(fileobj_read_dd, v_fieldname, v_data, edge_lengths=dim_sizes_v) + endif + else + call read_data(fileobj_read_dd, u_fieldname, u_data, edge_lengths=dim_sizes_u) + call read_data(fileobj_read_dd, v_fieldname, v_data, edge_lengths=dim_sizes_v) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_read_dd)) call fms2_close_file(fileobj_read_dd) + endif + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) + u_data(is:ie,js:je,:) = scale*u_data(is:ie,js:je,:) + call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) + v_data(is:ie,js:je,:) = scale*v_data(is:ie,js:je,:) + endif ; endif + if (allocated(dim_names_u)) deallocate(dim_names_u) + if (allocated(dim_names_v)) deallocate(dim_names_v) + if (allocated(dim_sizes_u)) deallocate(dim_sizes_u) + if (allocated(dim_sizes_v)) deallocate(dim_sizes_v) + if (allocated(units_u)) deallocate(units_u) + if (allocated(units_v)) deallocate(units_v) +end subroutine MOM_read_vector_3d_fms2 + +!> apply a scale factor to a 1d array +subroutine scale_data_1d(data, scale_factor) + real, dimension(:), intent(inout) :: data !< The 1-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + + if (scale_factor /= 1.0) then + data(:) = scale_factor*data(:) + endif +end subroutine scale_data_1d + +!> apply a scale factor to a 2d array +subroutine scale_data_2d(data, scale_factor, MOM_domain) + real, dimension(:,:), intent(inout) :: data !< The 2-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je) = scale_factor*data(is:ie,js:je) + else + data(:,:) = scale_factor*data(:,:) + endif + endif +end subroutine scale_data_2d + +!> apply a scale factor to a 3d array +subroutine scale_data_3d(data, scale_factor, MOM_domain) + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:) = scale_factor*data(is:ie,js:je,:) + else + data(:,:,:) = scale_factor*data(:,:,:) + endif + endif +end subroutine scale_data_3d + +!> apply a scale factor to a 4d array +subroutine scale_data_4d(data, scale_factor, MOM_domain) + real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:,:) = scale_factor*data(is:ie,js:je,:,:) + else + data(:,:,:,:) = scale_factor*data(:,:,:,:) + endif + endif +end subroutine scale_data_4d + +!> check that latitude or longitude units are valid CF-compliant values +!! return true or false and x_or_y character value corresponding to the axis direction +subroutine validate_lat_lon_units(unit_string, x_or_y, units_are_valid) +character(len=*), intent(in) :: unit_string !< string of units +character(len=1), intent(out) :: x_or_y !< "x" for longitude or "y" latitude +logical, intent(out) :: units_are_valid !< .true. if units match acceptable values; default is .false. + +select case (lowercase(trim(unit_string))) + case ("degrees_north"); units_are_valid = .true.; x_or_y = "y" + case ("degree_north"); units_are_valid = .true.; x_or_y = "y" + case ("degrees_n"); units_are_valid = .true.; x_or_y = "y" + case ("degree_n"); units_are_valid = .true.; x_or_y = "y" + case ("degreen"); units_are_valid = .true.; x_or_y = "y" + case ("degreesn"); units_are_valid = .true.; x_or_y = "y" + case ("degrees_east"); units_are_valid = .true.; x_or_y = "x" + case ("degree_east"); units_are_valid = .true.;x_or_y = "x" + case ("degreese"); units_are_valid = .true.; x_or_y = "x" + case ("degreee"); units_are_valid = .true.; x_or_y = "x" + case ("degree_e"); units_are_valid = .true.; x_or_y = "x" + case ("degrees_e"); units_are_valid = .true.; x_or_y = "x" + case default; units_are_valid = .false.; x_or_y = "" +end select + +end subroutine validate_lat_lon_units + +end module MOM_read_data_fms2 diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index ed29b99b55..07c054351a 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -2,34 +2,47 @@ module MOM_restart ! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_domains, only : pe_here, num_PEs 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_string_functions, only : lowercase +use MOM_string_functions, only : lowercase, append_substring use MOM_grid, only : ocean_grid_type use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file -use MOM_io, only : MOM_read_data, read_data, get_filename_appendix +use MOM_io, only : MOM_read_data, read_data, get_filename_appendix ! NOTE get_filename_appendix is not in fms2-io use MOM_io, only : get_file_info, get_file_atts, get_file_fields, get_file_times use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE +use MOM_axis, only : get_time_units, convert_checksum_to_string +use MOM_axis, only : axis_data_type, MOM_get_diagnostic_axis_data +use MOM_axis, only : MOM_register_diagnostic_axis, get_var_dimension_metadata use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : days_in_month, get_date, set_date use MOM_transform_FMS, only : mpp_chksum => rotated_mpp_chksum use MOM_transform_FMS, only : write_field => rotated_write_field use MOM_verticalGrid, only : verticalGrid_type +use mpp_domains_mod, only: mpp_define_io_domain, mpp_get_domain_npes, mpp_get_io_domain +use mpp_domains_mod, only: mpp_get_compute_domain, mpp_get_global_domain use mpp_io_mod, only : mpp_attribute_exist, mpp_get_atts -use mpp_mod, only : mpp_pe - +use mpp_mod, only: mpp_pe, mpp_max +! fms2-io interfaces +use fms2_io_mod, only : fms2_register_restart_field => register_restart_field +use fms2_io_mod, only : check_if_open, is_dimension_registered, register_field, register_axis +use fms2_io_mod, only : register_variable_attribute, read_data, read_restart, write_restart +use fms2_io_mod, only : write_data, fms2_close_file=>close_file, fms2_open_file=>open_file +use fms2_io_mod, only : global_att_exists, get_global_attribute, get_global_io_domain_indices +use fms2_io_mod, only : get_dimension_names, get_dimension_size, get_num_dimensions, variable_exists +use fms2_io_mod, only : dimension_exists, FmsNetcdfDomainFile_t, unlimited, get_variable_size +use fms2_io_mod, only : get_variable_num_dimensions, get_variable_dimension_names + +use platform_mod implicit none ; private public restart_init, restart_end, restore_state, register_restart_field public save_restart, query_initialized, restart_init_end, vardesc public restart_files_exist, determine_is_new_run, is_new_run public register_restart_field_as_obsolete +public write_initial_conditions public register_restart_pair - !> A type for making arrays of pointers to 4-d arrays type p4d real, dimension(:,:,:,:), pointer :: p => NULL() !< A pointer to a 4d array @@ -848,8 +861,32 @@ function query_initialized_4d_name(f_ptr, name, CS) result(query_initialized) end function query_initialized_4d_name +!> wrapper routine for save_restart_old, save_restart_fms2, and write_initial_conditions_file +subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, use_fms2, write_ic) + character(len=*), intent(in) :: directory !< The directory where the restart files + !! are to be written + type(time_type), intent(in) :: time !< The current model time + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + logical, optional, intent(in) :: time_stamped !< If present and true, add time-stamp + !! to the restart file names. + character(len=*), optional, intent(in) :: filename !< A filename that overrides the name in CS%restartfile. + type(verticalGrid_type), optional, intent(in) :: GV !< The ocean's vertical grid structure + logical, optional, intent(in) :: use_fms2 !< flag to call save_restart_fms2 + logical, optional, intent(in) :: write_ic !< flag to call write_initial_conditions + + if (present(write_ic) .and. write_ic) then + call write_initial_conditions(directory, time, G, CS, time_stamped=time_stamped, filename=filename, GV=GV) + elseif (present(use_fms2) .and. use_fms2) then + call save_restart_fms2(directory, time, G, CS, time_stamped=time_stamped, filename=filename, GV=GV) + else + call save_restart_old(directory, time, G, CS, time_stamped=time_stamped, filename=filename, GV=GV) + endif +end subroutine save_restart + !> save_restart saves all registered variables to restart files. -subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) +subroutine save_restart_old(directory, time, G, CS, time_stamped, filename, GV) character(len=*), intent(in) :: directory !< The directory where the restart files !! are to be written type(time_type), intent(in) :: time !< The current model time @@ -1056,12 +1093,488 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) num_files = num_files+1 enddo -end subroutine save_restart +end subroutine save_restart_old + +!> save all registered variables to a restart file using fms2-io +subroutine save_restart_fms2(directory, time, G, CS, time_stamped, filename, GV) + character(len=*), intent(in) :: directory !< The directory where the restart files + !! are to be written + type(time_type), intent(in) :: time !< The current model time + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + logical, optional, intent(in) :: time_stamped !< If present and true, add time-stamp + !! to the restart file names. + character(len=*), optional, intent(in) :: filename !< A filename that overrides the name in CS%restartfile. + type(verticalGrid_type), optional, intent(in) :: GV !< The ocean's vertical grid structure + + ! Local variables + type(vardesc) :: vars(CS%max_fields) ! Descriptions of the fields that + ! are to be read from the restart file. + type(fieldtype) :: fields(CS%max_fields) ! + type(FmsNetcdfDomainFile_t) :: fileObjWrite ! netcdf file object returned by a call to open_file + character(len=1024) :: restartpath ! The restart file path (dir/file). + character(len=512) :: restartname ! The restart file name (no dir). + character(len=700) :: restartpath_temp ! temporary location for the restart file path (dir/file). + character(len=600) :: restartname_temp ! temporary location for restart name + character(len=512) :: base_file_name ! Temporary location for restart file name (no dir) + character(len=8) :: suffix ! A suffix (like _2) that is appended + ! to the name of files after the first. + integer(kind=8) :: var_sz, size_in_file ! The size in bytes of each variable + ! and the variables already in a file. + integer(kind=8) :: max_file_size = 2147483647_8 ! The maximum size in bytes + ! for any one file. With NetCDF3, + ! this should be 2 Gb or less. + integer :: start_var, next_var ! The starting variables of the + ! current and next files. + integer :: unit ! The mpp unit of the open file. + integer :: m, nz, i, k, num_files + integer :: seconds, days, year, month, hour, minute + character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. + character(len=8) :: t_grid_read + character(len=64) :: var_name ! A variable's name. + character(len=256) :: date_appendix ! date string to append to a file name if desired + character(len=64) :: dim_names(4) ! Array to hold up to 4 strings for the variable axis names + integer, dimension(4) :: dim_lengths ! Array of integer lengths corresponding to the name(s) in axis_names + integer :: name_length + integer(kind=8) :: check_val(CS%max_fields,1) + integer :: is, ie + integer :: substring_index + integer :: horgrid_position + integer :: num_dims, total_axes + integer :: var_periods + logical :: fileOpenSuccess ! true if netcdf file is opened + real :: restart_time + character(len=32) :: filename_appendix = '' ! fms appendix to filename for ensemble runs + character(len=16) :: restart_time_units + character(len=64) :: checksum_char + character(len=64) :: units + character(len=256) :: longname + real, dimension(:), allocatable :: data_temp + type(axis_data_type) :: axis_data_CS + integer :: isL, ieL, jsL, jeL, pos + integer :: turns + + turns = CS%turns + + if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & + "save_restart_fms2: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) + + ! With parallel read & write, it is possible to disable the following... + + ! The maximum file size is 4294967292, according to the NetCDF documentation. + if (CS%large_file_support) max_file_size = 4294967292_8 + + horgrid_position = 1 + name_length = 0 + num_files = 0 + restartname = "" + base_file_name = "" + restartname_temp = "" + date_appendix = "" + restart_time_units = "" + + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(G%domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(G%domain%mpp_domain))) & + call mpp_define_io_domain(G%domain%mpp_domain, (/1,1/)) + endif + ! get the number of vertical levels + nz = 1 ; if (present(GV)) nz = GV%ke + + if (present(filename)) then + base_file_name = trim(filename) + else + base_file_name=trim(CS%restartfile) + endif + ! append a time stamp to the file name if time_stamp is specified + if (PRESENT(time_stamped)) then + if (time_stamped) then + call get_date(time,year,month,days,hour,minute,seconds) + ! Compute the year-day, because I don't like months. - RWH + do m=1,month-1 + days = days + days_in_month(set_date(year,m,2,0,0,0)) + enddo + seconds = seconds + 60*minute + 3600*hour + if (year <= 9999) then + write(date_appendix,'("_Y",I4.4,"_D",I3.3,"_S",I5.5)') year, days, seconds + elseif (year <= 99999) then + write(date_appendix,'("_Y",I5.5,"_D",I3.3,"_S",I5.5)') year, days, seconds + else + write(date_appendix,'("_Y",I10.10,"_D",I3.3,"_S",I5.5)') year, days, seconds + endif + restartname_temp = trim(base_file_name)//trim(date_appendix) + endif + else + restartname_temp = trim(base_file_name) + endif + + ! get the restart time units + restart_time = time_type_to_real(time) / 86400.0 + restart_time_units = "days" + next_var = 1 + do while (next_var <= CS%novars ) + start_var = next_var + ! get variable sizes in bytes + size_in_file = 8*(2*G%Domain%niglobal+2*G%Domain%njglobal+2*nz+1000) + + do m=start_var,CS%novars + call query_vardesc(CS%restart_field(m)%vars, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid, caller="save_restart") + if (hor_grid == '1') then + var_sz = 8 + else + var_sz = 8*(G%Domain%niglobal+1)*(G%Domain%njglobal+1) + endif + select case (z_grid) + case ('L') ; var_sz = var_sz * nz + case ('i') ; var_sz = var_sz * (nz+1) + end select + t_grid = adjustl(t_grid) + if (t_grid(1:1) == 'p') then + if (len_trim(t_grid(2:8)) > 0) then + var_periods = -1 + t_grid_read = adjustl(t_grid(2:8)) + read(t_grid_read,*) var_periods + if (var_periods > 1) var_sz = var_sz * var_periods + endif + endif + + if ((m==start_var) .OR. (size_in_file < max_file_size-var_sz)) then + size_in_file = size_in_file + var_sz + else ; exit + endif + + enddo + next_var = m + + restartpath = "" + restartpath_temp = "" + suffix = "" + + !query fms_io if there is a filename_appendix (for ensemble runs) + ! TODO move filename_appendix functionality to fms2-io or MOM6 framework + name_length = len_trim(restartname_temp) + call get_filename_appendix(filename_appendix) + if (len_trim(filename_appendix) > 0) then + if (restartname_temp(name_length-2:name_length) == '.nc') then + restartname = restartname_temp(1:name_length-3)//'.'//trim(filename_appendix)//'.nc' + else + if (trim(filename_appendix) .ne. " ") then + restartname = restartname_temp(1:name_length) //'.'//trim(filename_appendix) + else + restartname(1:name_length) = trim(restartname_temp) + endif + endif + else + restartname(1:name_length) = trim(restartname_temp) + endif + + if (num_files < 10) then + write(suffix,'("_",I1)') num_files + else + write(suffix,'("_",I2)') num_files + endif + + if (num_files .gt. 0) then + name_length = len_trim(directory//restartname//suffix) + restartpath_temp = trim(directory)//trim(restartname)//trim(suffix) + else + name_length = len_trim(directory//restartname) + restartpath_temp = trim(directory)//trim(restartname) + endif + ! append '.nc' to the restart file path if it is missing + substring_index = index(trim(restartpath_temp), ".nc") + if (substring_index <= 0) then + restartpath = append_substring(restartpath_temp,".nc") + else + restartpath(1:len_trim(restartpath_temp)) = trim(restartpath_temp) + endif + ! create the file and register and write the global axes to the file + if (present(GV)) then + call create_file(trim(restartpath), fileObjWrite, CS%restart_field%vars, CS%novars, register_time=.true., & + G=G, GV=GV, is_restart=.true.) + else + call create_file(trim(restartpath), fileObjWrite, CS%restart_field%vars, CS%novars, register_time=.true., & + G=G, is_restart=.true.) + endif + ! register the time data + if (.not. variable_exists(fileObjWrite, "Time")) then + call register_field(fileObjWrite, "Time", "double", dimensions=(/"Time"/)) + call register_variable_attribute(fileObjWrite, "Time", "units", restart_time_units) + endif + + do m=start_var,next_var-1 + vars(m-start_var+1) = CS%restart_field(m)%vars + enddo + + call query_vardesc(vars(1), t_grid=t_grid, hor_grid=hor_grid, caller="save_restart") + + t_grid = adjustl(t_grid) + if (t_grid(1:1) /= 'p') & + call modify_vardesc(vars(1), t_grid='s', caller="save_restart") + select case (hor_grid) + case ('q') ; pos = CORNER + case ('h') ; pos = CENTER + case ('u') ; pos = EAST_FACE + case ('v') ; pos = NORTH_FACE + case ('Bu') ; pos = CORNER + case ('T') ; pos = CENTER + case ('Cu') ; pos = EAST_FACE + case ('Cv') ; pos = NORTH_FACE + case ('1') ; pos = 0 + case default ; pos = 0 + end select + !Prepare the checksum of the restart fields to be written to restart files + if (modulo(turns, 2) /= 0) then + call get_checksum_loop_ranges(G, pos, jsL, jeL, isL, ieL) + else + call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) + endif + + do m=start_var,next_var-1 + if (associated(CS%var_ptr3d(m)%p)) then + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) + elseif (associated(CS%var_ptr2d(m)%p)) then + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) + elseif (associated(CS%var_ptr4d(m)%p)) then + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) + elseif (associated(CS%var_ptr1d(m)%p)) then + check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr1d(m)%p) + elseif (associated(CS%var_ptr0d(m)%p)) then + check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr0d(m)%p,pelist=(/mpp_pe()/)) + endif + enddo + + do m=start_var,next_var-1 + longname = "" + num_dims = 0 + units = "" + dim_names(:) = "" + if (.not.(variable_exists(fileObjWrite, CS%restart_field(m)%var_name))) then + call query_vardesc(vars(m-start_var+1), hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid, longname=longname, & + units=units, caller="save_restart") + + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, & + dim_names, dim_lengths, num_dims, G=G, GV=GV) + ! register the restart variables to the file + if (associated(CS%var_ptr3d(m)%p)) then + call fms2_register_restart_field(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + CS%var_ptr3d(m)%p, dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr2d(m)%p)) then + call fms2_register_restart_field(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + CS%var_ptr2d(m)%p, dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr4d(m)%p)) then + call fms2_register_restart_field(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + CS%var_ptr4d(m)%p, dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr1d(m)%p)) then + ! need to pass dim_names argument as a 1-D array + call fms2_register_restart_field(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + CS%var_ptr1d(m)%p, dimensions=(/dim_names(1:num_dims)/)) + elseif (associated(CS%var_ptr0d(m)%p)) then + ! need to pass dim_names argument as a 1-D array + call fms2_register_restart_field(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + CS%var_ptr0d(m)%p, dimensions=(/dim_names(1:num_dims)/)) + endif + ! convert the checksum to a string + checksum_char = '' + checksum_char = convert_checksum_to_string(check_val(m,1)) + !! register the variable attributes + !call register_variable_attribute(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + ! 'checksum', trim(checksum_char)) + call register_variable_attribute(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + 'units', units) + call register_variable_attribute(fileObjWrite, CS%restart_field(m-start_var+1)%var_name, & + 'long_name', longname) + endif + enddo + ! write the time data + call write_data(fileObjWrite, "Time", (/restart_time/)) + ! write the restart file + call write_restart(fileObjWrite) + ! close the file + if (check_if_open(fileObjWrite)) call fms2_close_file(fileObjWrite) + + if (associated(axis_data_CS%axis)) deallocate(axis_data_CS%axis) + if (associated(axis_data_CS%data)) deallocate(axis_data_CS%data) + + num_files = num_files+1 + enddo + +end subroutine save_restart_fms2 + +!> write initial condition fields to a netCDF file +subroutine write_initial_conditions(directory, time, G, CS, time_stamped, filename, GV) + character(len=*), intent(in) :: directory !< The directory where the restart files + !! are to be written + type(time_type), intent(in) :: time !< The current model time + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + logical, optional, intent(in) :: time_stamped !< If present and true, add time-stamp + !! to the restart file names. + character(len=*), optional, intent(in) :: filename !< A filename that overrides the name in CS%restartfile. + type(verticalGrid_type), optional, intent(in) :: GV !< The ocean's vertical grid structure + ! local + type(vardesc) :: vd ! structure for variable metadata + type(FmsNetcdfDomainFile_t) :: fileObjWrite ! netCDF file object returned by call to open_file + type(axis_data_type) :: axis_data_CS ! structure for coordinate variable metadata + integer :: substring_index + integer :: name_length + integer :: num_dims ! counter for variable dimensions + integer :: total_axes ! counter for all coordinate axes in file + integer :: i, is, ie, k, m, isc, jsc, iec, jec, isg, jsg, ieg, jeg + integer :: var_periods + integer, dimension(4) :: dim_lengths + integer, allocatable :: pos(:),first(:,:), last(:,:) + logical :: fileOpenSuccess ! .true. if netcdf file is opened + character(len=200) :: base_file_name + character(len=200) :: dim_names(4) + character(len=20) :: time_units + character(len=64) :: units + character(len=256) :: longname + character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. + real :: ic_time + real, dimension(:), allocatable :: data_temp + + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(G%domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(G%domain%mpp_domain))) & + call mpp_define_io_domain(G%domain%mpp_domain, (/1,1/)) + endif + ! append '.nc' to the restart file name if it is missing + ! TODO: require users to specify full file path including the file name appendix + ! in calls to open_file + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + base_file_name = append_substring(trim(directory)//trim(filename),".nc") + else + name_length = len(trim(directory)//trim(filename)) + base_file_name(1:name_length) = trim(directory)//trim(filename) + endif + ! get the time units + ic_time = time_type_to_real(time) / 86400.0 + time_units = get_time_units(ic_time*86400.0) + ! create the file and register and write the global axes to the file + if (present(GV)) then + call create_file(trim(base_file_name), fileObjWrite, CS%restart_field%vars, CS%novars, register_time=.true., & + G=G, GV=GV) + else + call create_file(trim(base_file_name), fileObjWrite, CS%restart_field%vars, CS%novars, register_time=.true., G=G) + endif + ! register the time data + if (.not. variable_exists(fileObjWrite, "Time")) then + call register_field(fileObjWrite, "Time", "double", dimensions=(/"Time"/)) + call register_variable_attribute(fileObjWrite, "Time", "units", time_units) + endif + ! allocate position indices for x- and y-dimensions associated with variables + allocate(pos(CS%novars)) + allocate(first(CS%novars,2)); allocate(last(CS%novars,2)); + first(:,:) = 0; last(:,:) = 0 + pos(:) = CENTER + ! register and write the field variables to the initial conditions file + do m=1,CS%novars + longname = "" + num_dims = 0 + units = "" + dim_names(:) = "" + + call query_vardesc(CS%restart_field(m)%vars, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid, longname=longname, & + units=units, caller="save_restart") + + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, & + dim_names, dim_lengths, num_dims, G=G, GV=GV) + select case (hor_grid) + case ('q') ; pos(m) = CORNER + case ('h') ; pos(m) = CENTER + case ('u') ; pos(m) = EAST_FACE + case ('v') ; pos(m) = NORTH_FACE + case ('Bu') ; pos(m) = CORNER + case ('T') ; pos(m) = CENTER + case ('Cu') ; pos(m) = EAST_FACE + case ('Cv') ; pos(m) = NORTH_FACE + case ('1') ; pos(m) = 0 + case default ; pos(m)= 0 + end select + ! register the variables + if (associated(CS%var_ptr3d(m)%p)) then + call register_field(fileObjWrite, CS%restart_field(m)%var_name, "double", & + dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr2d(m)%p)) then + call register_field(fileObjWrite, CS%restart_field(m)%var_name, "double", & + dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr4d(m)%p)) then + call register_field(fileObjWrite, CS%restart_field(m)%var_name, "double", & + dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr1d(m)%p)) then + ! need to explicitly define dim_names array for 1-D variable + call register_field(fileObjWrite, CS%restart_field(m)%var_name, "double", & + dimensions=(/dim_names(1)/)) + elseif (associated(CS%var_ptr0d(m)%p)) then + ! need to explicitly define dim_names array for scalar variable + call register_field(fileObjWrite, CS%restart_field(m)%var_name, "double", & + dimensions=(/dim_names(1)/)) + endif + ! register the variable attributes + call register_variable_attribute(fileObjWrite, CS%restart_field(m)%var_name, "units", units) + call register_variable_attribute(fileObjWrite, CS%restart_field(m)%var_name, "long_name", longname) + enddo + + do m=1,CS%novars + if (associated(CS%var_ptr3d(m)%p)) then + call write_data(fileObjWrite, CS%restart_field(m)%var_name, CS%var_ptr3d(m)%p, & + unlim_dim_level=1) + elseif (associated(CS%var_ptr2d(m)%p)) then + call write_data(fileObjWrite, CS%restart_field(m)%var_name, CS%var_ptr2d(m)%p, & + unlim_dim_level=1) + elseif (associated(CS%var_ptr4d(m)%p)) then + call write_data(fileObjWrite, CS%restart_field(m)%var_name, CS%var_ptr4d(m)%p, & + unlim_dim_level=1) + elseif (associated(CS%var_ptr1d(m)%p)) then + call write_data(fileObjWrite, CS%restart_field(m)%var_name, CS%var_ptr1d(m)%p, & + unlim_dim_level=1) + elseif (associated(CS%var_ptr0d(m)%p)) then + call write_data(fileObjWrite, CS%restart_field(m)%var_name,CS%var_ptr0d(m)%p) + endif + enddo + ! write the time data + call write_data(fileObjWrite, "Time", (/ic_time/)) + ! close the IC file and deallocate the allocatable arrays + if(check_if_open(fileObjWrite)) call fms2_close_file(fileObjWrite) + + if (associated(axis_data_CS%axis)) deallocate(axis_data_CS%axis) + if (associated(axis_data_CS%data)) deallocate(axis_data_CS%data) + deallocate(pos); deallocate(first); deallocate(last) +end subroutine write_initial_conditions + +!> wrapper routine for restore_state_old and restore_state_fms2 +subroutine restore_state(filename, directory, day, G, CS, use_fms2) + character(len=*), intent(in) :: filename !< The list of restart file names or a single + !! character 'r' to read automatically named files. + character(len=*), intent(in) :: directory !< The directory in which to find restart files + type(time_type), intent(out) :: day !< The time of the restarted run + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + logical, optional, intent(in) :: use_fms2 !< if .true., call restore_state_fms2 + + if (present(use_fms2) .and. use_fms2) then + call restore_state_fms2(filename, directory, day, G, CS) + else + call restore_state_old(filename, directory, day, G, CS) + endif +end subroutine restore_state !> restore_state reads the model state from previously generated files. All !! restart variables are read from the first file in the input filename list !! in which they are found. -subroutine restore_state(filename, directory, day, G, CS) +subroutine restore_state_old(filename, directory, day, G, CS) character(len=*), intent(in) :: filename !< The list of restart file names or a single !! character 'r' to read automatically named files. character(len=*), intent(in) :: directory !< The directory in which to find restart files @@ -1282,9 +1795,202 @@ subroutine restore_state(filename, directory, day, G, CS) endif enddo -end subroutine restore_state +end subroutine restore_state_old + +!> restore_state_fms2 reads the model state from previously generated files using fms2-io. All +!! restart variables are read from the first file in the input filename list +!! in which they are found. +subroutine restore_state_fms2(filename, directory, day, G, CS) + character(len=*), intent(in) :: filename !< The list of restart file names or a single + !! character 'r' to read automatically named files. + character(len=*), intent(in) :: directory !< The directory in which to find restart files + type(time_type), intent(out) :: day !< The time of the restarted run + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + + ! This subroutine reads the model state from previously + ! generated files. All restart variables are read from the first + ! file in the input filename list in which they are found. + ! Local variables + character(len=200) :: filepath ! The path (dir/file) to the file being opened. + character(len=80) :: fname ! The name of the current file. + character(len=8) :: suffix ! A suffix (like "_2") that is added to any + ! additional restart files. + character(len=512) :: mesg ! A message for warnings. + character(len=80) :: varname ! A variable's name. + integer :: i, m, n + integer :: isL, ieL, jsL, jeL, is0, js0 + integer :: ntime, pos + character(len=200) :: unit_path(CS%max_fields) ! The file names. + logical :: unit_is_global(CS%max_fields) ! True if the file is global. + character(len=200) :: base_file_name + character(len=1024) :: temp_file_name + character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. + real :: t1, t2 ! Two times. + real, allocatable :: time_vals(:) + logical :: check_exist, is_there_a_checksum + integer(l8_kind),dimension(3) :: checksum_file + integer(kind=8) :: checksum_data + integer :: missing_fields + logical :: fileOpenSuccess ! .true. if netcdf file object is opened + type(FmsNetcdfDomainFile_t) :: fileObjRead ! netcdf file object returned by open_file + integer :: str_index, num_file, is,ie,js,je + character(len=64) :: checksum_char, time_units + character(len=20), dimension(:), allocatable :: axis_names + character(len=32) :: dim_names(4) + integer :: dim_lengths(4), num_dims, dim_size + + if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & + "restore_state: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) + ! define the io domain if using 1 pe and the io domain is not set + if (mpp_get_domain_npes(G%domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(G%domain%mpp_domain))) & + call mpp_define_io_domain(G%domain%mpp_domain, (/1,1/)) + endif + + str_index = 0 + ! get the base restart file name + temp_file_name='' + if ((LEN_TRIM(filename) == 1 .and. filename(1:1) == 'F') .or. (trim(filename)=='r')) then + temp_file_name = trim(CS%restartfile) + else + temp_file_name = trim(filename) + endif + ! append '.nc.' to the file name if it is missing + base_file_name = "" + str_index = INDEX(temp_file_name, ".nc") + if (str_index <=0) then + base_file_name = trim(append_substring(temp_file_name, ".nc")) + else + base_file_name = trim(temp_file_name) + endif + + num_file = get_num_restart_files(temp_file_name, directory, G, CS, file_paths=unit_path) + CS%restart_field(:)%initialized = .false. + ! Read each variable from the first file in which it is found. + do n=1,num_file + ! Open the restart file. + if (.not.(check_if_open(fileObjRead))) & + fileOpenSuccess=fms2_open_file(fileObjRead, trim(unit_path(n)), "read", & + G%domain%mpp_domain, is_restart=.true.) + if (fileOpenSuccess) & + call MOM_error(NOTE, "MOM_restart_fms2: MOM run restarted using : "//trim(unit_path(n))) + + call get_dimension_size(fileObjRead, "Time", ntime) + + if (ntime .lt. 1) then + call MOM_error(NOTE, "MOM_restart_fms2: time is scalar.") + ntime=1 + endif + allocate(time_vals(ntime)) + call read_data(fileObjRead, "Time", time_vals) + t1 = time_vals(1) + deallocate(time_vals) + t2 = t1 + call mpp_max(t2) + if (t1 .ne. t2) then + call MOM_error(FATAL, "times are different in different restart files.") + endif + + day = real_to_time(t1*86400.0) + ! Register the horizontal axes that correspond to x and y of the domain. + num_dims=get_num_dimensions(fileObjRead) + allocate(axis_names(num_dims)) + axis_names(:)= "" + call get_dimension_names(fileObjRead, axis_names) + do i = 1,num_dims + call get_dimension_size(fileObjRead, trim(axis_names(i)), dim_size) + call MOM_register_diagnostic_axis(fileObjRead, trim(axis_names(i)), dim_size) + enddo + ! Read in each variable from the restart files. + missing_fields = 0 + do m = 1, CS%novars + varname = '' + varname = trim(CS%restart_field(m)%var_name) + ! Check for obsolete fields + do i = 1,CS%num_obsolete_vars + if (adjustl(lowercase(trim(varname))) .eq. adjustl(lowercase(trim(CS%restart_obsolete(i)%field_name)))) then + call MOM_error(FATAL, "MOM_restart:restore_state_fms2: Attempting to use obsolete restart field "//& + trim(varname)//" - the new corresponding restart field is "//& + trim(CS%restart_obsolete(i)%replacement_name)) + endif + enddo + + if (CS%restart_field(m)%initialized) cycle + call query_vardesc(CS%restart_field(m)%vars, hor_grid=hor_grid, & + caller="restore_state_fms2") + select case (hor_grid) + case ('q') ; pos = CORNER + case ('h') ; pos = CENTER + case ('u') ; pos = EAST_FACE + case ('v') ; pos = NORTH_FACE + case ('Bu') ; pos = CORNER + case ('T') ; pos = CENTER + case ('Cu') ; pos = EAST_FACE + case ('Cv') ; pos = NORTH_FACE + case ('1') ; pos = 0 + case default ; pos = 0 + end select + + call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) + ! Check if the variable is mandatory and present in the restart file(s) + if (.not. variable_exists(fileObjRead, trim(varname))) then + if (CS%restart_field(m)%mand_var) then + call MOM_error(WARNING, "MOM_restart_fms2: Unable to find mandatory variable " & + //trim(varname)//" in restart file "//trim(directory)//trim(base_file_name)) + missing_fields = missing_fields+1 + cycle + endif + endif + ! Get the variable's "domain position." + num_dims = 0 + dim_names(:) = "" + num_dims=get_variable_num_dimensions(fileobjRead, trim(varname)) + call get_variable_dimension_names(fileObjRead, trim(varname), dim_names(1:num_dims)) + ! Register the restart fields and compute the checksums. + if (associated(CS%var_ptr1d(m)%p)) then + call fms2_register_restart_field(fileObjRead, trim(varname), CS%var_ptr1d(m)%p, & + dimensions=(/dim_names(1)/)) + elseif (associated(CS%var_ptr0d(m)%p)) then + call fms2_register_restart_field(fileObjRead, trim(varname), CS%var_ptr0d(m)%p) + elseif (associated(CS%var_ptr2d(m)%p)) then + call fms2_register_restart_field(fileObjRead, trim(varname), CS%var_ptr2d(m)%p, & + dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr3d(m)%p)) then + call fms2_register_restart_field(fileObjRead, trim(varname), CS%var_ptr3d(m)%p, & + dimensions=dim_names(1:num_dims)) + elseif (associated(CS%var_ptr4d(m)%p)) then + call fms2_register_restart_field(fileObjRead, trim(varname), CS%var_ptr4d(m)%p, & + dimensions=dim_names(1:num_dims)) + else + call MOM_error(FATAL, "MOM_restart restore_state_fms2: No pointers set for "//trim(varname)) + endif + CS%restart_field(m)%initialized = .true. + enddo ! m=CS%novars + ! Read in restart data and then close the file. + call read_restart(fileObjRead, unlim_dim_level=1) + ! close the file + if (check_if_open(fileObjRead)) call fms2_close_file(fileObjRead) + if (allocated(axis_names)) deallocate(axis_names) + if (missing_fields == 0) exit + enddo + + do m=1,CS%novars + if (.not.(CS%restart_field(m)%initialized)) then + CS%restart = .false. + if (CS%restart_field(m)%mand_var) then + call MOM_error(FATAL,"MOM_restart: Unable to find mandatory variable " & + //trim(CS%restart_field(m)%var_name)//" in restart files.") + endif + endif + enddo + +end subroutine restore_state_fms2 !> restart_files_exist determines whether any restart files exist. +! TODO remove this function when fms2-io is fully implemented function restart_files_exist(filename, directory, G, CS) character(len=*), intent(in) :: filename !< The list of restart file names or a single !! character 'r' to read automatically named files. @@ -1497,6 +2203,136 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & end function open_restart_units +!> get_num_restart_files determines the number of existing restart files and returns paths +!! and whether the files are global or spatially decomposed. +function get_num_restart_files(filename, directory, G, CS, file_paths) result(num_files) + character(len=*), intent(in) :: filename !< The list of restart file names or a single + !! character 'r' to read automatically named files. + character(len=*), intent(in) :: directory !< The directory in which to find restart files + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous + !! call to restart_init. + character(len=*), dimension(:), & + optional, intent(out) :: file_paths !< The full paths to open files. + !logical, dimension(:), & + ! optional, intent(out) :: global_files !< True if a file is global. + + integer :: num_files !< The number of files (both automatically named restart + !! files and others explicitly in filename) that have been opened. + +! This subroutine reads the model state from previously +! generated files. All restart variables are read from the first +! file in the input filename list in which they are found. + + ! Local variables + character(len=256) :: filepath ! The path (dir/file) to the file being opened. + character(len=256) :: fname ! The name of the current file. + character(len=8) :: suffix ! A suffix (like "_2") that is added to any + ! additional restart files + integer :: num_restart ! The number of restart files that have already + ! been opened. + integer :: start_char ! The location of the starting character in the + ! current file name. + integer :: f, n, m, err, length, str_index + logical :: fexists + character(len=32) :: filename_appendix = '' !fms appendix to filename for ensemble runs + character(len=80) :: restartname + character(len=240) :: filepath_temp, filepath_temp2 + + if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & + "get_num_restart_files: Module must be initialized before it is used.") + ! Determine the file name + num_restart = 0 ; n=0; start_char = 1; str_index=0 + if (present(file_paths)) file_paths(:) = "" + do while (start_char <= len_trim(filename) ) + do m=start_char,len_trim(filename) + if (filename(m:m) == ' ') exit + enddo + fname = filename(start_char:m-1) + start_char = m + do while (start_char <= len_trim(filename)) + if (filename(start_char:start_char) == ' ') then + start_char = start_char + 1 + else + exit + endif + enddo + + err = 0 + if (num_restart > 0) err = 1 ! Avoid going through the file list twice. + do while (err == 0) + restartname = trim(CS%restartfile) + ! query fms_io if there is a filename_appendix (for ensemble runs) + ! TODO add support to fms2-io, or move to MOM6 framework + call get_filename_appendix(filename_appendix) + if (len_trim(filename_appendix) > 0 .and. trim(filename_appendix) .ne. " ") then + length = len_trim(restartname) + if (restartname(length-2:length) == '.nc') then + restartname = restartname(1:length-3)//'.'//trim(filename_appendix)//'.nc' + else + restartname = restartname(1:length) //'.'//trim(filename_appendix) + endif + endif + filepath = trim(directory) // trim(restartname) + + if (num_restart < 10) then + write(suffix,'("_",I1)') num_restart + else + write(suffix,'("_",I2)') num_restart + endif + if (num_restart > 0) filepath = trim(filepath) // suffix + + filepath_temp = trim(filepath)//".nc" + if (file_exists(trim(filepath_temp),.true.) .or. file_exists(trim(filepath_temp)//".0000",.true.)) then + n = n+1 + if (present(file_paths)) file_paths(n) = trim(filepath_temp) + num_restart = num_restart + 1 + call MOM_error(NOTE, "MOM_restart:get_num_restart_files: Found restart file : "//trim(filepath)) + endif + ! search for files with "res_#" in the name + str_index = index(filepath_temp,".res.nc") + if (str_index .gt. 0) then + f = 0 + do while (f .le. n) + f=f+1 + filepath_temp2="" + ! check for names with extra .res.nc added by fms2-io + if ( f .lt. 10) then + write(filepath_temp2,'(A,I1,A)') trim(filepath_temp(1:str_index-1))//".res_",f,".res.nc" + elseif (f .ge. 10 .and. f .lt. 100) then + write(filepath_temp2,'(A,I2,A)') trim(filepath_temp(1:str_index-1))//".res_",f,".res.nc" + endif + if (file_exists(trim(filepath_temp2),.true.) .or. file_exists(trim(filepath_temp2)//".0000",.true.)) then + call MOM_error(NOTE, "MOM_restart:get_num_restart_files: Found restart file : "//trim(filepath_temp2)) + num_restart=num_restart+1 + n=n+1 + if (present(file_paths)) file_paths(n) = trim(filepath_temp2) + else + ! check for fms-io-style name + filepath_temp2="" + if ( f .lt. 10) then + write(filepath_temp2,'(A,I1,A)') trim(filepath_temp(1:str_index-1))//".res_",f,".nc" + elseif (f .ge. 10 .and. f .lt. 100) then + write(filepath_temp2,'(A,I2,A)') trim(filepath_temp(1:str_index-1))//".res_",f,".nc" + endif + if (file_exists(trim(filepath_temp2),.true.) .or. file_exists(trim(filepath_temp2)//".0000",.true.)) then + call MOM_error(NOTE, "MOM_restart:get_num_restart_files: Found restart file : "//trim(filepath_temp2)) + num_restart=num_restart+1 + n=n+1 + if (present(file_paths)) file_paths(n) = trim(filepath_temp2) + else + exit + endif + endif + enddo ! while (f .le. n-1) + endif + err = 1 ; exit + enddo ! while (err == 0) loop + enddo ! while (start_char < strlen(filename)) loop + num_files = n + +end function get_num_restart_files + !> Initialize this module and set up a restart control structure. subroutine restart_init(param_file, CS, restart_root) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters @@ -1650,4 +2486,41 @@ subroutine get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) end subroutine get_checksum_loop_ranges +!> get the size of a variable in bytes +function get_variable_byte_size(hor_grid, z_grid, t_grid, G, num_zlevels) result(var_sz) + character(len=*), intent(in) :: hor_grid !< horizontal grid string + character(len=*), intent(in) :: z_grid !< vertical grid string + character(len=*), intent(in) :: t_grid !< time string + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure; + integer, intent(in) :: num_zlevels !< number of vertical levels + ! local + integer(kind=8) :: var_sz !< The size in bytes of each variable + integer :: var_periods + character(len=8) :: t_grid_read='' + + var_periods = 0 + + if (trim(hor_grid) == '1') then + var_sz = 8 + else + var_sz = 8*(G%Domain%niglobal+1)*(G%Domain%njglobal+1) + endif + + select case (trim(z_grid)) + case ('L') ; var_sz = var_sz * num_zlevels + case ('i') ; var_sz = var_sz * (num_zlevels+1) + end select + + if (adjustl(t_grid(1:1)) == 'p') then + if (len_trim(t_grid(2:8)) > 0) then + var_periods = -1 + t_grid_read = adjustl(t_grid(2:8)) + read(t_grid_read,*) var_periods + if (var_periods > 1) var_sz = var_sz * var_periods + endif + endif + +end function get_variable_byte_size + + end module MOM_restart diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index 1293499930..309a839750 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -17,6 +17,7 @@ module MOM_string_functions public extract_real public remove_spaces public slasher +public append_substring contains @@ -419,6 +420,34 @@ function slasher(dir) endif end function slasher +!> append a string (substring) to another string (string_in) and return the +!! concatenated string (string_out) +function append_substring(string_in, substring) result(string_out) + character(len=*), intent(in) :: string_in !< input string + character(len=*), intent(in) :: substring !< string to append string_in + ! local + character(len=1024) :: string_out + character(len=1024) :: string_joined + integer :: string_in_length + integer :: substring_length + + string_out = '' + string_joined = '' + string_in_length = 0 + substring_length = 0 + + string_in_length = len_trim(string_in) + substring_length = len_trim(substring) + + if (string_in_length > 0) then + if (substring_length > 0) then + string_joined = trim(string_in)//trim(substring) + string_out(1:len_trim(string_joined)) = trim(string_joined) + endif + endif + +end function append_substring + !> \namespace mom_string_functions !! !! By Alistair Adcroft and Robert Hallberg, last updated Sept. 2013. diff --git a/src/framework/MOM_write_field_fms2.F90 b/src/framework/MOM_write_field_fms2.F90 new file mode 100644 index 0000000000..2bfda13c9a --- /dev/null +++ b/src/framework/MOM_write_field_fms2.F90 @@ -0,0 +1,1663 @@ +!> This module contains wrapper functions to write data to netcdf files +module MOM_write_field_fms2 + +! This file is part of MOM6. See LICENSE.md for the license. + + +use MOM_axis, only : MOM_get_diagnostic_axis_data, MOM_register_diagnostic_axis +use MOM_axis, only : axis_data_type, get_time_index, get_var_dimension_metadata +use MOM_axis, only : get_time_units, convert_checksum_to_string +use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING +use MOM_domains, only : MOM_domain_type +use MOM_domains, only : get_simple_array_i_ind, get_simple_array_j_ind +use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_string_functions, only : lowercase, append_substring +use MOM_verticalGrid, only : verticalGrid_type +use mpp_mod, only : mpp_pe, mpp_npes +use mpp_domains_mod, only : domain2d +use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST +use mpp_domains_mod, only : mpp_get_domain_npes, mpp_define_io_domain, mpp_get_io_domain +use netcdf +! fms2_io +use fms2_io_mod, only : check_if_open, get_dimension_size +use fms2_io_mod, only : get_num_dimensions, get_num_variables, get_variable_names +use fms2_io_mod, only : get_unlimited_dimension_name, get_variable_dimension_names +use fms2_io_mod, only : get_variable_num_dimensions, get_variable_size, get_variable_units +use fms2_io_mod, only : get_variable_unlimited_dimension_index, is_dimension_unlimited +use fms2_io_mod, only : is_dimension_registered, register_axis +use fms2_io_mod, only : register_field, register_variable_attribute, fms2_open_file => open_file +use fms2_io_mod, only : fms2_close_file => close_file, write_data, variable_exists +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited + +implicit none; private + +public write_field + +! CAUTION: The following variables are saved by default, and are only necessary for consecutive calls to +! write_field with the same file name. The user should ensure that fms2_close_file on +! the fileobj_write_field structures are called at every requisite time step at after the last +! variable is written to the file by omitting the optional leave_file_open argument, or setting it to .false. + +!> netCDF non-domain-decomposed file object returned by call to open_file in write_field calls +type(FmsNetcdfFile_t), private :: fileobj_write_field + +!> netCDF domain-decomposed file object returned by call to open_file in write_field calls +type(FmsNetcdfDomainFile_t), private :: fileobj_write_field_dd + +!> index of the time_level value that is written to netCDF file by the write_field routines +integer, private :: write_field_time_index + +!> interface to write data to a netcdf file generated by create_file +interface write_field + module procedure write_field_4d_DD + module procedure write_field_3d_DD + module procedure write_field_2d_DD + module procedure write_field_1d_DD + module procedure write_scalar + module procedure write_field_4d_noDD + module procedure write_field_3d_noDD + module procedure write_field_2d_noDD + module procedure write_field_1d_noDD +end interface + +!> interface to apply a scale factor to an array after reading in a field +interface scale_data + module procedure scale_data_4d + module procedure scale_data_3d + module procedure scale_data_2d + module procedure scale_data_1d +end interface + +contains +!> This function uses the fms_io function write_data to write a 1-D domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_1d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is + !! the variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< if .true., leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension(:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: num_dims, substring_index + integer :: dim_unlim_size! size of the unlimited dimension + integer, dimension(1) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=1024) :: filename_temp + character(len=48), dimension(2) :: dim_names !< variable dimension names (or name, in the 1-D case); 1 extra + !! dimension in case appending along the time axis + integer, dimension(2) :: dim_lengths !< variable dimension lengths (or length, in the 1-D case) + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size=0 + dim_unlim_name="" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified + ! and not assumed from the t_grid value + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + ! define the start and edge_length arguments + start(:) = 1 + nwrite(:) = dim_lengths(1) + if (present(start_index)) then + start(1) = max(1, start_index(1)) + endif + + if (present(edge_lengths)) then + nwrite(1) = max(dim_lengths(1),edge_lengths(1)) + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + + if (.not.(check_if_open(fileobj_write_field_dd))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_1d_DD:mode argument must be write, overwrite, or append") + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & + domain%mpp_domain, is_restart=.false.) + ! register the diagnostic axis associated with the variable + call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(1)), dim_lengths(1)) + endif + ! register and write the time_level + if (present(time_level)) then + if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time_level if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & + corner=(/write_field_time_index/), edge_lengths=(/1/)) + endif + endif + ! register the field if it is not already in the file + if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then + call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the variable + if (present(time_level)) then + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) + write_field_time_index = 0 + endif + nullify(data_tmp) +end subroutine write_field_1d_DD + +!> This function uses the fms_io function write_data to write a 2-D domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_2d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is + !! the variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer :: data_tmp(:,:) => null() + integer :: i, is, ie, js, je, j, ndims, num_dims, substring_index + integer, allocatable, dimension(:) :: x_inds, y_inds + integer :: dim_unlim_size ! size of the unlimited dimension + integer :: file_dim_length + integer, dimension(2) :: start, nwrite ! indices for starting points and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(3) :: dim_names ! variable dimension names; 1 extra dimension in case appending + ! along the time axis + character(len=48), allocatable, dimension(:) :: file_dim_names + integer, dimension(3) :: dim_lengths ! variable dimension lengths + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_lengths(:) = 0 + dim_names(:) = "" + dim_unlim_size = 0 + dim_unlim_name = "" + ndims = 2 + num_dims = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension + ! is user-specified rather than derived from the t_grid value + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + start(:) = 1 + nwrite(:) = dim_lengths(1:ndims) + + if (present(start_index)) then + do i=1,ndims + start(i) = max(1,start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i),edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + + if (.not.(check_if_open(fileobj_write_field_dd))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_2d_DD:mode argument must be write, overwrite, or append") + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & + domain%mpp_domain, is_restart=.false.) + endif + ! register the horizontal diagnostic axes associated with the variable + do i=1,num_dims + if (.not.(is_dimension_registered(fileobj_write_field_dd, trim(dim_names(i))))) & + call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) + enddo + ! register and write the time_level + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field_dd, dim_unlim_name) + call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time_level if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & + corner=(/write_field_time_index/), edge_lengths=(/1/)) + endif + endif + ! register the variable if it is not already in the file + if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then + call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the variable + if (present(time_level)) then + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) + write_field_time_index=0 + if (allocated(file_dim_names)) deallocate(file_dim_names) + endif + nullify(data_tmp) +end subroutine write_field_2d_DD + +!> This function uses the fms_io function write_data to write a 3-D domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_3d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is + !! the variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension(:,:,:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: i, is, ie, js, je, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(3) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(4) :: dim_names !< variable dimension names; 1 extra dimension in case appending + !! along the time axis + integer, dimension(4) :: dim_lengths !< variable dimension lengths + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size = 0 + dim_unlim_name = "" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified + ! and not assumed from the t_grid value + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + ndims = 3 + start(:) = 1 + nwrite(:) = dim_lengths(1:3) + if (present(start_index)) then + do i=1,ndims + start(i) = max(1,start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i), edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + ! open the file + if (.not.(check_if_open(fileobj_write_field_dd))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_3d_DD:mode argument must be write, overwrite, or append") + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & + domain%mpp_domain, is_restart=.false.) + ! register the horizontal and vertical diagnostic axes associated with the variable + do i=1,ndims + call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) + enddo + endif + ! register and write the time_level + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field_dd ,dim_unlim_name) + call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time_level if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size ) & + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & + corner=(/write_field_time_index/), edge_lengths=(/1/)) + endif + endif + ! register the field if it is not already in the file + if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then + call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the data + if (present(time_level)) then + call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) + write_field_time_index=0 + endif + nullify(data_tmp) + +end subroutine write_field_3d_DD + +!> This function uses the fms_io function write_data to write a 4-D domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_4d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, t_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + character(len=*), intent(in) :: t_grid !< time descriptor + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is + !! the variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension(:,:,:,:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + real :: file_time ! most recent time currently written to file + integer :: i, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(4) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(4) :: dim_names ! variable dimension names + integer, dimension(4) :: dim_lengths ! variable dimension lengths + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + num_dims = 0 + dim_unlim_size = 0 + dim_unlim_name = "" + dim_names(:) = "" + dim_lengths(:) = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + ndims = 4 + start(:) = 1 + nwrite(:) = dim_lengths(:) + if (present(start_index)) then + do i=1,ndims + start(i) = max(1,start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i), edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + ! open the file + if (.not.(check_if_open(fileobj_write_field_dd))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_4d_DD:mode argument must be write, overwrite, or append") + ! define the io domain for 1-pe jobs because it is required to write domain-decomposed files + if (mpp_get_domain_npes(domain%mpp_domain) .eq. 1 ) then + if (.not. associated(mpp_get_io_domain(domain%mpp_domain))) & + call mpp_define_io_domain(domain%mpp_domain, (/1,1/)) + endif + ! get the index of the corresponding time_level the first time the file is opened + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & + domain%mpp_domain, is_restart=.false.) + ! register the horizontal and vertical diagnostic axes associated with the variable + do i=1,ndims + call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) + enddo + endif + ! register the time dimension and write the time_level + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field_dd, dim_unlim_name) + call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time_level if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & + corner=(/write_field_time_index/), edge_lengths=(/1/)) + endif + endif + ! register the variable if it is not already in the file + if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then + call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the data + if (present(time_level)) then + call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field_dd, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) + write_field_time_index=0 + endif + nullify(data_tmp) +end subroutine write_field_4d_DD + +!> This routine uses the fms_io function write_data to write a scalar variable named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_scalar(filename, fieldname, data, mode, time_level, time_units, leave_file_open, units, longname) + 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, intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=48), dimension(1) :: dim_names ! variable dimension names + integer :: i, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + real, allocatable, dimension(:) :: file_times + integer, dimension(1) :: dim_lengths ! variable dimension lengths + integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file + + dim_unlim_size = 0 + dim_unlim_name= "" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + + if (.not.(check_if_open(fileobj_write_field))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_scaler:mode argument must be write, overwrite, or append") + ! get the index of the corresponding time_level the first time the file is opened + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + if (.not.(allocated(pelist))) then + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), trim(mode), is_restart=.false., & + pelist=pelist) + endif + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field, dim_unlim_name) + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + ! write the time value if it is not already written to the file + if (.not.(variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/)) + else + ! write the next time value if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & + edge_lengths=(/1/)) + endif + endif + ! register the variable + if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then + if (present(time_level)) then + call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=(/trim(dim_unlim_name)/)) + else + call register_field(fileobj_write_field, trim(fieldname), "double") + endif + if (present(units)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) + endif + ! write the data + if (present(time_level)) then + call write_data(fileobj_write_field, trim(fieldname), data, unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field, trim(fieldname), data) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) + if (allocated(pelist)) deallocate(pelist) + write_field_time_index=0 + endif +end subroutine write_scalar + +!> This function uses the fms_io function write_data to write a 1-D non-domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_1d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the + !! variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< if .true., leave file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + real, pointer, dimension(:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: i, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(1) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(2) :: dim_names ! variable dimension names (up to 2 if appended at time level) + integer, dimension(2) :: dim_lengths ! variable dimension lengths + integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size = 0 + dim_unlim_name= "Time" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified + ! and not assumed from the t_grid value. + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + start(:) = 1 + nwrite(:) = dim_lengths(1) + if (present(start_index)) then + start(1) = max(1,start_index(1)) + endif + + if (present(edge_lengths)) then + nwrite(1) = max(dim_lengths(1),edge_lengths(1)) + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + + if (.not.(check_if_open(fileobj_write_field))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_1d_noDD:mode argument must be write, overwrite, or append") + ! get the index of the corresponding time_level the first time the file is opened + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + if (.not.(allocated(pelist))) then + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & + is_restart=.false., pelist=pelist) + endif + ! write the data, and the time value if it is not already written to the file + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time value if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & + edge_lengths=(/1/)) + endif + endif + ! register the field if it is not already in the file + if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then + call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = '' + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the variable to the file + if (present(time_level)) then + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) + if (allocated(pelist)) deallocate(pelist) + write_field_time_index = 0 + endif + nullify(data_tmp) + +end subroutine write_field_1d_noDD + +!> This function uses the fms_io function write_data to write a scalar variable named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_2d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:), intent(in) :: data !< The 2-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(2), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the + !! variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success ! .true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension (:,:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: i, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(2) :: start, nwrite ! indices for starting points and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(3) :: dim_names ! variable dimension names + integer, dimension(3) :: dim_lengths ! variable dimension lengths + integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size = 0 + dim_unlim_name = "" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified + ! and not assumed from the t_grid value + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + + ! set the start (start_index) and nwrite (edge_lengths) values + ndims=2 + start(:) = 1 + nwrite(:) = dim_lengths(1:2) + if (present(start_index)) then + do i=1,ndims + start(i) = max(1,start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i),edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + + if (.not.(check_if_open(fileobj_write_field))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_2d_noDD:mode argument must be write, overwrite, or append") + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + if(.not.(allocated(pelist))) then + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & + is_restart=.false., pelist=pelist) + endif + + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time value if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & + edge_lengths=(/1/)) + endif + endif + + ! register the variable to the file + if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then + call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the variable to the file + if (present(time_level)) then + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) + if (allocated(pelist)) deallocate(pelist) + write_field_time_index=0 + endif + nullify(data_tmp) +end subroutine write_field_2d_noDD + +!> This function uses the fms_io function write_data to write a 3-D non-domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_3d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:,:), intent(in) :: data !< The 3-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + integer, dimension(4), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the + !! variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension(:,:,:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: i, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(3) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time_units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(4) :: dim_names ! variable dimension names + integer, dimension(4) :: dim_lengths ! variable dimension lengths + integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size = 0 + dim_unlim_name = "" + dim_names(:) = "" + dim_lengths(:) = 0 + num_dims = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + ! get the dimension names and lengths + ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified + ! and not assumed from the t_grid value + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, dG=dG) + endif + + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + ndims = 3 + start(:) = 1 + nwrite(:) = dim_lengths(1:3) + if (present(start_index)) then + do i=1,ndims + start(i) = max(1,start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i), edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + ! open the file + if (.not.(check_if_open(fileobj_write_field))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_io:write_3d_noDD:mode argument must be write, overwrite, or append") + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + if (.not.(allocated(pelist))) then + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & + is_restart=.false., pelist=pelist) + endif + ! register and write the time_level + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + + if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + ! write the time_level if it is larger than the most recent file time + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & + edge_lengths=(/1/)) + endif + endif + ! register the field if it is not already in the file + if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then + call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) + endif + endif + + if (present(time_level)) then + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) + if (allocated(pelist)) deallocate(pelist) + write_field_time_index=0 + endif + nullify(data_tmp) +end subroutine write_field_3d_noDD + +!> This function uses the fms_io function write_data to write a 4-D non-domain-decomposed data field named "fieldname" +!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM +!! file write procedure. +subroutine write_field_4d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, t_grid, & + start_index, edge_lengths, time_level, time_units, scale, & + checksums, G, dG, GV, leave_file_open, units, longname) + 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, target, dimension(:,:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data + character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" + character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor + character(len=*), intent(in) :: z_grid !< vertical grid descriptor + character(len=*), intent(in) :: t_grid !< time descriptor + integer, dimension(4), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 + integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the + !! variable size + real, optional, intent(in) :: time_level !< time value to write + real, optional, intent(in) :: time_units !< length of the units for time [s]. The + !! default value is 86400.0, for 1 day. + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied by before they are written. + integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is + !! required if the new file uses any + !! vertical grid axes. + logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open + character(len=*), optional, intent(in) :: units !< variable units + character(len=*), optional, intent(in) :: longname !< long name variable attribute + ! local + logical :: file_open_success !.true. if call to open_file is successful + logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. + real, pointer, dimension(:,:,:,:) :: data_tmp => null() ! enables data to be passed to functions as intent(inout) + integer :: i, ndims, num_dims, substring_index + integer :: dim_unlim_size ! size of the unlimited dimension + integer, dimension(4) :: start, nwrite ! indices for first data value and number of values to write + character(len=20) :: t_units ! time units + character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file + character(len=1024) :: filename_temp + character(len=64) :: checksum_char ! checksum character array created from checksum argument + character(len=48), dimension(4) :: dim_names ! variable dimension names + integer, dimension(4) :: dim_lengths ! variable dimension lengths + integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file + + close_the_file = .true. + if (present(leave_file_open)) close_the_file = .not.(leave_file_open) + + dim_unlim_size = 0 + dim_unlim_name = "" + dim_names(:) = "" + dim_lengths(:) = 0 + ndims = 4 + num_dims = 0 + ! append '.nc' to the file name if it is missing + filename_temp = "" + substring_index = 0 + substring_index = index(trim(filename), ".nc") + if (substring_index <= 0) then + filename_temp = append_substring(filename,".nc") + else + filename_temp = filename + endif + + ! get the dimension names and lengths + if (present(G)) then + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, G=G) + elseif(present(dG)) then + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, dG=dG) + endif + if (present(GV)) & + call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & + dim_lengths, num_dims, GV=GV) + ! set the start (start_index) and nwrite (edge_lengths) values + start(:) = 1 + nwrite(:) = dim_lengths(:) + if (present(start_index)) then + do i=1,ndims + start(i) = max(1, start_index(i)) + enddo + endif + + if (present(edge_lengths)) then + do i=1,ndims + nwrite(i) = max(dim_lengths(i), edge_lengths(i)) + enddo + endif + + data_tmp => data + ! scale the data + if (present(scale)) then ; if (scale /= 1.0) then + call scale_data(data_tmp,scale) + endif ; endif + ! open the file + if (.not.(check_if_open(fileobj_write_field))) then + if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & + (lowercase(trim(mode)) .ne. "overwrite")) & + call MOM_error(FATAL,"MOM_write_field_fms2:write_4d_noDD:mode argument must be write, overwrite, or append") + ! get the time_level index + if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) + ! get the pes associated with the file. + !>\note this is required so that only pe(1) is identified as the root pe to create the file + !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure + if (.not.(allocated(pelist))) then + allocate(pelist(mpp_npes())) + pelist(:) = 0 + do i=1,size(pelist) + pelist(i) = i-1 + enddo + endif + ! open the file in write or append mode + file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & + is_restart=.false., pelist=pelist) + endif + ! register and write the time_level + if (present(time_level)) then + call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) + call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) + num_dims=num_dims+1 + dim_names(num_dims) = trim(dim_unlim_name) + ! write the time value if it is not already written to the file + if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then + ! set the time units + t_units = "" + if (present(time_units)) then + t_units = get_time_units(time_units) + else + t_units = "days" + endif + + call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) + call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) + else + if (write_field_time_index .gt. dim_unlim_size) & + call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & + edge_lengths=(/1/)) + endif + endif + ! register the variable + if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then + call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) + if (present(units)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) + if (present(longname)) & + call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) + ! write the checksum attribute + if (present(checksums)) then + ! convert the checksum to a string + checksum_char = "" + checksum_char = convert_checksum_to_string(checksums(1,1)) + call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) + endif + endif + ! write the variable to the file + if (present(time_level)) then + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite, & + unlim_dim_level=write_field_time_index) + else + call write_data(fileobj_write_field, trim(fieldname), data_tmp, corner=start, edge_lengths=nwrite) + endif + ! close the file + if (close_the_file) then + if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) + deallocate(pelist) + write_field_time_index=0 + endif + nullify(data_tmp) +end subroutine write_field_4d_nodd + +!> apply a scale factor to a 1d array +subroutine scale_data_1d(data, scale_factor) + real, dimension(:), intent(inout) :: data !< The 1-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + + if (scale_factor /= 1.0) then + data(:) = scale_factor*data(:) + endif +end subroutine scale_data_1d + +!> apply a scale factor to a 2d array +subroutine scale_data_2d(data, scale_factor, MOM_domain) + real, dimension(:,:), intent(inout) :: data !< The 2-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je) = scale_factor*data(is:ie,js:je) + else + data(:,:) = scale_factor*data(:,:) + endif + endif +end subroutine scale_data_2d + +!> apply a scale factor to a 3d array +subroutine scale_data_3d(data, scale_factor, MOM_domain) + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:) = scale_factor*data(is:ie,js:je,:) + else + data(:,:,:) = scale_factor*data(:,:,:) + endif + endif +end subroutine scale_data_3d + +!> apply a scale factor to a 4d array +subroutine scale_data_4d(data, scale_factor, MOM_domain) + real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional data array + real, intent(in) :: scale_factor !< Scale factor + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< The domain that describes the decomposition + ! local + integer :: is, ie, js, je + + if (scale_factor /= 1.0) then + if (present(MOM_domain)) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:,:) = scale_factor*data(is:ie,js:je,:,:) + else + data(:,:,:,:) = scale_factor*data(:,:,:,:) + endif + endif +end subroutine scale_data_4d + + +end module MOM_write_field_fms2 diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 66fd873f67..fe9a5bc75f 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1512,8 +1512,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl elseif (.not.new_sim) then ! This line calls a subroutine that reads the initial conditions from a restart file. call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.") - call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & - G, CS%restart_CSp) + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface + call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & + G, CS%restart_CSp, use_fms2=.true.) if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= US%m_to_Z)) then Z_rescale = US%m_to_Z / US%m_to_Z_restart @@ -1587,7 +1588,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (save_IC .and. .not.((dirs%input_filename(1:1) == 'r') .and. & (LEN_TRIM(dirs%input_filename) == 1))) then call save_restart(dirs%output_directory, CS%Time, G, & - CS%restart_CSp, filename=IC_file) + CS%restart_CSp, filename=IC_file, write_ic=.true.) endif @@ -1780,8 +1781,8 @@ subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_su if (present(directory)) then ; restart_dir = directory else ; restart_dir = CS%restart_output_dir ; endif - - call save_restart(restart_dir, Time, CS%grid, CS%restart_CSp, time_stamped) + ! NOTE: first use_fms2=.true. routes routine to fms2 IO interface + call save_restart(restart_dir, Time, CS%grid, CS%restart_CSp, time_stamped, use_fms2=.true.) end subroutine ice_shelf_save_restart diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 9f505325bf..983c008473 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -481,8 +481,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (.not.new_sim) then ! This block restores the state from a restart file. ! This line calls a subroutine that reads the initial conditions ! from a previously generated file. + + ! NOTE: use_fms2=.true. routes routine to fms2 IO interface call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & - G, restart_CS) + G, restart_CS, use_fms2=.true.) if (present(Time_in)) Time = Time_in if ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then H_rescale = GV%m_to_H / GV%m_to_H_restart