diff --git a/config_src/ice_solo_driver/ice_shelf_driver.F90 b/config_src/ice_solo_driver/ice_shelf_driver.F90 index b1323a5485..bd64050a6f 100644 --- a/config_src/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/ice_solo_driver/ice_shelf_driver.F90 @@ -38,9 +38,9 @@ program Shelf_main use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_hor_index, only : hor_index_type, hor_index_init - use MOM_io, only : MOM_io_init, file_exists, open_file, close_file + use MOM_io, only : MOM_io_init, file_exists, open_ASCII_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end - use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE + use MOM_io, only : APPEND_FILE, READONLY_FILE, SINGLE_FILE use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : save_restart use MOM_string_functions,only : uppercase @@ -176,7 +176,7 @@ program Shelf_main if (file_exists('input.nml')) then ! Provide for namelist specification of the run length and calendar data. - call open_file(unit, 'input.nml', form=ASCII_FILE, action=READONLY_FILE) + call open_ASCII_file(unit, 'input.nml', action=READONLY_FILE) read(unit, ice_solo_nml, iostat=io_status) call close_file(unit) ierr = check_nml_error(io_status,'ice_solo_nml') @@ -187,15 +187,14 @@ program Shelf_main ! Read ocean_solo restart, which can override settings from the namelist. if (file_exists(trim(dirs%restart_input_dir)//'ice_solo.res')) then - call open_file(unit,trim(dirs%restart_input_dir)//'ice_solo.res', & - form=ASCII_FILE,action=READONLY_FILE) + call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ice_solo.res', action=READONLY_FILE) read(unit,*) calendar_type read(unit,*) date_init read(unit,*) date call close_file(unit) else calendar = uppercase(calendar) - if (calendar(1:6) == 'JULIAN') then ; calendar_type = JULIAN + if (calendar(1:6) == 'JULIAN') then ; calendar_type = JULIAN elseif (calendar(1:9) == 'GREGORIAN') then ; calendar_type = GREGORIAN elseif (calendar(1:6) == 'NOLEAP') then ; calendar_type = NOLEAP elseif (calendar(1:10)=='THIRTY_DAY') then ; calendar_type = THIRTY_DAY_MONTHS @@ -341,15 +340,14 @@ program Shelf_main call diag_mediator_close_registration(diag) ! Write out a time stamp file. - if (calendar_type /= NO_CALENDAR) then - call open_file(unit, 'time_stamp.out', form=ASCII_FILE, action=APPEND_FILE, & - threading=SINGLE_FILE) + if (is_root_pe() .and. (calendar_type /= NO_CALENDAR)) then + call open_ASCII_file(unit, 'time_stamp.out', action=APPEND_FILE) call get_date(Time, date(1), date(2), date(3), date(4), date(5), date(6)) month = month_name(date(2)) - if (is_root_pe()) write(unit,'(6i4,2x,a3)') date, month(1:3) + write(unit,'(6i4,2x,a3)') date, month(1:3) call get_date(Time_end, date(1), date(2), date(3), date(4), date(5), date(6)) month = month_name(date(2)) - if (is_root_pe()) write(unit,'(6i4,2x,a3)') date, month(1:3) + write(unit,'(6i4,2x,a3)') date, month(1:3) call close_file(unit) endif @@ -428,19 +426,19 @@ program Shelf_main dirs%restart_output_dir) ! Write ice shelf solo restart file. - call open_file(unit, trim(dirs%restart_output_dir)//'shelf.res', nohdrs=.true.) if (is_root_pe())then - write(unit, '(i6,8x,a)') calendar_type, & - '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - - call get_date(Start_time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Model start time: year, month, day, hour, minute, second' - call get_date(Time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Current model time: year, month, day, hour, minute, second' + call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res') + write(unit, '(i6,8x,a)') calendar_type, & + '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' + + call get_date(Start_time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Model start time: year, month, day, hour, minute, second' + call get_date(Time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Current model time: year, month, day, hour, minute, second' + call close_file(unit) endif - call close_file(unit) endif if (is_root_pe()) then diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 584282f27f..9c222bb0bb 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -48,9 +48,9 @@ program MOM_main use MOM_ice_shelf, only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces use MOM_interpolate, only : time_interp_external_init - use MOM_io, only : file_exists, open_file, close_file + use MOM_io, only : file_exists, open_ASCII_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end - use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE + use MOM_io, only : APPEND_FILE, READONLY_FILE use MOM_restart, only : MOM_restart_CS, save_restart use MOM_string_functions,only : uppercase use MOM_surface_forcing, only : set_forcing, forcing_save_restart @@ -238,7 +238,7 @@ program MOM_main if (file_exists('input.nml')) then ! Provide for namelist specification of the run length and calendar data. - call open_file(unit, 'input.nml', form=ASCII_FILE, action=READONLY_FILE) + call open_ASCII_file(unit, 'input.nml', action=READONLY_FILE) read(unit, ocean_solo_nml, iostat=io_status) call close_file(unit) ierr = check_nml_error(io_status,'ocean_solo_nml') @@ -252,15 +252,14 @@ program MOM_main ! Read ocean_solo restart, which can override settings from the namelist. if (file_exists(trim(dirs%restart_input_dir)//'ocean_solo.res')) then - call open_file(unit,trim(dirs%restart_input_dir)//'ocean_solo.res', & - form=ASCII_FILE,action=READONLY_FILE) + call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ocean_solo.res', action=READONLY_FILE) read(unit,*) calendar_type read(unit,*) date_init read(unit,*) date call close_file(unit) else calendar = uppercase(calendar) - if (calendar(1:6) == 'JULIAN') then ; calendar_type = JULIAN + if (calendar(1:6) == 'JULIAN') then ; calendar_type = JULIAN elseif (calendar(1:9) == 'GREGORIAN') then ; calendar_type = GREGORIAN elseif (calendar(1:6) == 'NOLEAP') then ; calendar_type = NOLEAP elseif (calendar(1:10)=='THIRTY_DAY') then ; calendar_type = THIRTY_DAY_MONTHS @@ -432,15 +431,14 @@ program MOM_main call diag_mediator_close_registration(diag) ! Write out a time stamp file. - if (calendar_type /= NO_CALENDAR) then - call open_file(unit, 'time_stamp.out', form=ASCII_FILE, action=APPEND_FILE, & - threading=SINGLE_FILE) + if (is_root_pe() .and. (calendar_type /= NO_CALENDAR)) then + call open_ASCII_file(unit, 'time_stamp.out', action=APPEND_FILE) call get_date(Time, date(1), date(2), date(3), date(4), date(5), date(6)) month = month_name(date(2)) - if (is_root_pe()) write(unit,'(6i4,2x,a3)') date, month(1:3) + write(unit,'(6i4,2x,a3)') date, month(1:3) call get_date(Time_end, date(1), date(2), date(3), date(4), date(5), date(6)) month = month_name(date(2)) - if (is_root_pe()) write(unit,'(6i4,2x,a3)') date, month(1:3) + write(unit,'(6i4,2x,a3)') date, month(1:3) call close_file(unit) endif @@ -618,19 +616,19 @@ program MOM_main if (use_ice_shelf) call ice_shelf_save_restart(ice_shelf_CSp, Time, & dirs%restart_output_dir) ! Write ocean solo restart file. - call open_file(unit, trim(dirs%restart_output_dir)//'ocean_solo.res', nohdrs=.true.) - if (is_root_pe())then - write(unit, '(i6,8x,a)') calendar_type, & - '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - - call get_date(Start_time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Model start time: year, month, day, hour, minute, second' - call get_date(Time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Current model time: year, month, day, hour, minute, second' + if (is_root_pe()) then + call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'ocean_solo.res') + write(unit, '(i6,8x,a)') calendar_type, & + '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' + + call get_date(Start_time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Model start time: year, month, day, hour, minute, second' + call get_date(Time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Current model time: year, month, day, hour, minute, second' + call close_file(unit) endif - call close_file(unit) endif if (is_root_pe()) then diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index fe1563232f..4523f029bf 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -22,7 +22,7 @@ module MOM_ALE use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, param_file_type, log_param use MOM_io, only : vardesc, var_desc, fieldtype, SINGLE_FILE -use MOM_io, only : create_file, write_field, close_file +use MOM_io, only : create_file, write_field, close_file, file_type use MOM_interface_heights,only : find_eta use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S @@ -1273,7 +1273,7 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory ) character(len=240) :: filepath type(vardesc) :: vars(2) type(fieldtype) :: fields(2) - integer :: unit + type(file_type) :: IO_handle ! The I/O handle of the fileset real :: ds(GV%ke), dsi(GV%ke+1) filepath = trim(directory) // trim("Vertical_coordinate") @@ -1287,14 +1287,13 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory ) vars(2) = var_desc('ds_interface', getCoordinateUnits( CS%regridCS ), & 'Layer Center Coordinate Separation','1','i','1') - call create_file(unit, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) - call write_field(unit, fields(1), ds) - call write_field(unit, fields(2), dsi) - call close_file(unit) + call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) + call write_field(IO_handle, fields(1), ds) + call write_field(IO_handle, fields(2), dsi) + call close_file(IO_handle) end subroutine ALE_writeCoordinateFile - !> Set h to coordinate values for fixed coordinate systems subroutine ALE_initThicknessToCoord( CS, G, GV, h ) type(ALE_CS), intent(inout) :: CS !< module control structure diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index 303809ac06..45b08cc799 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -15,8 +15,7 @@ module MOM_PointAccel use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type -use MOM_io, only : open_file -use MOM_io, only : APPEND_FILE, ASCII_FILE, MULTIPLE, SINGLE_FILE +use MOM_io, only : open_ASCII_file, APPEND_FILE, MULTIPLE, SINGLE_FILE use MOM_time_manager, only : time_type, get_time, get_date, set_date, operator(-) use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs @@ -120,8 +119,8 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp ! Open up the file for output if this is the first call. if (CS%u_file < 0) then if (len_trim(CS%u_trunc_file) < 1) return - call open_file(CS%u_file, trim(CS%u_trunc_file), action=APPEND_FILE, & - form=ASCII_FILE, threading=MULTIPLE, fileset=SINGLE_FILE) + call open_ASCII_file(CS%u_file, trim(CS%u_trunc_file), action=APPEND_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) if (CS%u_file < 0) then call MOM_error(NOTE, 'Unable to open file '//trim(CS%u_trunc_file)//'.') return @@ -453,8 +452,8 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp ! Open up the file for output if this is the first call. if (CS%v_file < 0) then if (len_trim(CS%v_trunc_file) < 1) return - call open_file(CS%v_file, trim(CS%v_trunc_file), action=APPEND_FILE, & - form=ASCII_FILE, threading=MULTIPLE, fileset=SINGLE_FILE) + call open_ASCII_file(CS%v_file, trim(CS%v_trunc_file), action=APPEND_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) if (CS%v_file < 0) then call MOM_error(NOTE, 'Unable to open file '//trim(CS%v_trunc_file)//'.') return diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 14a4ceabe3..b2e0275ea8 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -12,10 +12,10 @@ module MOM_sum_output use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_interface_heights, only : find_eta -use MOM_io, only : create_file, fieldtype, flush_file, open_file, reopen_file, stdout +use MOM_io, only : create_file, file_type, fieldtype, flush_file, reopen_file use MOM_io, only : file_exists, slasher, vardesc, var_desc, write_field, get_filename_appendix -use MOM_io, only : field_size, read_variable, read_attribute -use MOM_io, only : APPEND_FILE, ASCII_FILE, SINGLE_FILE, WRITEONLY_FILE +use MOM_io, only : field_size, read_variable, read_attribute, open_ASCII_file, stdout +use MOM_io, only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>) @@ -53,15 +53,15 @@ module MOM_sum_output !> A list of depths and corresponding globally integrated ocean area at each !! depth and the ocean volume below each depth. type :: Depth_List - real :: depth !< A depth [Z ~> m]. - real :: area !< The cross-sectional area of the ocean at that depth [L2 ~> m2]. - real :: vol_below !< The ocean volume below that depth [Z m2 ~> m3]. + integer :: listsize !< length of the list <= niglobal*njglobal + 1 + real, allocatable, dimension(:) :: depth !< A list of depths [Z ~> m] + real, allocatable, dimension(:) :: area !< The cross-sectional area of the ocean at that depth [L2 ~> m2] + real, allocatable, dimension(:) :: vol_below !< The ocean volume below that depth [Z m2 ~> m3] end type Depth_List !> The control structure for the MOM_sum_output module type, public :: sum_output_CS ; private - type(Depth_List), pointer, dimension(:) :: DL => NULL() !< The sorted depth list. - integer :: list_size !< length of sorting vector <= niglobal*njglobal + type(Depth_List) :: DL !< The sorted depth list. integer, allocatable, dimension(:) :: lH !< This saves the entry in DL with a volume just @@ -122,7 +122,7 @@ module MOM_sum_output !! to stdout when the energy files are written. integer :: previous_calls = 0 !< The number of times write_energy has been called. integer :: prev_n = 0 !< The value of n from the last call. - integer :: fileenergy_nc !< NetCDF id of the energy file. + type(file_type) :: fileenergy_nc !< The file handle for the netCDF version of the energy file. integer :: fileenergy_ascii !< The unit number of the ascii version of the energy file. type(fieldtype), dimension(NUM_FIELDS+MAX_FIELDS_) :: & fields !< fieldtype variables for the output fields. @@ -251,9 +251,9 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, & endif allocate(CS%lH(GV%ke)) - call depth_list_setup(G, GV, US, CS) + call depth_list_setup(G, GV, US, CS%DL, CS) else - CS%list_size = 0 + CS%DL%listsize = 1 endif call get_param(param_file, mdl, "TIMEUNIT", Time_unit, & @@ -287,7 +287,8 @@ subroutine MOM_sum_output_end(CS) !! previous call to MOM_sum_output_init. if (associated(CS)) then if (CS%do_APE_calc) then - deallocate(CS%lH, CS%DL) + deallocate(CS%DL%depth, CS%DL%area, CS%DL%vol_below) + deallocate(CS%lH) endif deallocate(CS) @@ -398,8 +399,8 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ integer :: num_nc_fields ! The number of fields that will actually go into ! the NetCDF file. integer :: i, j, k, is, ie, js, je, ns, nz, m, Isq, Ieq, Jsq, Jeq, isr, ier, jsr, jer - integer :: l, lbelow, labove ! indices of deep_area_vol, used to find Z_0APE. - ! lbelow & labove are lower & upper limits for l + integer :: li, lbelow, labove ! indices of deep_area_vol, used to find Z_0APE. + ! lbelow & labove are lower & upper limits for li ! in the search for the entry in lH to use. integer :: start_of_day, num_days real :: reday, var @@ -583,11 +584,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ ! Reopen or create a text output file, with an explanatory header line. if (is_root_pe()) then if (day > CS%Start_time) then - call open_file(CS%fileenergy_ascii, trim(CS%energyfile), & - action=APPEND_FILE, form=ASCII_FILE, nohdrs=.true.) + call open_ASCII_file(CS%fileenergy_ascii, trim(CS%energyfile), action=APPEND_FILE) else - call open_file(CS%fileenergy_ascii, trim(CS%energyfile), & - action=WRITEONLY_FILE, form=ASCII_FILE, nohdrs=.true.) + call open_ASCII_file(CS%fileenergy_ascii, trim(CS%energyfile), action=WRITEONLY_FILE) if (abs(CS%timeunit - 86400.0) < 1.0) then if (CS%use_temperature) then write(CS%fileenergy_ascii,'(" Step,",7x,"Day, Truncs, & @@ -647,23 +646,23 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ lbelow = 1 ; volbelow = 0.0 do k=nz,1,-1 volbelow = volbelow + vol_lay(k) - if ((volbelow >= CS%DL(CS%lH(k))%vol_below) .and. & - (volbelow < CS%DL(CS%lH(k)+1)%vol_below)) then - l = CS%lH(k) + if ((volbelow >= CS%DL%vol_below(CS%lH(k))) .and. & + (volbelow < CS%DL%vol_below(CS%lH(k)+1))) then + li = CS%lH(k) else - labove=CS%list_size+1 - l = (labove + lbelow) / 2 - do while (l > lbelow) - if (volbelow < CS%DL(l)%vol_below) then ; labove = l - else ; lbelow = l ; endif - l = (labove + lbelow) / 2 + labove=CS%DL%listsize + li = (labove + lbelow) / 2 + do while (li > lbelow) + if (volbelow < CS%DL%vol_below(li)) then ; labove = li + else ; lbelow = li ; endif + li = (labove + lbelow) / 2 enddo - CS%lH(k) = l + CS%lH(k) = li endif - lbelow = l - Z_0APE(K) = CS%DL(l)%depth - (volbelow - CS%DL(l)%vol_below) / CS%DL(l)%area + lbelow = li + Z_0APE(K) = CS%DL%depth(li) - (volbelow - CS%DL%vol_below(li)) / CS%DL%area(li) enddo - Z_0APE(nz+1) = CS%DL(2)%depth + Z_0APE(nz+1) = CS%DL%depth(2) ! Calculate the Available Potential Energy integrated over each interface. With a nonlinear ! equation of state or with a bulk mixed layer this calculation is only approximate. @@ -1091,41 +1090,53 @@ end subroutine accumulate_net_input !! cross sectional areas at each depth and the volume of fluid deeper !! than each depth. This might be read from a previously created file !! or it might be created anew. (For now only new creation occurs. -subroutine depth_list_setup(G, GV, US, CS) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(Sum_output_CS), pointer :: CS !< The control structure returned by a - !! previous call to MOM_sum_output_init. +subroutine depth_list_setup(G, GV, US, DL, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes to set up + type(Sum_output_CS), pointer :: CS !< The control structure returned by a + !! previous call to MOM_sum_output_init. ! Local variables + logical :: valid_DL_read integer :: k if (CS%read_depth_list) then if (file_exists(CS%depth_list_file)) then - call read_depth_list(G, US, CS, CS%depth_list_file) + if (CS%update_depth_list_chksum) then + call read_depth_list(G, US, DL, CS%depth_list_file, & + require_chksum=CS%require_depth_list_chksum, file_matches=valid_DL_read) + else + call read_depth_list(G, US, DL, CS%depth_list_file, require_chksum=CS%require_depth_list_chksum) + valid_DL_read = .true. ! Otherwise there would have been a fatal error. + endif else if (is_root_pe()) call MOM_error(WARNING, "depth_list_setup: "// & trim(CS%depth_list_file)//" does not exist. Creating a new file.") - call create_depth_list(G, CS) + valid_DL_read = .false. + endif - call write_depth_list(G, US, CS, CS%depth_list_file, CS%list_size+1) + if (.not.valid_DL_read) then + call create_depth_list(G, DL, CS%D_list_min_inc) + call write_depth_list(G, US, DL, CS%depth_list_file) endif else - call create_depth_list(G, CS) + call create_depth_list(G, DL, CS%D_list_min_inc) endif do k=1,GV%ke - CS%lH(k) = CS%list_size + CS%lH(k) = DL%listsize-1 enddo end subroutine depth_list_setup !> create_depth_list makes an ordered list of depths, along with the cross !! sectional areas at each depth and the volume of fluid deeper than each depth. -subroutine create_depth_list(G, CS) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(Sum_output_CS), pointer :: CS !< The control structure set up in MOM_sum_output_init, - !! in which the ordered depth list is stored. +subroutine create_depth_list(G, DL, min_depth_inc) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes to create + real, intent(in) :: min_depth_inc !< The minimum increment bewteen depths in the list [Z ~> m] + ! Local variables real, dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: & Dlist, & !< The global list of bottom depths [Z ~> m]. @@ -1196,14 +1207,14 @@ subroutine create_depth_list(G, CS) D_list_prev = Dlist(indx2(mls)) list_size = 2 do k=mls-1,1,-1 - if (Dlist(indx2(k)) < D_list_prev-CS%D_list_min_inc) then + if (Dlist(indx2(k)) < D_list_prev-min_depth_inc) then list_size = list_size + 1 D_list_prev = Dlist(indx2(k)) endif enddo - CS%list_size = list_size - allocate(CS%DL(CS%list_size+1)) + DL%listsize = list_size+1 + allocate(DL%depth(DL%listsize), DL%area(DL%listsize), DL%vol_below(DL%listsize)) vol = 0.0 ; area = 0.0 Dprev = Dlist(indx2(mls)) @@ -1218,42 +1229,41 @@ subroutine create_depth_list(G, CS) add_to_list = .false. if ((kl == 0) .or. (k==1)) then add_to_list = .true. - elseif (Dlist(indx2(k-1)) < D_list_prev-CS%D_list_min_inc) then + elseif (Dlist(indx2(k-1)) < D_list_prev-min_depth_inc) then add_to_list = .true. D_list_prev = Dlist(indx2(k-1)) endif if (add_to_list) then kl = kl+1 - CS%DL(kl)%depth = Dlist(i) - CS%DL(kl)%area = area - CS%DL(kl)%vol_below = vol + DL%depth(kl) = Dlist(i) + DL%area(kl) = area + DL%vol_below(kl) = vol endif Dprev = Dlist(i) enddo - do while (kl < list_size) + do while (kl+1 < DL%listsize) ! I don't understand why this is needed... RWH kl = kl+1 - CS%DL(kl)%vol_below = CS%DL(kl-1)%vol_below * 1.000001 - CS%DL(kl)%area = CS%DL(kl-1)%area - CS%DL(kl)%depth = CS%DL(kl-1)%depth + DL%vol_below(kl) = DL%vol_below(kl-1) * 1.000001 + DL%area(kl) = DL%area(kl-1) + DL%depth(kl) = DL%depth(kl-1) enddo - CS%DL(CS%list_size+1)%vol_below = CS%DL(CS%list_size)%vol_below * 1000.0 - CS%DL(CS%list_size+1)%area = CS%DL(CS%list_size)%area - CS%DL(CS%list_size+1)%depth = CS%DL(CS%list_size)%depth + DL%vol_below(DL%listsize) = DL%vol_below(DL%listsize-1) * 1000.0 + DL%area(DL%listsize) = DL%area(DL%listsize-1) + DL%depth(DL%listsize) = DL%depth(DL%listsize-1) end subroutine create_depth_list !> This subroutine writes out the depth list to the specified file. -subroutine write_depth_list(G, US, CS, filename, list_size) +subroutine write_depth_list(G, US, DL, filename) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(Sum_output_CS), pointer :: CS !< The control structure returned by a - !! previous call to MOM_sum_output_init. + type(Depth_List), intent(in) :: DL !< The list of depths, areas and volumes to write character(len=*), intent(in) :: filename !< The path to the depth list file to write. - integer, intent(in) :: list_size !< The size of the depth list. + ! Local variables real, allocatable :: tmp(:) integer :: ncid, dimid(1), Did, Aid, Vid, status, k @@ -1264,7 +1274,7 @@ subroutine write_depth_list(G, US, CS, filename, list_size) if (.not.is_root_pe()) return - allocate(tmp(list_size)) ; tmp(:) = 0.0 + allocate(tmp(DL%listsize)) ; tmp(:) = 0.0 status = NF90_CREATE(filename, 0, ncid) if (status /= NF90_NOERR) then @@ -1272,7 +1282,7 @@ subroutine write_depth_list(G, US, CS, filename, list_size) return endif - status = NF90_DEF_DIM(ncid, "list", list_size, dimid(1)) + status = NF90_DEF_DIM(ncid, "list", DL%listsize, dimid(1)) if (status /= NF90_NOERR) call MOM_error(WARNING, & trim(filename)//trim(NF90_STRERROR(status))) @@ -1319,17 +1329,17 @@ subroutine write_depth_list(G, US, CS, filename, list_size) if (status /= NF90_NOERR) call MOM_error(WARNING, & trim(filename)//trim(NF90_STRERROR(status))) - do k=1,list_size ; tmp(k) = US%Z_to_m*CS%DL(k)%depth ; enddo + do k=1,DL%listsize ; tmp(k) = US%Z_to_m*DL%depth(k) ; enddo status = NF90_PUT_VAR(ncid, Did, tmp) if (status /= NF90_NOERR) call MOM_error(WARNING, & trim(filename)//" depth "//trim(NF90_STRERROR(status))) - do k=1,list_size ; tmp(k) = US%L_to_m**2*CS%DL(k)%area ; enddo + do k=1,DL%listsize ; tmp(k) = US%L_to_m**2*DL%area(k) ; enddo status = NF90_PUT_VAR(ncid, Aid, tmp) if (status /= NF90_NOERR) call MOM_error(WARNING, & trim(filename)//" area "//trim(NF90_STRERROR(status))) - do k=1,list_size ; tmp(k) = US%Z_to_m*US%L_to_m**2*CS%DL(k)%vol_below ; enddo + do k=1,DL%listsize ; tmp(k) = US%Z_to_m*US%L_to_m**2*DL%vol_below(k) ; enddo status = NF90_PUT_VAR(ncid, Vid, tmp) if (status /= NF90_NOERR) call MOM_error(WARNING, & trim(filename)//" vol_below "//trim(NF90_STRERROR(status))) @@ -1340,14 +1350,18 @@ subroutine write_depth_list(G, US, CS, filename, list_size) end subroutine write_depth_list -!> This subroutine reads in the depth list to the specified file -!! and allocates and sets up CS%DL and CS%list_size . -subroutine read_depth_list(G, US, CS, filename) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(Sum_output_CS), pointer :: CS !< The control structure returned by a - !! previous call to MOM_sum_output_init. - character(len=*), intent(in) :: filename !< The path to the depth list file to read. +!> This subroutine reads in the depth list from the specified file +!! and allocates the memory within and sets up DL. +subroutine read_depth_list(G, US, DL, filename, require_chksum, file_matches) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(Depth_List), intent(inout) :: DL !< The list of depths, areas and volumes + character(len=*), intent(in) :: filename !< The path to the depth list file to read. + logical, intent(in) :: require_chksum !< If true, missing or mismatched depth + !! and area checksums result in a fatal error. + logical, optional, intent(out) :: file_matches !< If present, this indicates whether the file + !! has been read with matching depth and area checksums + ! Local variables character(len=240) :: var_msg real, allocatable :: tmp(:) @@ -1361,30 +1375,27 @@ subroutine read_depth_list(G, US, CS, filename) call read_attribute(filename, area_chksum_attr, area_file_chksum, found=area_att_found) if ((.not.depth_att_found) .or. (.not.area_att_found)) then - var_msg = trim(CS%depth_list_file) // " checksums are missing;" - if (CS%require_depth_list_chksum) then + var_msg = trim(filename) // " checksums are missing;" + if (require_chksum) then call MOM_error(FATAL, trim(var_msg) // " aborting.") - elseif (CS%update_depth_list_chksum) then + elseif (present(file_matches)) then call MOM_error(WARNING, trim(var_msg) // " updating file.") - call create_depth_list(G, CS) - call write_depth_list(G, US, CS, CS%depth_list_file, CS%list_size+1) + file_matches = .false. return else - call MOM_error(WARNING, & - trim(var_msg) // " some diagnostics may not be reproducible.") + call MOM_error(WARNING, trim(var_msg) // " some diagnostics may not be reproducible.") endif else call get_depth_list_checksums(G, depth_grid_chksum, area_grid_chksum) if ((trim(depth_grid_chksum) /= trim(depth_file_chksum)) .or. & (trim(area_grid_chksum) /= trim(area_file_chksum)) ) then - var_msg = trim(CS%depth_list_file) // " checksums do not match;" - if (CS%require_depth_list_chksum) then + var_msg = trim(filename) // " checksums do not match;" + if (require_chksum) then call MOM_error(FATAL, trim(var_msg) // " aborting.") - elseif (CS%update_depth_list_chksum) then + elseif (present(file_matches)) then call MOM_error(WARNING, trim(var_msg) // " updating file.") - call create_depth_list(G, CS) - call write_depth_list(G, US, CS, CS%depth_list_file, CS%list_size+1) + file_matches = .false. return else call MOM_error(WARNING, trim(var_msg) // " some diagnostics may not be reproducible.") @@ -1400,20 +1411,14 @@ subroutine read_depth_list(G, US, CS, filename) trim(filename)//" has too many or too few dimensions.") list_size = sizes(1) - CS%list_size = list_size-1 - allocate(CS%DL(list_size)) - allocate(tmp(list_size)) - - call read_variable(filename, "depth", tmp) - do k=1,list_size ; CS%DL(k)%depth = US%m_to_Z*tmp(k) ; enddo - - call read_variable(filename, "area", tmp) - do k=1,list_size ; CS%DL(k)%area = US%m_to_L**2*tmp(k) ; enddo + DL%listsize = list_size + allocate(DL%depth(list_size), DL%area(list_size), DL%vol_below(list_size)) - call read_variable(filename, "vol_below", tmp) - do k=1,list_size ; CS%DL(k)%vol_below = US%m_to_Z*US%m_to_L**2*tmp(k) ; enddo + call read_variable(filename, "depth", DL%depth, scale=US%m_to_Z) + call read_variable(filename, "area", DL%area, scale=US%m_to_L**2) + call read_variable(filename, "vol_below", DL%vol_below, scale=US%m_to_Z*US%m_to_L**2) - deallocate(tmp) + if (present(file_matches)) file_matches = .true. end subroutine read_depth_list diff --git a/src/framework/MOM_array_transform.F90 b/src/framework/MOM_array_transform.F90 index 179bd6550e..d524f618a3 100644 --- a/src/framework/MOM_array_transform.F90 +++ b/src/framework/MOM_array_transform.F90 @@ -13,9 +13,8 @@ module MOM_array_transform -implicit none +implicit none ; private -private public rotate_array public rotate_array_pair public rotate_vector diff --git a/src/framework/MOM_domain_infra.F90 b/src/framework/MOM_domain_infra.F90 index 68385b1c6d..86e85e60a6 100644 --- a/src/framework/MOM_domain_infra.F90 +++ b/src/framework/MOM_domain_infra.F90 @@ -39,7 +39,7 @@ module MOM_domain_infra ! These interfaces are actually implemented or have explicit interfaces in this file. public :: create_MOM_domain, clone_MOM_domain, get_domain_components, get_domain_extent public :: deallocate_MOM_domain, get_global_shape, compute_block_extent -public :: pass_var, pass_vector, fill_symmetric_edges +public :: pass_var, pass_vector, fill_symmetric_edges, rescale_comp_data public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete public :: create_group_pass, do_group_pass, start_group_pass, complete_group_pass public :: redistribute_array, broadcast_domain, global_field @@ -98,6 +98,11 @@ module MOM_domain_infra ! module procedure fill_scalar_symmetric_edges_2d, fill_scalar_symmetric_edges_3d end interface fill_symmetric_edges +!> Rescale the values of an array in its computational domain by a constant factor +interface rescale_comp_data + module procedure rescale_comp_data_4d, rescale_comp_data_3d, rescale_comp_data_2d +end interface rescale_comp_data + !> Pass an array from one MOM domain to another interface redistribute_array module procedure redistribute_array_3d, redistribute_array_2d @@ -1228,6 +1233,57 @@ subroutine redistribute_array_3d(Domain1, array1, Domain2, array2, complete) end subroutine redistribute_array_3d +!> Rescale the values of a 4-D array in its computational domain by a constant factor +subroutine rescale_comp_data_4d(domain, array, scale) + type(MOM_domain_type), intent(in) :: domain !< MOM domain from which to extract information + real, dimension(:,:,:,:), intent(inout) :: array !< The array which is having the data in its + !! computational domain rescaled + real, intent(in) :: scale !< A scaling factor by which to multiply the + !! values in the computational domain of array + integer :: is, ie, js, je + + if (scale == 1.0) return + + call get_simple_array_i_ind(domain, size(array,1), is, ie) + call get_simple_array_j_ind(domain, size(array,2), js, je) + array(is:ie,js:je,:,:) = scale*array(is:ie,js:je,:,:) + +end subroutine rescale_comp_data_4d + +!> Rescale the values of a 3-D array in its computational domain by a constant factor +subroutine rescale_comp_data_3d(domain, array, scale) + type(MOM_domain_type), intent(in) :: domain !< MOM domain from which to extract information + real, dimension(:,:,:), intent(inout) :: array !< The array which is having the data in its + !! computational domain rescaled + real, intent(in) :: scale !< A scaling factor by which to multiply the + !! values in the computational domain of array + integer :: is, ie, js, je + + if (scale == 1.0) return + + call get_simple_array_i_ind(domain, size(array,1), is, ie) + call get_simple_array_j_ind(domain, size(array,2), js, je) + array(is:ie,js:je,:) = scale*array(is:ie,js:je,:) + +end subroutine rescale_comp_data_3d + +!> Rescale the values of a 2-D array in its computational domain by a constant factor +subroutine rescale_comp_data_2d(domain, array, scale) + type(MOM_domain_type), intent(in) :: domain !< MOM domain from which to extract information + real, dimension(:,:), intent(inout) :: array !< The array which is having the data in its + !! computational domain rescaled + real, intent(in) :: scale !< A scaling factor by which to multiply the + !! values in the computational domain of array + integer :: is, ie, js, je + + if (scale == 1.0) return + + call get_simple_array_i_ind(domain, size(array,1), is, ie) + call get_simple_array_j_ind(domain, size(array,2), js, je) + array(is:ie,js:je) = scale*array(is:ie,js:je) + +end subroutine rescale_comp_data_2d + !> create_MOM_domain creates and initializes a MOM_domain_type variables, based on the information !! provided in arguments. subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar_N, layout, io_layout, & @@ -1743,7 +1799,6 @@ subroutine get_domain_extent_d2D(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed) end subroutine get_domain_extent_d2D - !> Return the (potentially symmetric) computational domain i-bounds for an array !! passed without index specifications (i.e. indices start at 1) based on an array size. subroutine get_simple_array_i_ind(domain, size, is, ie, symmetric) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 88415c6782..b25a934b97 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -15,7 +15,7 @@ module MOM_domains use MOM_domain_infra, only : pass_vector_start, pass_vector_complete use MOM_domain_infra, only : create_group_pass, do_group_pass use MOM_domain_infra, only : start_group_pass, complete_group_pass -use MOM_domain_infra, only : global_field, redistribute_array, broadcast_domain +use MOM_domain_infra, only : rescale_comp_data, global_field, redistribute_array, broadcast_domain use MOM_domain_infra, only : MOM_thread_affinity_set, set_MOM_thread_affinity use MOM_domain_infra, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM use MOM_domain_infra, only : CORNER, CENTER, NORTH_FACE, EAST_FACE @@ -28,22 +28,24 @@ module MOM_domains implicit none ; private public :: MOM_infra_init, MOM_infra_end -! Domain types and creation and destruction routines +! Domain types and creation and destruction routines public :: MOM_domain_type, domain2D, domain1D public :: MOM_domains_init, create_MOM_domain, clone_MOM_domain, deallocate_MOM_domain public :: MOM_thread_affinity_set, set_MOM_thread_affinity -! Domain query routines +! Domain query routines public :: get_domain_extent, get_domain_components, compute_block_extent, get_global_shape public :: PE_here, root_PE, num_PEs -! Single call communication routines +! Single call communication routines public :: pass_var, pass_vector, fill_symmetric_edges, broadcast -! Non-blocking communication routines +! Non-blocking communication routines public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete -! Multi-variable group communication routines and type +! Multi-variable group communication routines and type public :: create_group_pass, do_group_pass, group_pass_type, start_group_pass, complete_group_pass -! Global reduction routines +! Global reduction routines public :: sum_across_PEs, min_across_PEs, max_across_PEs public :: global_field, redistribute_array, broadcast_domain +! Simple index-convention-invariant array manipulation routine +public :: rescale_comp_data !> These encoding constants are used to indicate the staggering of scalars and vectors public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR !> These encoding constants are used to indicate the discretization position of a variable diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 050f24db27..1159ac87e1 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -5,7 +5,7 @@ module MOM_io use MOM_array_transform, only : allocate_rotated_array, rotate_array use MOM_domains, only : MOM_domain_type, domain1D, broadcast, get_domain_components -use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE +use MOM_domains, only : rescale_comp_data, AGRID, BGRID_NE, CGRID_NE use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_ensemble_manager, only : get_ensemble_id use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING, is_root_PE @@ -13,10 +13,10 @@ module MOM_io use MOM_grid, only : ocean_grid_type use MOM_io_infra, only : MOM_read_data, MOM_read_vector, read_field_chksum use MOM_io_infra, only : read_data=>MOM_read_data ! read_data will be removed soon. -use MOM_io_infra, only : file_exists, get_file_info, get_file_fields, get_field_atts -use MOM_io_infra, only : open_file, close_file, get_field_size, fieldtype, field_exists -use MOM_io_infra, only : flush_file, get_filename_suffix -use MOM_io_infra, only : get_file_times, axistype, get_axis_data +use MOM_io_infra, only : file_type, file_exists, get_file_info, get_file_fields +use MOM_io_infra, only : open_file, open_ASCII_file, close_file, flush_file, file_is_open +use MOM_io_infra, only : get_field_size, fieldtype, field_exists, get_field_atts +use MOM_io_infra, only : get_file_times, axistype, get_axis_data, get_filename_suffix use MOM_io_infra, only : write_field, write_metadata, write_version use MOM_io_infra, only : MOM_namelist_file, check_namelist_error, io_infra_init, io_infra_end use MOM_io_infra, only : APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE @@ -40,8 +40,8 @@ module MOM_io public :: get_var_sizes, verify_variable_units, num_timelevels, read_variable, read_attribute public :: open_file_to_read, close_file_to_read ! The following are simple pass throughs of routines from MOM_io_infra or other modules. -public :: file_exists, open_file, close_file, flush_file, get_filename_appendix -public :: get_file_info, field_exists, get_file_fields, get_file_times +public :: file_exists, open_file, open_ASCII_file, close_file, flush_file, file_type +public :: get_file_info, field_exists, get_file_fields, get_file_times, get_filename_appendix public :: fieldtype, field_size, get_field_atts public :: axistype, get_axis_data public :: MOM_read_data, MOM_read_vector, read_field_chksum @@ -101,12 +101,12 @@ module MOM_io contains -!> Routine creates a new NetCDF file. It also sets up +!> Routine creates a new NetCDF file. It also sets up fieldtype !! 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) - integer, intent(out) :: unit !< unit id of an open file or -1 on a - !! nonwriting PE with single file output +!! later be written to this file. +subroutine create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, G, dG, GV, checksums) + type(file_type), intent(inout) :: IO_handle !< Handle for a file or fileset that is to be + !! opened or reopened for writing character(len=*), intent(in) :: filename !< full path to the file to create type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename integer, intent(in) :: novars !< number of fields written to filename @@ -173,9 +173,9 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, threading=thread) + call open_file(IO_handle, filename, OVERWRITE_FILE, threading=thread) else - call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, MOM_domain=Domain) + call open_file(IO_handle, filename, OVERWRITE_FILE, MOM_domain=Domain) endif ! Define the coordinates. @@ -245,28 +245,28 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit "create_file: A vertical grid type is required to create a file with a vertical coordinate.") if (use_lath) & - call write_metadata(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & + call write_metadata(IO_handle, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & cartesian='Y', domain=y_domain, data=gridLatT(jsg:jeg)) if (use_lonh) & - call write_metadata(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & + call write_metadata(IO_handle, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & cartesian='X', domain=x_domain, data=gridLonT(isg:ieg)) if (use_latq) & - call write_metadata(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & + call write_metadata(IO_handle, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & cartesian='Y', domain=y_domain, data=gridLatB(JsgB:JegB)) if (use_lonq) & - call write_metadata(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & + call write_metadata(IO_handle, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & cartesian='X', domain=x_domain, data=gridLonB(IsgB:IegB)) if (use_layer) & - call write_metadata(unit, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & + call write_metadata(IO_handle, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & longname="Layer "//trim(GV%zAxisLongName), cartesian='Z', & sense=1, data=GV%sLayer(1:GV%ke)) if (use_int) & - call write_metadata(unit, axis_int, name="Interface", units=trim(GV%zAxisUnits), & + call write_metadata(IO_handle, axis_int, name="Interface", units=trim(GV%zAxisUnits), & longname="Interface "//trim(GV%zAxisLongName), cartesian='Z', & sense=1, data=GV%sInterface(1:GV%ke+1)) @@ -286,9 +286,9 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit write(time_units,'(es8.2," s")') timeunit endif - call write_metadata(unit, axis_time, name="Time", units=time_units, longname="Time", cartesian='T') + call write_metadata(IO_handle, axis_time, name="Time", units=time_units, longname="Time", cartesian='T') else - call write_metadata(unit, axis_time, name="Time", units="days", longname="Time", cartesian= 'T') + call write_metadata(IO_handle, axis_time, name="Time", units="days", longname="Time", cartesian= 'T') endif ; endif if (use_periodic) then @@ -297,7 +297,7 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit ! Define a periodic axis with unit labels. allocate(period_val(num_periods)) do k=1,num_periods ; period_val(k) = real(k) ; enddo - call write_metadata(unit, axis_periodic, name="Period", units="nondimensional", & + call write_metadata(IO_handle, axis_periodic, name="Period", units="nondimensional", & longname="Periods for cyclical varaiables", cartesian='T', data=period_val) deallocate(period_val) endif @@ -338,21 +338,21 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit pack = 1 if (present(checksums)) then - call write_metadata(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + call write_metadata(IO_handle, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & vars(k)%longname, pack=pack, checksum=checksums(k,:)) else - call write_metadata(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + call write_metadata(IO_handle, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & vars(k)%longname, pack=pack) endif enddo - if (use_lath) call write_field(unit, axis_lath) - if (use_latq) call write_field(unit, axis_latq) - if (use_lonh) call write_field(unit, axis_lonh) - if (use_lonq) call write_field(unit, axis_lonq) - if (use_layer) call write_field(unit, axis_layer) - if (use_int) call write_field(unit, axis_int) - if (use_periodic) call write_field(unit, axis_periodic) + if (use_lath) call write_field(IO_handle, axis_lath) + if (use_latq) call write_field(IO_handle, axis_latq) + if (use_lonh) call write_field(IO_handle, axis_lonh) + if (use_lonq) call write_field(IO_handle, axis_lonq) + if (use_layer) call write_field(IO_handle, axis_layer) + if (use_int) call write_field(IO_handle, axis_int) + if (use_periodic) call write_field(IO_handle, axis_periodic) end subroutine create_file @@ -361,9 +361,9 @@ end subroutine create_file !! does not find the file, a new file is created. It also sets up !! structures that describe this file and the variables that will !! later be written to this file. -subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV) - integer, intent(out) :: unit !< unit id of an open file or -1 on a - !! nonwriting PE with single file output +subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, timeunit, G, dG, GV) + type(file_type), intent(inout) :: IO_handle !< Handle for a file or fileset that is to be + !! opened or reopened for writing character(len=*), intent(in) :: filename !< full path to the file to create type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename integer, intent(in) :: novars !< number of fields written to filename @@ -397,7 +397,7 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit inquire(file=check_name,EXIST=exists) if (.not.exists) then - call create_file(unit, filename, vars, novars, fields, threading, timeunit, & + call create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, & G=G, dG=dG, GV=GV) else @@ -412,26 +412,26 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, threading=thread) + call open_file(IO_handle, filename, APPEND_FILE, threading=thread) else - call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, MOM_domain=Domain) + call open_file(IO_handle, filename, APPEND_FILE, MOM_domain=Domain) endif - if (unit < 0) return + if (.not.file_is_open(IO_handle)) return - call get_file_info(unit, ndim, nvar, natt, ntime) + call get_file_info(IO_handle, ndim, nvar, natt, ntime) if (nvar == -1) then write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,& " variables. Clobbering and creating file with ",novars," instead." call MOM_error(WARNING,"MOM_io: "//mesg) - call create_file(unit, filename, vars, novars, fields, threading, timeunit, G=G, GV=GV) + call create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, G=G, GV=GV) elseif (nvar /= novars) then write (mesg,*) "Reopening file ",trim(filename)," with ",novars,& " variables instead of ",nvar,"." call MOM_error(FATAL,"MOM_io: "//mesg) endif - if (nvar > 0) call get_file_fields(unit, fields(1:nvar)) + if (nvar > 0) call get_file_fields(IO_handle, fields(1:nvar)) ! Check for inconsistent field names... ! do i=1,nvar @@ -608,12 +608,14 @@ end subroutine read_var_sizes !> Read a real scalar variable from a netCDF file with the root PE, and broadcast the !! results to all the other PEs. -subroutine read_variable_0d(filename, varname, var, ncid_in) +subroutine read_variable_0d(filename, varname, var, ncid_in, scale) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: varname !< The variable name of the data in the file real, intent(inout) :: var !< The scalar into which to read the data integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the - !! file is opened and closed within this routine. + !! file is opened and closed within this routine + real, optional, intent(in) :: scale !< A scaling factor that the variable is + !! multiplied by before it is returned integer :: varid, ncid, rc character(len=256) :: hdr @@ -634,6 +636,8 @@ subroutine read_variable_0d(filename, varname, var, ncid_in) " Difficulties reading "//trim(varname)//" from "//trim(filename)) if (.not.present(ncid_in)) call close_file_to_read(ncid, filename) + + if (present(scale)) var = scale * var endif call broadcast(var, blocking=.true.) @@ -641,12 +645,14 @@ end subroutine read_variable_0d !> Read a 1-d real variable from a netCDF file with the root PE, and broadcast the !! results to all the other PEs. -subroutine read_variable_1d(filename, varname, var, ncid_in) +subroutine read_variable_1d(filename, varname, var, ncid_in, scale) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: varname !< The variable name of the data in the file real, dimension(:), intent(inout) :: var !< The 1-d array into which to read the data integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the - !! file is opened and closed within this routine. + !! file is opened and closed within this routine + real, optional, intent(in) :: scale !< A scaling factor that the variable is + !! multiplied by before it is returned integer :: varid, ncid, rc character(len=256) :: hdr @@ -667,6 +673,10 @@ subroutine read_variable_1d(filename, varname, var, ncid_in) " Difficulties reading "//trim(varname)//" from "//trim(filename)) if (.not.present(ncid_in)) call close_file_to_read(ncid, filename) + + if (present(scale)) then ; if (scale /= 1.0) then + var(:) = scale * var(:) + endif ; endif endif call broadcast(var, size(var), blocking=.true.) @@ -1331,9 +1341,9 @@ end subroutine query_vardesc !> Write a 4d field to an output file, potentially with rotation -subroutine MOM_write_field_4d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & - fill_value, turns) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine MOM_write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns, scale) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Unrotated field to write @@ -1341,28 +1351,34 @@ subroutine MOM_write_field_4d(io_unit, field_md, MOM_domain, field, tstamp, tile integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written - real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units + real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units or rescaled + real :: scale_fac ! A scaling factor to use before writing the array integer :: qturns ! The number of quarter turns through which to rotate field + integer :: is, ie, js, je ! The extent of the computational domain qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale - if (qturns == 0) then - call write_field(io_unit, field_md, MOM_domain, field, tstamp=tstamp, & + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) else call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, MOM_domain, field_rot, tstamp=tstamp, & + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif end subroutine MOM_write_field_4d !> Write a 3d field to an output file, potentially with rotation -subroutine MOM_write_field_3d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & - fill_value, turns) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine MOM_write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns, scale) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Unrotated field to write @@ -1370,28 +1386,34 @@ subroutine MOM_write_field_3d(io_unit, field_md, MOM_domain, field, tstamp, tile integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written - real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units + real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units or rescaled + real :: scale_fac ! A scaling factor to use before writing the array integer :: qturns ! The number of quarter turns through which to rotate field + integer :: is, ie, js, je ! The extent of the computational domain qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale - if (qturns == 0) then - call write_field(io_unit, field_md, MOM_domain, field, tstamp=tstamp, & + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) else call allocate_rotated_array(field, [1,1,1], qturns, field_rot) call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, MOM_domain, field_rot, tstamp=tstamp, & + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif end subroutine MOM_write_field_3d !> Write a 2d field to an output file, potentially with rotation -subroutine MOM_write_field_2d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & - fill_value, turns) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine MOM_write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns, scale) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Unrotated field to write @@ -1399,45 +1421,75 @@ subroutine MOM_write_field_2d(io_unit, field_md, MOM_domain, field, tstamp, tile integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units + real :: scale_fac ! A scaling factor to use before writing the array integer :: qturns ! The number of quarter turns through which to rotate field + integer :: is, ie, js, je ! The extent of the computational domain - qturns = 0 - if (present(turns)) qturns = modulo(turns, 4) + qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale - if (qturns == 0) then - call write_field(io_unit, field_md, MOM_domain, field, tstamp=tstamp, & + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) else call allocate_rotated_array(field, [1,1], qturns, field_rot) call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, MOM_domain, field_rot, tstamp=tstamp, & + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, & tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif end subroutine MOM_write_field_2d !> Write a 1d field to an output file -subroutine MOM_write_field_1d(io_unit, field_md, field, tstamp, fill_value) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine MOM_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, scale) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp real, optional, intent(in) :: fill_value !< Missing data fill value + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + + real, dimension(:), allocatable :: array ! A rescaled copy of field + real :: scale_fac ! A scaling factor to use before writing the array + integer :: i - call write_field(io_unit, field_md, field, tstamp=tstamp) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + + if (scale_fac == 1.0) then + call write_field(IO_handle, field_md, field, tstamp=tstamp) + else + allocate(array(size(field))) + array(:) = scale_fac * field(:) + if (present(fill_value)) then + do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo + endif + call write_field(IO_handle, field_md, array, tstamp=tstamp) + deallocate(array) + endif end subroutine MOM_write_field_1d !> Write a 0d field to an output file -subroutine MOM_write_field_0d(io_unit, field_md, field, tstamp, fill_value) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine MOM_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, scale) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp real, optional, intent(in) :: fill_value !< Missing data fill value + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + real :: scaled_val ! A rescaled copy of field + + scaled_val = field + if (present(scale)) scaled_val = scale*field + if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif - call write_field(io_unit, field_md, field, tstamp=tstamp) + call write_field(IO_handle, field_md, scaled_val, tstamp=tstamp) end subroutine MOM_write_field_0d !> Given filename and fieldname, this subroutine returns the size of the field in the file diff --git a/src/framework/MOM_io_infra.F90 b/src/framework/MOM_io_infra.F90 index 978ca71c8a..3ea201235a 100644 --- a/src/framework/MOM_io_infra.F90 +++ b/src/framework/MOM_io_infra.F90 @@ -3,8 +3,7 @@ module MOM_io_infra ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_domain_infra, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE -use MOM_domain_infra, only : get_simple_array_i_ind, get_simple_array_j_ind +use MOM_domain_infra, only : MOM_domain_type, rescale_comp_data, AGRID, BGRID_NE, CGRID_NE use MOM_domain_infra, only : domain2d, domain1d, CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING @@ -19,17 +18,17 @@ module MOM_io_infra use mpp_io_mod, only : mpp_get_info, mpp_get_times use mpp_io_mod, only : mpp_io_init ! These are encoding constants. -use mpp_io_mod, only : APPEND_FILE=>MPP_APPEND, ASCII_FILE=>MPP_ASCII -use mpp_io_mod, only : MULTIPLE=>MPP_MULTI, NETCDF_FILE=>MPP_NETCDF +use mpp_io_mod, only : APPEND_FILE=>MPP_APPEND, WRITEONLY_FILE=>MPP_WRONLY use mpp_io_mod, only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY -use mpp_io_mod, only : SINGLE_FILE=>MPP_SINGLE, WRITEONLY_FILE=>MPP_WRONLY +use mpp_io_mod, only : NETCDF_FILE=>MPP_NETCDF, ASCII_FILE=>MPP_ASCII +use mpp_io_mod, only : MULTIPLE=>MPP_MULTI, SINGLE_FILE=>MPP_SINGLE use iso_fortran_env, only : int64 implicit none ; private ! These interfaces are actually implemented or have explicit interfaces in this file. -public :: open_file, close_file, flush_file, file_exists, get_filename_suffix -public :: get_file_info, get_file_fields, get_file_times +public :: open_file, open_ASCII_file, file_is_open, close_file, flush_file, file_exists +public :: get_file_info, get_file_fields, get_file_times, get_filename_suffix public :: MOM_read_data, MOM_read_vector, write_metadata, write_field public :: field_exists, get_field_atts, get_field_size, get_axis_data, read_field_chksum public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version @@ -37,8 +36,8 @@ module MOM_io_infra ! information about fields and axes, respectively, and are opaque to this module. public :: fieldtype, axistype ! These are encoding constant parmeters. -public :: APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE -public :: READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE +public :: ASCII_FILE, NETCDF_FILE, SINGLE_FILE, MULTIPLE +public :: APPEND_FILE, READONLY_FILE, OVERWRITE_FILE, WRITEONLY_FILE public :: CENTER, CORNER, NORTH_FACE, EAST_FACE !> Indicate whether a file exists, perhaps with domain decomposition @@ -47,6 +46,11 @@ module MOM_io_infra module procedure MOM_file_exists end interface +!> Open a file (or fileset) for parallel or single-file I/). +interface open_file + module procedure open_file_type, open_file_unit +end interface open_file + !> Read a data field from a file interface MOM_read_data module procedure MOM_read_data_4d @@ -77,6 +81,25 @@ module MOM_io_infra module procedure write_metadata_axis, write_metadata_field end interface write_metadata +!> Close a file (or fileset). If the file handle does not point to an open file, +!! close_file simply returns without doing anything. +interface close_file + module procedure close_file_type, close_file_unit +end interface close_file + +!> Ensure that the output stream associated with a file handle is fully sent to disk +interface flush_file + module procedure flush_file_type, flush_file_unit +end interface flush_file + +!> Type for holding a handle to an open file and related information +type, public :: file_type ; private + integer :: unit = -1 !< The framework identfier or netCDF unit number of an output file + character(len=:), allocatable :: filename !< The path to this file, if it is open + logical :: open_to_read = .false. !< If true, this file or fileset can be read + logical :: open_to_write = .false. !< If true, this file or fileset can be written to +end type file_type + contains !> Reads the checksum value for a field that was recorded in a file, along with a flag indicating @@ -99,47 +122,67 @@ subroutine read_field_chksum(field, chksum, valid_chksum) end subroutine read_field_chksum !> Returns true if the named file or its domain-decomposed variant exists. -function MOM_file_exists(filename, MOM_Domain) +logical function MOM_file_exists(filename, MOM_Domain) character(len=*), intent(in) :: filename !< The name of the file being inquired about type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition ! This function uses the fms_io function file_exist to determine whether ! a named file (or its decomposed variant) exists. - logical :: MOM_file_exists - MOM_file_exists = file_exist(filename, MOM_Domain%mpp_domain) end function MOM_file_exists !> Returns true if the named file or its domain-decomposed variant exists. -function FMS_file_exists(filename, domain, no_domain) +logical function FMS_file_exists(filename, domain, no_domain) character(len=*), intent(in) :: filename !< The name of the file being inquired about type(domain2d), optional, intent(in) :: domain !< The mpp domain2d that describes the decomposition logical, optional, intent(in) :: no_domain !< This file does not use domain decomposition ! This function uses the fms_io function file_exist to determine whether ! a named file (or its decomposed variant) exists. - logical :: FMS_file_exists - FMS_file_exists = file_exist(filename, domain, no_domain) end function FMS_file_exists -!> close_file closes a file (or fileset). If the file handle does not point to an open file, -!! close_file simply returns without doing anything. -subroutine close_file(unit) - integer, intent(out) :: unit !< The I/O unit for the file to be closed +!> indicates whether an I/O handle is attached to an open file +logical function file_is_open(IO_handle) + type(file_type), intent(in) :: IO_handle !< Handle to a file to inquire about + + file_is_open = (IO_handle%unit >= 0) +end function file_is_open + +!> closes a file (or fileset). If the file handle does not point to an open file, +!! close_file_type simply returns without doing anything. +subroutine close_file_type(IO_handle) + type(file_type), intent(inout) :: IO_handle !< The I/O handle for the file to be closed + + call mpp_close(IO_handle%unit) + if (allocated(IO_handle%filename)) deallocate(IO_handle%filename) + IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .false. +end subroutine close_file_type + +!> closes a file. If the unit does not point to an open file, +!! close_file_unit simply returns without doing anything. +subroutine close_file_unit(unit) + integer, intent(inout) :: unit !< The I/O unit for the file to be closed call mpp_close(unit) -end subroutine close_file +end subroutine close_file_unit -!> Ensure that the output stream associated with a unit is fully sent to dis. -subroutine flush_file(unit) - integer, intent(out) :: unit !< The I/O unit for the file to flush +!> Ensure that the output stream associated with a file handle is fully sent to disk. +subroutine flush_file_type(file) + type(file_type), intent(in) :: file !< The I/O handle for the file to flush + + call mpp_flush(file%unit) +end subroutine flush_file_type + +!> Ensure that the output stream associated with a unit is fully sent to disk. +subroutine flush_file_unit(unit) + integer, intent(in) :: unit !< The I/O unit for the file to flush call mpp_flush(unit) -end subroutine flush_file +end subroutine flush_file_unit !> Initialize the underlying I/O infrastructure subroutine io_infra_init(maxunits) @@ -179,9 +222,9 @@ subroutine write_version(version, tag, unit) end subroutine write_version !> open_file opens a file for parallel or single-file I/O. -subroutine open_file(unit, file, action, form, threading, fileset, nohdrs, domain, MOM_domain) +subroutine open_file_unit(unit, filename, action, form, threading, fileset, nohdrs, domain, MOM_domain) integer, intent(out) :: unit !< The I/O unit for the opened file - character(len=*), intent(in) :: file !< The name of the file being opened + character(len=*), intent(in) :: filename !< The name of the file being opened integer, optional, intent(in) :: action !< A flag indicating whether the file can be read !! or written to and how to handle existing files. integer, optional, intent(in) :: form !< A flag indicating the format of a new file. The @@ -198,13 +241,68 @@ subroutine open_file(unit, file, action, form, threading, fileset, nohdrs, domai type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition if (present(MOM_Domain)) then - call mpp_open(unit, file, action=action, form=form, threading=threading, fileset=fileset, & + call mpp_open(unit, filename, action=action, form=form, threading=threading, fileset=fileset, & nohdrs=nohdrs, domain=MOM_Domain%mpp_domain) else - call mpp_open(unit, file, action=action, form=form, threading=threading, fileset=fileset, & + call mpp_open(unit, filename, action=action, form=form, threading=threading, fileset=fileset, & nohdrs=nohdrs, domain=domain) endif -end subroutine open_file +end subroutine open_file_unit + +!> open_file opens a file for parallel or single-file I/O. +subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fileset) + type(file_type), intent(inout) :: IO_handle !< The handle for the opened file + character(len=*), intent(in) :: filename !< The path name of the file being opened + integer, optional, intent(in) :: action !< A flag indicating whether the file can be read + !! or written to and how to handle existing files. + !! The default is WRITE_ONLY. + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition + integer, optional, intent(in) :: threading !< A flag indicating whether one (SINGLE_FILE) + !! or multiple PEs (MULTIPLE) participate in I/O. + !! With the default, the root PE does I/O. + integer, optional, intent(in) :: fileset !< A flag indicating whether multiple PEs doing I/O due + !! to threading=MULTIPLE write to the same file (SINGLE_FILE) + !! or to one file per PE (MULTIPLE, the default). + + if (present(MOM_Domain)) then + call mpp_open(IO_handle%unit, filename, action=action, form=NETCDF_FILE, threading=threading, & + fileset=fileset, domain=MOM_Domain%mpp_domain) + else + call mpp_open(IO_handle%unit, filename, action=action, form=NETCDF_FILE, threading=threading, & + fileset=fileset) + endif + IO_handle%filename = trim(filename) + if (present(action)) then + if (action == READONLY_FILE) then + IO_handle%open_to_read = .true. ; IO_handle%open_to_write = .false. + else + IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .true. + endif + else + IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .true. + endif + +end subroutine open_file_type + +!> open_file opens an ascii file for parallel or single-file I/O using Fortran read and write calls. +subroutine open_ASCII_file(unit, file, action, threading, fileset) + integer, intent(out) :: unit !< The I/O unit for the opened file + character(len=*), intent(in) :: file !< The name of the file being opened + integer, optional, intent(in) :: action !< A flag indicating whether the file can be read + !! or written to and how to handle existing files. + integer, optional, intent(in) :: threading !< A flag indicating whether one (SINGLE_FILE) + !! or multiple PEs (MULTIPLE) participate in I/O. + !! With the default, the root PE does I/O. + integer, optional, intent(in) :: fileset !< A flag indicating whether multiple PEs doing I/O due + !! to threading=MULTIPLE write to the same file (SINGLE_FILE) + !! or to one file per PE (MULTIPLE, the default). + + call mpp_open(unit, file, action=action, form=ASCII_FILE, threading=threading, fileset=fileset, & + nohdrs=.true.) + +end subroutine open_ASCII_file + !> Provide a string to append to filenames, to differentiate ensemble members, for example. subroutine get_filename_suffix(suffix) @@ -216,8 +314,8 @@ end subroutine get_filename_suffix !> Get information about the number of dimensions, variables, global attributes and time levels !! in the file associated with an open file unit -subroutine get_file_info(unit, ndim, nvar, natt, ntime) - integer, intent(in) :: unit !< The I/O unit for the open file +subroutine get_file_info(IO_handle, ndim, nvar, natt, ntime) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O integer, optional, intent(out) :: ndim !< The number of dimensions in the file integer, optional, intent(out) :: nvar !< The number of variables in the file integer, optional, intent(out) :: natt !< The number of global attributes in the file @@ -226,7 +324,7 @@ subroutine get_file_info(unit, ndim, nvar, natt, ntime) ! Local variables integer :: ndims, nvars, natts, ntimes - call mpp_get_info( unit, ndims, nvars, natts, ntimes ) + call mpp_get_info(IO_handle%unit, ndims, nvars, natts, ntimes ) if (present(ndim)) ndim = ndims if (present(nvar)) nvar = nvars @@ -238,29 +336,29 @@ end subroutine get_file_info !> Get the times of records from a file !### Modify this to also convert to time_type, using information about the dimensions? -subroutine get_file_times(unit, time_values, ntime) - integer, intent(in) :: unit !< The I/O unit for the open file +subroutine get_file_times(IO_handle, time_values, ntime) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O real, allocatable, dimension(:), intent(inout) :: time_values !< The real times for the records in file. integer, optional, intent(out) :: ntime !< The number of time levels in the file integer :: ntimes if (allocated(time_values)) deallocate(time_values) - call get_file_info(unit, ntime=ntimes) + call get_file_info(IO_handle, ntime=ntimes) if (present(ntime)) ntime = ntimes if (ntimes > 0) then allocate(time_values(ntimes)) - call mpp_get_times(unit, time_values) + call mpp_get_times(IO_handle%unit, time_values) endif end subroutine get_file_times !> Set up the field information (e.g., names and metadata) for all of the variables in a file. The !! argument fields must be allocated with a size that matches the number of variables in a file. -subroutine get_file_fields(unit, fields) - integer, intent(in) :: unit !< The I/O unit for the open file +subroutine get_file_fields(IO_handle, fields) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O type(fieldtype), dimension(:), intent(inout) :: fields !< Field-type descriptions of all of !! the variables in a file. - call mpp_get_fields(unit, fields) + call mpp_get_fields(IO_handle%unit, fields) end subroutine get_file_fields !> Extract information from a field type, as stored or as found in a file @@ -378,15 +476,11 @@ subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. - integer :: is, ie, js, je - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & timelevel=timelevel, position=position) if (present(scale)) then ; if (scale /= 1.0) 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*data(is:ie,js:je) + call rescale_comp_data(MOM_Domain, data, scale) endif ; endif end subroutine MOM_read_data_2d @@ -420,8 +514,12 @@ subroutine MOM_read_data_2d_region(filename, fieldname, data, start, nread, MOM_ endif if (present(scale)) then ; if (scale /= 1.0) then - ! Dangerously rescale the whole array - data(:,:) = scale*data(:,:) + if (present(MOM_Domain)) then + call rescale_comp_data(MOM_Domain, data, scale) + else + ! Dangerously rescale the whole array + data(:,:) = scale*data(:,:) + endif endif ; endif end subroutine MOM_read_data_2d_region @@ -441,15 +539,11 @@ subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. - integer :: is, ie, js, je - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & timelevel=timelevel, position=position) if (present(scale)) then ; if (scale /= 1.0) 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*data(is:ie,js:je,:) + call rescale_comp_data(MOM_Domain, data, scale) endif ; endif end subroutine MOM_read_data_3d @@ -469,15 +563,11 @@ subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. - integer :: is, ie, js, je - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & timelevel=timelevel, position=position) if (present(scale)) then ; if (scale /= 1.0) 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*data(is:ie,js:je,:,:) + call rescale_comp_data(MOM_Domain, data, scale) endif ; endif end subroutine MOM_read_data_4d @@ -525,7 +615,6 @@ subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied !! by before they are returned. - integer :: is, ie, js, je integer :: u_pos, v_pos u_pos = EAST_FACE ; v_pos = NORTH_FACE @@ -541,12 +630,8 @@ subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data timelevel=timelevel, position=v_pos) 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) + call rescale_comp_data(MOM_Domain, u_data, scale) + call rescale_comp_data(MOM_Domain, v_data, scale) endif ; endif end subroutine MOM_read_vector_2d @@ -570,7 +655,6 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied !! by before they are returned. - integer :: is, ie, js, je integer :: u_pos, v_pos u_pos = EAST_FACE ; v_pos = NORTH_FACE @@ -586,20 +670,16 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data timelevel=timelevel, position=v_pos) 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,:) + call rescale_comp_data(MOM_Domain, u_data, scale) + call rescale_comp_data(MOM_Domain, v_data, scale) endif ; endif end subroutine MOM_read_vector_3d !> Write a 4d field to an output file. -subroutine write_field_4d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Field to write @@ -607,13 +687,13 @@ subroutine write_field_4d(io_unit, field_md, MOM_domain, field, tstamp, tile_cou integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & tile_count=tile_count, default_data=fill_value) end subroutine write_field_4d !> Write a 3d field to an output file. -subroutine write_field_3d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Field to write @@ -621,13 +701,13 @@ subroutine write_field_3d(io_unit, field_md, MOM_domain, field, tstamp, tile_cou integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & tile_count=tile_count, default_data=fill_value) end subroutine write_field_3d !> Write a 2d field to an output file. -subroutine write_field_2d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Field to write @@ -635,43 +715,43 @@ subroutine write_field_2d(io_unit, field_md, MOM_domain, field, tstamp, tile_cou integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & tile_count=tile_count, default_data=fill_value) end subroutine write_field_2d !> Write a 1d field to an output file. -subroutine write_field_1d(io_unit, field_md, field, tstamp) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine write_field_1d(IO_handle, field_md, field, tstamp) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp - call mpp_write(io_unit, field_md, field, tstamp=tstamp) + call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_1d !> Write a 0d field to an output file. -subroutine write_field_0d(io_unit, field_md, field, tstamp) - integer, intent(in) :: io_unit !< File I/O unit handle +subroutine write_field_0d(IO_handle, field_md, field, tstamp) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp - call mpp_write(io_unit, field_md, field, tstamp=tstamp) + call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_0d !> Write the data for an axis -subroutine MOM_write_axis(io_unit, axis) - integer, intent(in) :: io_unit !< File I/O unit handle - type(axistype), intent(in) :: axis !< An axis type variable with information to write +subroutine MOM_write_axis(IO_handle, axis) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(axistype), intent(in) :: axis !< An axis type variable with information to write - call mpp_write(io_unit, axis) + call mpp_write(IO_handle%unit, axis) end subroutine MOM_write_axis !> Store information about an axis in a previously defined axistype and write this !! information to the file indicated by unit. -subroutine write_metadata_axis( unit, axis, name, units, longname, cartesian, sense, domain, data, calendar) - integer, intent(in) :: unit !< The I/O unit for the file to write to +subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, data, calendar) + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(axistype), intent(inout) :: axis !< The axistype where this information is stored. character(len=*), intent(in) :: name !< The name in the file of this axis character(len=*), intent(in) :: units !< The units of this axis @@ -685,15 +765,15 @@ subroutine write_metadata_axis( unit, axis, name, units, longname, cartesian, se real, dimension(:), optional, intent(in) :: data !< The coordinate values of the points on this axis character(len=*), optional, intent(in) :: calendar !< The name of the calendar used with a time axis - call mpp_write_meta(unit, axis, name, units, longname, cartesian=cartesian, sense=sense, & + call mpp_write_meta(IO_handle%unit, axis, name, units, longname, cartesian=cartesian, sense=sense, & domain=domain, data=data, calendar=calendar) end subroutine write_metadata_axis !> Store information about an output variable in a previously defined fieldtype and write this !! information to the file indicated by unit. -subroutine write_metadata_field(unit, field, axes, name, units, longname, & +subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & min, max, fill, scale, add, pack, standard_name, checksum) - integer, intent(in) :: unit !< The I/O unit for the file to write to + type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(inout) :: field !< The fieldtype where this information is stored type(axistype), dimension(:), intent(in) :: axes !< Handles for the axis used for this variable character(len=*), intent(in) :: name !< The name in the file of this variable @@ -713,8 +793,8 @@ subroutine write_metadata_field(unit, field, axes, name, units, longname, & optional, intent(in) :: checksum !< Checksum values that can be used to verify reads. - call mpp_write_meta( unit, field, axes, name, units, longname, & - min=min, max=max, fill=fill, scale=scale, add=add, pack=pack, standard_name=standard_name, checksum=checksum) + call mpp_write_meta(IO_handle%unit, field, axes, name, units, longname, min=min, max=max, & + fill=fill, scale=scale, add=add, pack=pack, standard_name=standard_name, checksum=checksum) end subroutine write_metadata_field diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 29500875ca..0625177d77 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -8,11 +8,11 @@ module MOM_restart 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_grid, only : ocean_grid_type -use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file +use MOM_io, only : create_file, file_type, fieldtype, file_exists, open_file, close_file use MOM_io, only : MOM_read_data, read_data, MOM_write_field, read_field_chksum use MOM_io, only : get_file_info, get_file_fields, get_field_atts, get_file_times use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc, get_filename_appendix -use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE +use MOM_io, only : MULTIPLE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type, time_type_to_real, real_to_time @@ -874,7 +874,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ ! this should be 2 Gb or less. integer :: start_var, next_var ! The starting variables of the ! current and next files. - integer :: unit ! The I/O unit of the open file. + type(file_type) :: IO_handle ! The I/O handle of the open fileset integer :: m, nz, num_files, var_periods integer :: seconds, days, year, month, hour, minute character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. @@ -1020,33 +1020,31 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ enddo if (CS%parallel_restartfiles) then - call create_file(unit, trim(restartpath), vars, (next_var-start_var), & + call create_file(IO_handle, trim(restartpath), vars, (next_var-start_var), & fields, MULTIPLE, G=G, GV=GV, checksums=check_val) else - call create_file(unit, trim(restartpath), vars, (next_var-start_var), & + call create_file(IO_handle, trim(restartpath), vars, (next_var-start_var), & fields, SINGLE_FILE, G=G, GV=GV, checksums=check_val) endif do m=start_var,next_var-1 if (associated(CS%var_ptr3d(m)%p)) then - call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & CS%var_ptr3d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & CS%var_ptr2d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & CS%var_ptr4d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then - call MOM_write_field(unit, fields(m-start_var+1), CS%var_ptr1d(m)%p, & - restart_time) + call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr1d(m)%p, restart_time) elseif (associated(CS%var_ptr0d(m)%p)) then - call MOM_write_field(unit, fields(m-start_var+1), CS%var_ptr0d(m)%p, & - restart_time) + call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr0d(m)%p, restart_time) endif enddo - call close_file(unit) + call close_file(IO_handle) num_files = num_files+1 @@ -1086,7 +1084,7 @@ subroutine restore_state(filename, directory, day, G, CS) integer :: sizes(7) integer :: nvar, ntime, pos - integer :: unit(CS%max_fields) ! The I/O units of all open files. + type(file_type) :: IO_handles(CS%max_fields) ! The I/O units of all open files. character(len=200) :: unit_path(CS%max_fields) ! The file names. logical :: unit_is_global(CS%max_fields) ! True if the file is global. @@ -1104,10 +1102,10 @@ subroutine restore_state(filename, directory, day, G, CS) ! Get NetCDF ids for all of the restart files. if ((LEN_TRIM(filename) == 1) .and. (filename(1:1) == 'F')) then - num_file = open_restart_units('r', directory, G, CS, units=unit, & + num_file = open_restart_units('r', directory, G, CS, IO_handles=IO_handles, & file_paths=unit_path, global_files=unit_is_global) else - num_file = open_restart_units(filename, directory, G, CS, units=unit, & + num_file = open_restart_units(filename, directory, G, CS, IO_handles=IO_handles, & file_paths=unit_path, global_files=unit_is_global) endif @@ -1119,7 +1117,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! Get the time from the first file in the list that has one. do n=1,num_file - call get_file_times(unit(n), time_vals, ntime) + call get_file_times(IO_handles(n), time_vals, ntime) if (ntime < 1) cycle t1 = time_vals(1) @@ -1136,7 +1134,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! if they differ from the first time. if (is_root_pe()) then do m = n+1,num_file - call get_file_times(unit(n), time_vals, ntime) + call get_file_times(IO_handles(n), time_vals, ntime) if (ntime < 1) cycle t2 = time_vals(1) @@ -1153,10 +1151,10 @@ subroutine restore_state(filename, directory, day, G, CS) ! Read each variable from the first file in which it is found. do n=1,num_file - call get_file_info(unit(n), nvar=nvar) + call get_file_info(IO_handles(n), nvar=nvar) allocate(fields(nvar)) - call get_file_fields(unit(n), fields(1:nvar)) + call get_file_fields(IO_handles(n), fields(1:nvar)) do m=1, nvar call get_field_atts(fields(m), name=varname) @@ -1262,7 +1260,7 @@ subroutine restore_state(filename, directory, day, G, CS) enddo do n=1,num_file - call close_file(unit(n)) + call close_file(IO_handles(n)) enddo ! Check whether any mandatory fields have not been found. @@ -1355,7 +1353,7 @@ end function is_new_run !> open_restart_units determines the number of existing restart files and optionally opens !! them and returns unit ids, paths and whether the files are global or spatially decomposed. -function open_restart_units(filename, directory, G, CS, units, file_paths, & +function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, & global_files) result(num_files) character(len=*), intent(in) :: filename !< The list of restart file names or a single !! character 'r' to read automatically named files. @@ -1363,8 +1361,8 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & 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. - integer, dimension(:), & - optional, intent(out) :: units !< The I/O units of all opened files. + type(file_type), dimension(:), & + optional, intent(out) :: IO_handles !< The I/O handles of all opened files. character(len=*), dimension(:), & optional, intent(out) :: file_paths !< The full paths to open files. logical, dimension(:), & @@ -1444,22 +1442,22 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & num_restart = num_restart + 1 inquire(file=filepath, exist=fexists) if (fexists) then - if (present(units)) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & + if (present(IO_handles)) & + call open_file(IO_handles(n), trim(filepath), READONLY_FILE, & threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(n) = .true. elseif (CS%parallel_restartfiles) then ! Look for decomposed files using the I/O Layout. fexists = file_exists(filepath, G%Domain) - if (fexists .and. (present(units))) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, MOM_domain=G%Domain) + if (fexists .and. (present(IO_handles))) & + call open_file(IO_handles(n), trim(filepath), READONLY_FILE, MOM_domain=G%Domain) if (fexists .and. present(global_files)) global_files(n) = .false. endif if (fexists) then if (present(file_paths)) file_paths(n) = filepath n = n + 1 - if (is_root_pe() .and. (present(units))) & + if (is_root_pe() .and. (present(IO_handles))) & call MOM_error(NOTE, "MOM_restart: MOM run restarted using : "//trim(filepath)) else err = 1 ; exit @@ -1472,16 +1470,16 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & inquire(file=filepath, exist=fexists) if (fexists) then - if (present(units)) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & + if (present(IO_handles)) & + call open_file(IO_handles(n), trim(filepath), READONLY_FILE, & threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(n) = .true. if (present(file_paths)) file_paths(n) = filepath n = n + 1 - if (is_root_pe() .and. (present(units))) & + if (is_root_pe() .and. (present(IO_handles))) & call MOM_error(NOTE,"MOM_restart: MOM run restarted using : "//trim(filepath)) else - if (present(units)) & + if (present(IO_handles)) & call MOM_error(WARNING,"MOM_restart: Unable to find restart file : "//trim(filepath)) endif diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index 4ef00707fe..c9200cf41c 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -3,11 +3,11 @@ module MOM_write_cputime ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_coms, only : sum_across_PEs, num_pes +use MOM_coms, only : sum_across_PEs, num_pes use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe -use MOM_io, only : open_file, close_file, APPEND_FILE, ASCII_FILE, WRITEONLY_FILE -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_time_manager, only : time_type, get_time, operator(>) +use MOM_io, only : open_ASCII_file, close_file, APPEND_FILE, WRITEONLY_FILE +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_time_manager, only : time_type, get_time, operator(>) implicit none ; private @@ -181,11 +181,9 @@ subroutine write_cputime(day, n, CS, nmax, call_end) ! Reopen or create a text output file. if ((CS%previous_calls == 0) .and. (is_root_pe())) then if (day > CS%Start_time) then - call open_file(CS%fileCPU_ascii, trim(CS%CPUfile), & - action=APPEND_FILE, form=ASCII_FILE, nohdrs=.true.) + call open_ASCII_file(CS%fileCPU_ascii, trim(CS%CPUfile), action=APPEND_FILE) else - call open_file(CS%fileCPU_ascii, trim(CS%CPUfile), & - action=WRITEONLY_FILE, form=ASCII_FILE, nohdrs=.true.) + call open_ASCII_file(CS%fileCPU_ascii, trim(CS%CPUfile), action=WRITEONLY_FILE) endif endif diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index 454060414b..4f04fb285f 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -8,8 +8,9 @@ module MOM_coord_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type, log_version -use MOM_io, only : MOM_read_data, close_file, create_file, fieldtype, file_exists -use MOM_io, only : write_field, vardesc, var_desc, SINGLE_FILE, MULTIPLE +use MOM_io, only : close_file, create_file, file_type, fieldtype, file_exists +use MOM_io, only : MOM_read_data, MOM_write_field, vardesc, var_desc +use MOM_io, only : SINGLE_FILE, MULTIPLE use MOM_string_functions, only : slasher, uppercase use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs @@ -517,19 +518,19 @@ subroutine write_vertgrid_file(GV, US, param_file, directory) character(len=240) :: filepath type(vardesc) :: vars(2) type(fieldtype) :: fields(2) - integer :: unit + type(file_type) :: IO_handle ! The I/O handle of the fileset filepath = trim(directory) // trim("Vertical_coordinate") vars(1) = var_desc("R","kilogram meter-3","Target Potential Density",'1','L','1') vars(2) = var_desc("g","meter second-2","Reduced gravity",'1','L','1') - call create_file(unit, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) + call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) - call write_field(unit, fields(1), US%R_to_kg_m3*GV%Rlay(:)) - call write_field(unit, fields(2), US%L_T_to_m_s**2*US%m_to_Z*GV%g_prime(:)) + call MOM_write_field(IO_handle, fields(1), GV%Rlay, scale=US%R_to_kg_m3) + call MOM_write_field(IO_handle, fields(2), GV%g_prime, scale=US%L_T_to_m_s**2*US%m_to_Z) - call close_file(unit) + call close_file(IO_handle) end subroutine write_vertgrid_file diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 24c09a881c..70ef0768d5 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -11,9 +11,9 @@ module MOM_shared_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, param_file_type, log_version -use MOM_io, only : close_file, create_file, fieldtype, file_exists, field_size, stdout -use MOM_io, only : MOM_read_data, MOM_read_vector, read_variable, SINGLE_FILE, MULTIPLE -use MOM_io, only : open_file_to_read, close_file_to_read +use MOM_io, only : close_file, create_file, file_type, fieldtype, file_exists, field_size +use MOM_io, only : MOM_read_data, MOM_read_vector, read_variable, stdout +use MOM_io, only : open_file_to_read, close_file_to_read, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, MOM_write_field, var_desc use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type @@ -1189,42 +1189,30 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) integer, parameter :: nFlds=23 type(vardesc) :: vars(nFlds) type(fieldtype) :: fields(nFlds) - real :: Z_to_m_scale ! A unit conversion factor from Z to m. - real :: s_to_T_scale ! A unit conversion factor from T-1 to s-1. - real :: L_to_m_scale ! A unit conversion factor from L to m. - integer :: unit + real :: Z_to_m_scale ! A unit conversion factor from Z to m + real :: s_to_T_scale ! A unit conversion factor from T-1 to s-1 + real :: L_to_m_scale ! A unit conversion factor from L to m + type(file_type) :: IO_handle ! The I/O handle of the fileset integer :: file_threading integer :: nFlds_used - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB logical :: multiple_files - real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h - real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q - real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u - real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v call callTree_enter('write_ocean_geometry_file()') - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - Z_to_m_scale = 1.0 ; if (present(US)) Z_to_m_scale = US%Z_to_m s_to_T_scale = 1.0 ; if (present(US)) s_to_T_scale = US%s_to_T L_to_m_scale = 1.0 ; if (present(US)) L_to_m_scale = US%L_to_m -! vardesc is a structure defined in MOM_io.F90. The elements of -! this structure, in order, are: -! (1) the variable name for the NetCDF file -! (2) the variable's long name -! (3) a character indicating the horizontal grid, which may be '1' (column), -! 'h', 'q', 'u', or 'v', for the corresponding C-grid variable -! (4) a character indicating the vertical grid, which may be 'L' (layer), -! 'i' (interface), or '1' (no vertical location) -! (5) a character indicating the time levels of the field, which may be -! 's' (snap-shot), 'p' (periodic), or '1' (no time variation) -! (6) the variable's units + ! var_desc populates a type defined in MOM_io.F90. The arguments, in order, are: + ! (1) the variable name for the NetCDF file + ! (2) the units of the variable when output + ! (3) the variable's long name + ! (4) a character indicating the horizontal grid, which may be '1' (column), + ! 'h', 'q', 'u', or 'v', for the corresponding C-grid variable + ! (5) a character indicating the vertical grid, which may be 'L' (layer), + ! 'i' (interface), or '1' (no vertical location) + ! (6) a character indicating the time levels of the field, which may be + ! 's' (snap-shot), 'p' (periodic), or '1' (no time variation) vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1') vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1') vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1') @@ -1260,11 +1248,6 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) filepath = trim(directory) // "ocean_geometry" endif - out_h(:,:) = 0.0 - out_u(:,:) = 0.0 - out_v(:,:) = 0.0 - out_q(:,:) = 0.0 - call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, & "If true, each processor writes its own restart file, "//& "otherwise a single restart file is generated", & @@ -1272,67 +1255,40 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) file_threading = SINGLE_FILE if (multiple_files) file_threading = MULTIPLE - call create_file(unit, trim(filepath), vars, nFlds_used, fields, & - file_threading, dG=G) - - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = G%geoLatBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(1), G%Domain, out_q) - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = G%geoLonBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(2), G%Domain, out_q) - call MOM_write_field(unit, fields(3), G%Domain, G%geoLatT) - call MOM_write_field(unit, fields(4), G%Domain, G%geoLonT) - - do j=js,je ; do i=is,ie ; out_h(i,j) = Z_to_m_scale*G%bathyT(i,j) ; enddo ; enddo - call MOM_write_field(unit, fields(5), G%Domain, out_h) - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(i,J) = s_to_T_scale*G%CoriolisBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(6), G%Domain, out_q) - - ! I think that all of these copies are holdovers from a much earlier - ! ancestor code in which many of the metrics were macros that could have - ! had reduced dimensions, and that they are no longer needed in MOM6. -RWH - do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dxCv(i,J) ; enddo ; enddo - call MOM_write_field(unit, fields(7), G%Domain, out_v) - do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dyCu(I,j) ; enddo ; enddo - call MOM_write_field(unit, fields(8), G%Domain, out_u) - - do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dxCu(I,j) ; enddo ; enddo - call MOM_write_field(unit, fields(9), G%Domain, out_u) - do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dyCv(i,J) ; enddo ; enddo - call MOM_write_field(unit, fields(10), G%Domain, out_v) - - do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale*G%dxT(i,j); enddo ; enddo - call MOM_write_field(unit, fields(11), G%Domain, out_h) - do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale*G%dyT(i,j) ; enddo ; enddo - call MOM_write_field(unit, fields(12), G%Domain, out_h) - - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(i,J) = L_to_m_scale*G%dxBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(13), G%Domain, out_q) - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = L_to_m_scale*G%dyBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(14), G%Domain, out_q) - - do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale**2*G%areaT(i,j) ; enddo ; enddo - call MOM_write_field(unit, fields(15), G%Domain, out_h) - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = L_to_m_scale**2*G%areaBu(I,J) ; enddo ; enddo - call MOM_write_field(unit, fields(16), G%Domain, out_q) - - do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dx_Cv(i,J) ; enddo ; enddo - call MOM_write_field(unit, fields(17), G%Domain, out_v) - do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dy_Cu(I,j) ; enddo ; enddo - call MOM_write_field(unit, fields(18), G%Domain, out_u) - call MOM_write_field(unit, fields(19), G%Domain, G%mask2dT) + call create_file(IO_handle, trim(filepath), vars, nFlds_used, fields, file_threading, dG=G) + + call MOM_write_field(IO_handle, fields(1), G%Domain, G%geoLatBu) + call MOM_write_field(IO_handle, fields(2), G%Domain, G%geoLonBu) + call MOM_write_field(IO_handle, fields(3), G%Domain, G%geoLatT) + call MOM_write_field(IO_handle, fields(4), G%Domain, G%geoLonT) + + call MOM_write_field(IO_handle, fields(5), G%Domain, G%bathyT, scale=Z_to_m_scale) + call MOM_write_field(IO_handle, fields(6), G%Domain, G%CoriolisBu, scale=s_to_T_scale) + + call MOM_write_field(IO_handle, fields(7), G%Domain, G%dxCv, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(8), G%Domain, G%dyCu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(9), G%Domain, G%dxCu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(10), G%Domain, G%dyCv, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(11), G%Domain, G%dxT, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(12), G%Domain, G%dyT, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(13), G%Domain, G%dxBu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(14), G%Domain, G%dyBu, scale=L_to_m_scale) + + call MOM_write_field(IO_handle, fields(15), G%Domain, G%areaT, scale=L_to_m_scale**2) + call MOM_write_field(IO_handle, fields(16), G%Domain, G%areaBu, scale=L_to_m_scale**2) + + call MOM_write_field(IO_handle, fields(17), G%Domain, G%dx_Cv, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(18), G%Domain, G%dy_Cu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(19), G%Domain, G%mask2dT) if (G%bathymetry_at_vel) then - do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dblock_u(I,j) ; enddo ; enddo - call MOM_write_field(unit, fields(20), G%Domain, out_u) - do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dopen_u(I,j) ; enddo ; enddo - call MOM_write_field(unit, fields(21), G%Domain, out_u) - do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dblock_v(i,J) ; enddo ; enddo - call MOM_write_field(unit, fields(22), G%Domain, out_v) - do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dopen_v(i,J) ; enddo ; enddo - call MOM_write_field(unit, fields(23), G%Domain, out_v) + call MOM_write_field(IO_handle, fields(20), G%Domain, G%Dblock_u, scale=Z_to_m_scale) + call MOM_write_field(IO_handle, fields(21), G%Domain, G%Dopen_u, scale=Z_to_m_scale) + call MOM_write_field(IO_handle, fields(22), G%Domain, G%Dblock_v, scale=Z_to_m_scale) + call MOM_write_field(IO_handle, fields(23), G%Domain, G%Dopen_v, scale=Z_to_m_scale) endif - call close_file(unit) + call close_file(IO_handle) call callTree_leave('write_ocean_geometry_file()') end subroutine write_ocean_geometry_file diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 3e745dcafd..5050b6fce3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -190,8 +190,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call get_param(PF, mdl, "DEBUG", debug, default=.false.) call get_param(PF, mdl, "DEBUG_OBC", debug_obc, default=.false.) - new_sim = determine_is_new_run(dirs%input_filename, dirs%restart_input_dir, & - G, restart_CS) + new_sim = determine_is_new_run(dirs%input_filename, dirs%restart_input_dir, G, restart_CS) just_read = .not.new_sim call get_param(PF, mdl, "INPUTDIR", inputdir, & @@ -485,8 +484,7 @@ 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. - call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & - G, restart_CS) + call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, restart_CS) 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