io_restart.F90

MODULE io_RESTART
  use restart_file_group_module
  use restart_derivedtype_module
  use g_clock
  use g_config
  use o_arrays
  use g_backscatter
  use MOD_TRACER
  use MOD_ICE
  use MOD_DYN
  use MOD_MESH
  use MOD_PARTIT
  use MOD_PARSUP
  use fortran_utils
  use mpi
#if defined(__icepack)
  use icedrv_main
#endif 
#if defined (__cvmix)
  use g_cvmix_tke
  use g_cvmix_idemix
#endif  
#if defined(__recom)
  use recom_glovar
  use recom_config
  use recom_ciso
#endif
  
  implicit none
  public :: read_initial_conditions, write_initial_conditions, finalize_restart
  private

  integer,       save       :: globalstep=0 ! todo: remove this from module scope as it will mess things up if we use async read/write from the same process
  real(kind=WP)             :: ctime !current time in seconds from the beginning of the year

  type(restart_file_group) , save :: oce_files
  
  type(restart_file_group) , save :: ice_files
  
#if defined(__icepack)
  type(restart_file_group) , save, public :: icepack_files
#endif  

#if defined(__recom)
  type(restart_file_group) , save :: bio_files
#endif

  integer, parameter :: RAW_RESTART_METADATA_RANK = 0


  contains

!--------------------------------------------------------------------------------------------
! Helper functions for constructing restart file paths
!--------------------------------------------------------------------------------------------

! Build NetCDF restart file path
pure function nc_restart_path(component, year, root_path) result(path)
  implicit none
  character(len=*), intent(in) :: component, root_path
  integer, intent(in) :: year
  character(:), allocatable :: path
  character(4) :: cyear
  
  write(cyear, '(i4)') year
  path = trim(root_path) // trim(runid) // '.' // cyear // '.' // trim(component) // '.restart.nc'
end function nc_restart_path

! Build raw restart directory path
pure function build_raw_restart_dirpath(root_path) result(path)
  implicit none
  character(len=*), intent(in) :: root_path
  character(:), allocatable :: path
  
  path = trim(root_path) // trim(runid) // '_raw_restart'
end function build_raw_restart_dirpath

! Build raw restart info file path
pure function build_raw_restart_infopath(root_path) result(path)
  implicit none
  character(len=*), intent(in) :: root_path
  character(:), allocatable :: path
  
  path = trim(root_path) // trim(runid) // '_raw_restart'
end function build_raw_restart_infopath

! Build binary restart directory path
pure function build_bin_restart_dirpath(root_path) result(path)
  implicit none
  character(len=*), intent(in) :: root_path
  character(:), allocatable :: path
  
  path = trim(root_path) // trim(runid) // '_bin_restart'
end function build_bin_restart_dirpath

! Build binary restart info file path
pure function build_bin_restart_infopath(root_path) result(path)
  implicit none
  character(len=*), intent(in) :: root_path
  character(:), allocatable :: path
  
  path = trim(root_path) // trim(runid) // '_bin_restart'
end function build_bin_restart_infopath

!
!--------------------------------------------------------------------------------------------
! ini_ocean_io initializes ocean_file datatype which contains information of all variables need to be written into 
! the ocean restart file. This is the only place need to be modified if a new variable is added!
subroutine ini_ocean_io(dynamics, tracers, partit, mesh)
#ifdef ENABLE_NVHPC_WORKAROUNDS
  use nvfortran_subarray_workaround_module
#endif
  integer                   :: j, id
  character(500)            :: longname
  character(500)            :: trname, units
  type(t_mesh), target :: mesh
  type(t_partit), intent(inout), target :: partit
  type(t_tracer), target :: tracers
  type(t_dyn), target :: dynamics
  logical, save :: has_been_called = .false.

  if(has_been_called) return
  has_been_called = .true.

  !===========================================================================
  !===================== Definition part =====================================
  !===========================================================================
  !___Define the netCDF variables for 2D fields_______________________________
  !___SSH_____________________________________________________________________
  call oce_files%def_node_var('ssh', 'sea surface elevation', 'm',   dynamics%eta_n, mesh, partit)
  !___ALE related fields______________________________________________________
  call oce_files%def_node_var('hbar', 'ALE surface elevation', 'm',   mesh%hbar, mesh, partit)
!!PS   call oce_files%def_node_var('ssh_rhs', 'RHS for the elevation', '?',   ssh_rhs, mesh, partit)
  call oce_files%def_node_var('ssh_rhs_old', 'RHS for the elevation', '?',   dynamics%ssh_rhs_old, mesh, partit)
  call oce_files%def_node_var('hnode', 'nodal layer thickness', 'm',   mesh%hnode, mesh, partit)
  
  !___Define the netCDF variables for 3D fields_______________________________
#ifdef ENABLE_NVHPC_WORKAROUNDS
  dynamics_workaround => dynamics
#endif
  call oce_files%def_elem_var('u', 'zonal velocity',        'm/s', dynamics%uv(1,:,:), mesh, partit)
  call oce_files%def_elem_var('v', 'meridional velocity',   'm/s', dynamics%uv(2,:,:), mesh, partit)
  call oce_files%def_elem_var('urhs_AB', 'Adams-Bashforth for u (n-1 for AB2 and n-2 for AB3)', 'm/s', dynamics%uv_rhsAB(1,1,:,:), mesh, partit)
  call oce_files%def_elem_var('vrhs_AB', 'Adams-Bashforth for v (n-1 for AB2 and n-2 for AB3)', 'm/s', dynamics%uv_rhsAB(1,2,:,:), mesh, partit)
  if (dynamics%AB_order==3) then
       call oce_files%def_elem_var_optional('urhs_AB3', 'Adams-Bashforth for u (n-1) for AB3', 'm/s', dynamics%uv_rhsAB(2,1,:,:), mesh, partit)
       call oce_files%def_elem_var_optional('vrhs_AB3', 'Adams-Bashforth for v (n-1) for AB3', 'm/s', dynamics%uv_rhsAB(2,2,:,:), mesh, partit)
  end if
  
!___Save restart variables for TKE and IDEMIX_________________________________
!   if (trim(mix_scheme)=='cvmix_TKE' .or. trim(mix_scheme)=='cvmix_TKE+IDEMIX') then
#if defined (__cvmix)
  if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then
        call oce_files%def_node_var_optional('tke', 'Turbulent Kinetic Energy', 'm2/s2', tke(:,:), mesh, partit)
  endif
  if (mix_scheme_nmb==6 .or. mix_scheme_nmb==56) then
        call oce_files%def_elem_var_optional('iwe', 'Internal Wave Energy'    , 'm2/s2', iwe(:,:), mesh, partit)
  endif 
#endif  
  if (dynamics%opt_visc==8) then
        call oce_files%def_elem_var_optional('uke', 'unresolved kinetic energy', 'm2/s2', uke(:,:), mesh, partit)
        call oce_files%def_elem_var_optional('uke_rhs', 'unresolved kinetic energy rhs', 'm2/s2', uke_rhs(:,:), mesh, partit)
  endif
  
  do j=1,tracers%num_tracers
     id=tracers%data(j)%ID  !MB: Avoid hard-wired tracer assignments like SELECT CASE(j)
     SELECT CASE (id) 
       CASE(1)
         trname='temp'
         longname='potential temperature'
         units='degC'
       CASE(2)
         trname='salt'
         longname='salinity'
         units='psu'
       CASE(6)
         trname='sf6'
         longname='sulfur hexafluoride'
         units='mol / m**3'
       CASE(11)
         trname='cfc11'
         longname='chlorofluorocarbon CFC-11'
         units='mol / m**3'
       CASE(12)
         trname='cfc12'
         longname='chlorofluorocarbon CFC-12'
         units='mol / m**3'
       CASE(14)
         trname='r14c'
         longname='14C / C ratio of DIC'
         units='none'
       CASE(39)
         trname='r39ar'
         longname='39Ar / Ar ratio'
         units='none'
       CASE(101)
         trname='h2o18'
         longname='h2o18 concentration'
         units='kmol/m**3'
       CASE(102)
         trname='hDo16'
         longname='hDo16 concentration'
         units='kmol/m**3'
       CASE(103)
         trname='h2o16'
         longname='h2o16 concentration'
         units='kmol/m**3'
       CASE DEFAULT
         write(trname,'(A3,i4.4)') 'tra_', j
         write(longname,'(A15,i4.4)') 'passive tracer ', j
         units='none'
     END SELECT
     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103)) then
        call oce_files%def_node_var_optional(trim(trname), trim(longname), trim(units), tracers%data(j)%values(:,:), mesh, partit)
     else
        call oce_files%def_node_var(trim(trname), trim(longname), trim(units), tracers%data(j)%values(:,:), mesh, partit)
     endif
     longname=trim(longname)//', Adams-Bashforth'
     if ((tracers%data(j)%ID==101) .or. (tracers%data(j)%ID==102) .or. (tracers%data(j)%ID==103)) then
        call oce_files%def_node_var_optional(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:),    mesh, partit)
     else
        call oce_files%def_node_var(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:),    mesh, partit)
     endif
     call oce_files%def_node_var_optional(trim(trname)//'_M1', trim(longname), trim(units), tracers%data(j)%valuesold(1,:,:), mesh, partit)
     if (tracers%data(j)%AB_order==3) &
     call oce_files%def_node_var_optional(trim(trname)//'_M2', trim(longname), trim(units), tracers%data(j)%valuesold(2,:,:), mesh, partit)
  end do
  call oce_files%def_node_var('w', 'vertical velocity', 'm/s', dynamics%w, mesh, partit)
  call oce_files%def_node_var('w_expl', 'vertical velocity', 'm/s', dynamics%w_e, mesh, partit)
  call oce_files%def_node_var('w_impl', 'vertical velocity', 'm/s', dynamics%w_i, mesh, partit)
end subroutine ini_ocean_io
!
!--------------------------------------------------------------------------------------------
! ini_ice_io initializes ice_file datatype which contains information of all variables need to be written into
! the ice restart file. This is the only place need to be modified if a new variable is added!
subroutine ini_ice_io(ice, partit, mesh)
  type(t_mesh), intent(in) , target :: mesh
  type(t_partit), intent(inout), target :: partit
  type(t_ice), target :: ice
  logical, save :: has_been_called = .false.


  if(has_been_called) return
  has_been_called = .true.

  !===========================================================================
  !===================== Definition part =====================================
  !===========================================================================
  !___Define the netCDF variables for 2D fields_______________________________
  call ice_files%def_node_var('area', 'ice concentration [0 to 1]', '%',   ice%data(1)%values(:), mesh, partit)
  call ice_files%def_node_var('hice', 'effective ice thickness',    'm',   ice%data(2)%values(:), mesh, partit)
  call ice_files%def_node_var('hsnow', 'effective snow thickness',  'm',   ice%data(3)%values(:), mesh, partit)
  call ice_files%def_node_var('uice', 'zonal velocity',             'm/s', ice%uice, mesh, partit)
  call ice_files%def_node_var('vice', 'meridional velocity',        'm',   ice%vice, mesh, partit)
#if defined (__oifs) || defined (__ifsinterface)
  call ice_files%def_node_var_optional('ice_albedo', 'ice albedo',    '-',   ice%atmcoupl%ice_alb, mesh, partit)
  call ice_files%def_node_var_optional('ice_temp', 'ice surface temperature',  'K',   ice%data(4)%values, mesh, partit)
#endif /* (__oifs) */
#if defined (__oasis)
  !---wiso-code
  if (lwiso) then
    call ice_files%def_node_var_optional('h2o18_ice', 'h2o18 concentration in sea ice', 'kmol/m**3', tr_arr_ice(:,1), mesh, partit)
    call ice_files%def_node_var_optional('hDo16_ice', 'hDo16 concentration in sea ice', 'kmol/m**3', tr_arr_ice(:,2), mesh, partit)
    call ice_files%def_node_var_optional('h2o16_ice', 'h2o16 concentration in sea ice', 'kmol/m**3', tr_arr_ice(:,3), mesh, partit)
  end if
  !---wiso-code-end
#endif

end subroutine ini_ice_io
!
!--------------------------------------------------------------------------------------------
!

! ini_bio_io initializes bid datatype which contains information of all variables need to be written into
! the bio restart file. This is the only place need to be modified if a new variable is added!
#if defined(__recom)
subroutine ini_bio_io(tracers, partit, mesh)
    integer                   :: j
    character(500)            :: longname
    character(500)            :: trname, units

    type(t_mesh), intent(in) , target :: mesh
    type(t_partit), intent(inout), target :: partit
    type(t_tracer), target :: tracers
    logical, save :: has_been_called = .false.


    if(has_been_called) return
    has_been_called = .true.

    !===========================================================================
    !===================== Definition part =====================================
    !===========================================================================
    !___Define the netCDF variables for 2D fields_______________________________
    call bio_files%def_node_var('BenN',    'Benthos Nitrogen', 'mmol/m3',   Benthos(:,1), mesh, partit);
    call bio_files%def_node_var('BenC',    'Benthos Carbon',   'mmol/m3',   Benthos(:,2), mesh, partit);
    call bio_files%def_node_var('BenSi',   'Benthos Silicate', 'mmol/m3',   Benthos(:,3), mesh, partit);
    call bio_files%def_node_var('BenCalc', 'Benthos Calcite',  'mmol/m3',   Benthos(:,4), mesh, partit);
    call bio_files%def_node_var('HPlus',   'Conc. of H-plus ions in the surface water', 'mol/kg',   GloHplus, mesh, partit);

end subroutine ini_bio_io
#endif

!--------------------------------------------------------------------------------------------
! Separate subroutine for reading restart files (initial conditions)
subroutine read_initial_conditions(which_readr, ice, dynamics, tracers, partit, mesh)
  use fortran_utils
  implicit none
  
  ! Parameters
  type(t_mesh)  , intent(inout), target :: mesh
  type(t_partit), intent(inout), target :: partit
  type(t_tracer), intent(inout), target :: tracers
  type(t_dyn)   , intent(inout), target :: dynamics
  type(t_ice)   , intent(inout), target :: ice
  integer, intent(out) :: which_readr
  
  ! Local variables
  logical :: rawfiles_exist, binfiles_exist
  integer :: mpierr
  character(:), allocatable :: read_raw_dirpath, read_raw_infopath
  character(:), allocatable :: read_bin_dirpath, read_bin_infopath
  character(:), allocatable :: read_oce_path, read_ice_path, read_bio_path
  
  ! Build paths for reading using RestartInPath
  read_raw_dirpath = build_raw_restart_dirpath(RestartInPath)//"/np"//int_to_txt(partit%npes)
  read_raw_infopath = build_raw_restart_infopath(RestartInPath)//"/np"//int_to_txt(partit%npes)//".info"
  read_bin_dirpath = build_bin_restart_dirpath(RestartInPath)//"/np"//int_to_txt(partit%npes)
  read_bin_infopath = build_bin_restart_infopath(RestartInPath)//"/np"//int_to_txt(partit%npes)//".info"
  read_oce_path = nc_restart_path('oce', yearold, RestartInPath)
  read_ice_path = nc_restart_path('ice', yearold, RestartInPath)
  read_bio_path = nc_restart_path('bio', yearold, RestartInPath)

  ! Initialize file groups for reading
  call ini_ocean_io(dynamics, tracers, partit, mesh)
  if (use_ice) then
#if defined(__icepack)    
      call ini_icepack_io(yearold, partit, mesh)
#else
      call ini_ice_io  (ice, partit, mesh)
#endif        
  end if     
#if defined(__recom)
  if (REcoM_restart) call ini_bio_io(tracers, partit, mesh)
#endif
  
  ! Check for raw restart files
  if(partit%mype == RAW_RESTART_METADATA_RANK) then
    inquire(file=read_raw_infopath, exist=rawfiles_exist)
  end if
  call MPI_Bcast(rawfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, mpierr)
  
  ! Check for binary restart files
  if(partit%mype == RAW_RESTART_METADATA_RANK) then
    inquire(file=read_bin_infopath, exist=binfiles_exist)
  end if
  call MPI_Bcast(binfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, mpierr)

  ! Read restart files in order of preference
  if(rawfiles_exist) then
    ! Read raw/core dump restart
    which_readr = 1
    ! Note: This will need to be updated once we have read functions that accept paths
    call read_all_raw_restarts(read_raw_dirpath, read_raw_infopath, partit%MPI_COMM_FESOM, partit%mype)
    
  elseif(binfiles_exist .and. bin_restart_length_unit /= "off") then
    ! Read binary restart
    which_readr = 2
    if (use_ice) then 
        call read_all_bin_restarts(read_bin_dirpath, &
                                   partit   = partit,   &
                                   mesh     = mesh,     &
                                   ice      = ice,      &
                                   dynamics = dynamics, &
                                   tracers  = tracers   )
    else
        call read_all_bin_restarts(read_bin_dirpath, &
                                   partit   = partit,   &
                                   mesh     = mesh,     &                    
                                   dynamics = dynamics, &
                                   tracers  = tracers   )
    end if     
    
  else
    ! Read NetCDF restart
    which_readr = 0
    
    ! Read OCEAN restart
    if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ocean'//achar(27)//'[0m'
    call read_netcdf_restarts(read_oce_path, oce_files, partit%MPI_COMM_FESOM, partit%mype)
    
    ! Read ICE/ICEPACK restart
    if (use_ice) then
#if defined(__icepack)   
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: icepack'//achar(27)//'[0m'
        call read_netcdf_restarts(nc_restart_path('icepack', yearold, RestartInPath), icepack_files, partit%MPI_COMM_FESOM, partit%mype)
#else            
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ice'//achar(27)//'[0m'
        call read_netcdf_restarts(read_ice_path, ice_files, partit%MPI_COMM_FESOM, partit%mype)            
#endif
    end if 

#if defined(__recom)
    ! Read RECOM restarts
    if (REcoM_restart) then
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: bio'//achar(27)//'[0m'
        call read_netcdf_restarts(read_bio_path, bio_files, partit%MPI_COMM_FESOM, partit%mype)
    end if
#endif

  end if
  
end subroutine read_initial_conditions

!--------------------------------------------------------------------------------------------
! Separate subroutine for writing restart files
subroutine write_initial_conditions(istep, nstart, ntotal, which_readr, ice, dynamics, tracers, partit, mesh)
  use fortran_utils
  implicit none
  
  ! Parameters
  integer :: istep, nstart, ntotal
  type(t_mesh)  , intent(inout), target :: mesh
  type(t_partit), intent(inout), target :: partit
  type(t_tracer), intent(inout), target :: tracers
  type(t_dyn)   , intent(inout), target :: dynamics
  type(t_ice)   , intent(inout), target :: ice
  integer, intent(in) :: which_readr
  
  ! Local variables
  logical :: is_portable_restart_write, is_raw_restart_write, is_bin_restart_write
  logical, save :: initialized_raw = .false.
  logical, save :: initialized_bin = .false.
  logical, save :: initialized_io = .false.
  integer :: mpierr

  character(:), allocatable :: write_raw_dirpath, write_raw_infopath
  character(:), allocatable :: write_bin_dirpath, write_bin_infopath
  character(:), allocatable :: write_oce_path, write_ice_path
  character(:), allocatable :: write_icepack_path, write_bio_path


  ! Build paths for reading using RestartInPath
  write_raw_dirpath = build_raw_restart_dirpath(RestartOutPath)//"/np"//int_to_txt(partit%npes)
  write_raw_infopath = build_raw_restart_infopath(RestartOutPath)//"/np"//int_to_txt(partit%npes)//".info"
  write_bin_dirpath = build_bin_restart_dirpath(RestartOutPath)//"/np"//int_to_txt(partit%npes)
  write_bin_infopath = build_bin_restart_infopath(RestartOutPath)//"/np"//int_to_txt(partit%npes)//".info"
  write_oce_path = nc_restart_path('oce', yearnew, RestartOutPath)
  write_ice_path = nc_restart_path('ice', yearnew, RestartOutPath)
  write_icepack_path = nc_restart_path('icepack', yearnew, RestartOutPath)
  write_bio_path = nc_restart_path('bio', yearnew, RestartOutPath)
  
  !_____________________________________________________________________________
  ! Initialize output directories on first call
  if(.not. initialized_raw) then
    initialized_raw = .true.
    if(raw_restart_length_unit /= "off") then
      if(partit%mype == RAW_RESTART_METADATA_RANK) then
        call mkdir(build_raw_restart_dirpath(RestartOutPath))
        call mkdir(write_raw_dirpath)
      end if
      call MPI_Barrier(partit%MPI_COMM_FESOM, mpierr)
    end if
  end if

  if(.not. initialized_bin) then
    initialized_bin = .true.
    if(bin_restart_length_unit /= "off") then
        if(partit%mype == RAW_RESTART_METADATA_RANK) then
            call mkdir(build_bin_restart_dirpath(RestartOutPath))
            call mkdir(write_bin_dirpath)
        end if
        call MPI_Barrier(partit%MPI_COMM_FESOM, mpierr)
    end if
  end if
  
  ! Initialize file groups for writing on first call
  if(.not. initialized_io) then
    initialized_io = .true.
    call ini_ocean_io(dynamics, tracers, partit, mesh)
    if (use_ice) then
#if defined(__icepack)
        call ini_icepack_io(yearnew, partit, mesh)
#else        
        call ini_ice_io(ice, partit, mesh)        
#endif        
    end if     
#if defined(__recom)
    if (use_REcoM) call ini_bio_io(tracers, partit, mesh)
#endif
  end if
  
  ! Skip writing on step 0
  if (istep==0) return
  
  ! Calculate current time from clock (seconds from beginning of year)
  ctime = timeold + (dayold - 1.0_WP) * 86400.0_WP
  
  ! Check whether restart will be written
  is_portable_restart_write = is_due(trim(restart_length_unit), restart_length, istep)
  
  ! Should write core dump restart?
  if(is_portable_restart_write .and. (raw_restart_length_unit /= "off")) then
    is_raw_restart_write = .true. ! always write a raw restart together with the portable restart
  else
#if !defined __ifsinterface
    is_raw_restart_write = is_due(trim(raw_restart_length_unit), raw_restart_length, istep)
#else
    is_raw_restart_write = is_due(trim(raw_restart_length_unit), raw_restart_length, istep) .OR. (istep==ntotal)
#endif
  end if
  
  ! Should write derived type binary restart?
  if(is_portable_restart_write .and. (bin_restart_length_unit /= "off")) then
    is_bin_restart_write = .true. ! always write a binary restart together with the portable restart
  else
    is_bin_restart_write = is_due(trim(bin_restart_length_unit), bin_restart_length, istep)
  end if

  ! Write restart files
  if(is_portable_restart_write) then
    ! Write OCEAN restart
    if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ocean'//achar(27)//'[0m'
    call write_netcdf_restarts(write_oce_path, oce_files, istep)
    
    ! Write ICE/ICEPACK restart
    if(use_ice) then
#if defined(__icepack)        
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: icepack'//achar(27)//'[0m'
        call write_netcdf_restarts(write_icepack_path, icepack_files, istep)
#else
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ice'//achar(27)//'[0m'
        call write_netcdf_restarts(write_ice_path, ice_files, istep)
#endif 
    end if

#if defined(__recom)
    ! Write RECOM restart
    if (REcoM_restart .or. use_REcoM) then
        if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: bio'//achar(27)//'[0m'
        call write_netcdf_restarts(write_bio_path, bio_files, istep)
    end if
#endif
  end if

  ! Write core dump
  if(is_raw_restart_write) then
    call write_all_raw_restarts(write_raw_dirpath, write_raw_infopath, istep, partit%MPI_COMM_FESOM, partit%mype)
  end if

  ! Write derived type binary
  if(is_bin_restart_write) then
    call write_all_bin_restarts((/globalstep+istep, int(ctime), yearnew/), &
                                write_bin_dirpath,                   &
                                write_bin_infopath,                  &
                                partit,                              &
                                mesh,                                &                                
                                ice,                                 &
                                dynamics,                            &
                                tracers                              )
  end if

  ! Update clock file to latest restart point
  if (partit%mype==0) then
    if(is_portable_restart_write .or. is_raw_restart_write .or. is_bin_restart_write) then
        write(*,*) ' --> actualize clock file to latest restart point'
        call clock_finish
    end if
  end if

end subroutine write_initial_conditions
!
!
!_______________________________________________________________________________
subroutine write_netcdf_restarts(path, filegroup, istep)
  use fortran_utils
  character(len=*), intent(in) :: path
  type(restart_file_group), intent(inout) :: filegroup
  integer,  intent(in)          :: istep
  ! EO parameters
  integer cstep
  integer i
  character(:), allocatable :: dirpath
  character(:), allocatable :: filepath
  logical file_exists
  
  cstep = globalstep+istep
  
  ! Calculate current time from clock (seconds from beginning of year)
  ctime = timeold + (dayold - 1.0_WP) * 86400.0_WP
  
  do i=1, filegroup%nfiles
    call filegroup%files(i)%join() ! join the previous write (if required)

    if(filegroup%files(i)%is_iorank()) then
      if(filegroup%files(i)%is_attached()) call filegroup%files(i)%close_file() ! close the file from previous write
            
      dirpath = path(1:len(path)-3) ! chop of the ".nc" suffix
      filepath = dirpath//"/"//filegroup%files(i)%varname//".nc"

      if(filegroup%files(i)%path == "" .or. (.not. filegroup%files(i)%must_exist_on_read)) then
        ! the path to an existing restart file is not set in read_netcdf_restarts if we had a restart from a raw restart
        ! OR we might have skipped the file when reading restarts and it does not exist at all
        inquire(file=filepath, exist=file_exists)
        if(file_exists) then
          filegroup%files(i)%path = filepath
        else if(.not. filegroup%files(i)%must_exist_on_read) then
          filegroup%files(i)%path = ""
        end if
      end if
      if(filegroup%files(i)%path .ne. filepath) then
        ! execute_command_line with mkdir sometimes fails, use a custom implementation around mkdir from C instead
        call mkdir(dirpath)
        filegroup%files(i)%path = filepath
        call filegroup%files(i)%open_write_create(filegroup%files(i)%path)
      else
        call filegroup%files(i)%open_write_append(filegroup%files(i)%path) ! todo: keep the file open between writes
      end if

      write(*,*) 'writing restart record ', filegroup%files(i)%rec_count()+1, ' to ', filegroup%files(i)%path
      call filegroup%files(i)%write_var(filegroup%files(i)%iter_varindex, [filegroup%files(i)%rec_count()+1], [1], [cstep])
      ! todo: write time via the fesom_file_type
      call filegroup%files(i)%write_var(filegroup%files(i)%time_varindex(), [filegroup%files(i)%rec_count()+1], [1], [ctime])
    end if

    call filegroup%files(i)%async_gather_and_write_variables()
  end do
  
end subroutine write_netcdf_restarts
!
!
!_______________________________________________________________________________
subroutine write_all_raw_restarts(dirpath, infopath, istep, mpicomm, mype)
  character(len=*), intent(in) :: dirpath
  character(len=*), intent(in) :: infopath
  integer,  intent(in):: istep
  integer, intent(in) :: mpicomm
  integer, intent(in) :: mype
  ! EO parameters
  integer cstep
  integer fileunit

  open(newunit = fileunit, file = dirpath//'/'//mpirank_to_txt(mpicomm)//'.dump', form = 'unformatted')
  call write_raw_restart_group(oce_files, fileunit)
  if(use_ice) call write_raw_restart_group(ice_files, fileunit)
#if defined(__recom)
  call write_raw_restart_group(bio_files, fileunit)
#endif
  close(fileunit)

  if(mype == RAW_RESTART_METADATA_RANK) then
    print *,"writing raw restart to "//dirpath
    ! store metadata about the raw restart
    cstep = globalstep+istep
    open(newunit = fileunit, file = infopath)
    write(fileunit, '(g0)') cstep
    write(fileunit, '(g0)') ctime
    write(fileunit, '(2(g0))') "! year: ",yearnew
    write(fileunit, '(3(g0))') "! oce: ", oce_files%nfiles, " variables"
    if(use_ice) write(fileunit, '(3(g0))') "! ice: ", ice_files%nfiles, " variables"
#if defined(__recom)
    write(fileunit, '(3(g0))') "! bio: ", bio_files%nfiles, " variables"
#endif
    close(fileunit)
  end if
end subroutine write_all_raw_restarts
!
!
!_______________________________________________________________________________
subroutine write_raw_restart_group(filegroup, fileunit)
  type(restart_file_group), intent(inout) :: filegroup
  integer, intent(in) :: fileunit
  ! EO parameters
  integer i
  
  do i=1, filegroup%nfiles
    call filegroup%files(i)%write_variables_raw(fileunit)
  end do
end subroutine write_raw_restart_group

subroutine read_all_raw_restarts(dirpath, infopath, mpicomm, mype)
  character(len=*), intent(in) :: dirpath
  character(len=*), intent(in) :: infopath
  integer, intent(in) :: mpicomm
  integer, intent(in) :: mype
  ! EO parameters
  integer rstep
  real(kind=WP) rtime
  integer fileunit
  integer status
  integer mpierr

  if(mype == RAW_RESTART_METADATA_RANK) then
    ! read metadata info for the raw restart
    print *,"trying to open ", infopath
    open(newunit = fileunit, status = 'old', iostat = status, file = infopath)
    if(status == 0) then
      read(fileunit,*) rstep
      read(fileunit,*) rtime
      close(fileunit)
    else
      print *,"can not open ", infopath
      stop 1
    end if
    
    
    ! compare the restart time with our actual time
    if(int(ctime) /= int(rtime)) then
        write(*,*)
        print *, achar(27)//'[5,33m'  
        write(*,*) '____________________________________________________________________'
        write(*,*) "WARNING: raw restart time ",rtime," does not match current clock time",ctime
        write(*,*) "         If you restart with a different time step this might be ok!"
        write(*,*) "         If that is not the case, check your fesom.clock file and"
        write(*,*) "         the time information of your restart file, make sure they "
        write(*,*) "         are not messed up!!!"
        write(*,*) '____________________________________________________________________'
        print *, achar(27)//'[0m'
        write(*,*)
    end if
    globalstep = rstep
    print *,"reading raw restart from " // dirpath
  end if
  ! sync globalstep with the other processes to let all processes writing portable restart files know the globalstep
  call MPI_Bcast(globalstep, 1, MPI_INTEGER, RAW_RESTART_METADATA_RANK, mpicomm, mpierr)

  open(newunit = fileunit, status = 'old', iostat = status, file = dirpath//'/'//mpirank_to_txt(mpicomm)//'.dump', form = 'unformatted')
  if(status == 0) then
    call read_raw_restart_group(oce_files, fileunit)
    if(use_ice) call read_raw_restart_group(ice_files, fileunit)
#if defined(__recom)
    call read_raw_restart_group(bio_files, fileunit)
#endif
    close(fileunit)
  else
    print *,"can not open ",dirpath//'/'//mpirank_to_txt(mpicomm)//'.dump'
    stop 1
  end if
end subroutine read_all_raw_restarts
!
!
!_______________________________________________________________________________
subroutine read_raw_restart_group(filegroup, fileunit)
  type(restart_file_group), intent(inout) :: filegroup
  integer, intent(in) :: fileunit
  ! EO parameters
  integer i
  
  do i=1, filegroup%nfiles
    call filegroup%files(i)%read_variables_raw(fileunit)
  end do  
end subroutine read_raw_restart_group
!
!
!_______________________________________________________________________________
! join remaining threads and close all open files
subroutine finalize_restart()
  integer i

  ! join all previous writes
  ! close all restart files

  do i=1, oce_files%nfiles
    call oce_files%files(i)%join()
    if(oce_files%files(i)%is_iorank()) then
      if(oce_files%files(i)%is_attached()) call oce_files%files(i)%close_file()
    end if
  end do

  if(use_ice) then
    do i=1, ice_files%nfiles
      call ice_files%files(i)%join()
      if(ice_files%files(i)%is_iorank()) then
        if(ice_files%files(i)%is_attached()) call ice_files%files(i)%close_file()
      end if
    end do
  end if
#if defined(__recom)
  do i=1, bio_files%nfiles
    call bio_files%files(i)%join()
    if(bio_files%files(i)%is_iorank()) then
      if(bio_files%files(i)%is_attached()) call bio_files%files(i)%close_file()
    end if
  end do
#endif
end subroutine finalize_restart
!
!
!_______________________________________________________________________________
subroutine read_netcdf_restarts(path, filegroup, mpicomm, mype)
  character(len=*), intent(in) :: path
  type(restart_file_group), intent(inout) :: filegroup
  integer, intent(in) :: mpicomm
  integer, intent(in) :: mype
  ! EO parameters
  real(kind=WP) rtime
  integer i
  character(:), allocatable :: dirpath
  integer mpistatus(MPI_STATUS_SIZE)
  logical file_exists
  logical, allocatable :: skip_file(:)
  integer current_iorank_snd, current_iorank_rcv
  integer max_globalstep
  integer mpierr
  
  ! Calculate current time from clock (seconds from beginning of year)
  ctime = timeold + (dayold - 1.0_WP) * 86400.0_WP
  
  allocate(skip_file(filegroup%nfiles))
  skip_file = .false.
  
  do i=1, filegroup%nfiles
    current_iorank_snd = 0
    current_iorank_rcv = 0
    if( filegroup%files(i)%is_iorank() ) then
      dirpath = path(1:len(path)-3) ! chop of the ".nc" suffix
      if(filegroup%files(i)%path .ne. dirpath//"/"//filegroup%files(i)%varname//".nc") then
        filegroup%files(i)%path = dirpath//"/"//filegroup%files(i)%varname//".nc"

        ! determine if the file should be skipped
        if(.not. filegroup%files(i)%must_exist_on_read) then
          current_iorank_snd = mype
          inquire(file=filegroup%files(i)%path, exist=file_exists)
          if(.not. file_exists) skip_file(i) = .true.
        end if

        if(.not. skip_file(i)) then
          write(*,*) 'reading restart for ', filegroup%files(i)%varname, ' at ', filegroup%files(i)%path
        else
          write(*,*) 'skipping reading restart for ', filegroup%files(i)%varname, ' at ', filegroup%files(i)%path
        end if
        
        if(.not. skip_file(i)) call filegroup%files(i)%open_read(filegroup%files(i)%path) ! do we need to bother with read-only access?
        ! todo: print a reasonable error message if the file does not exist
      end if      
    end if

    ! iorank already knows if we skip the file, tell the others
    if(.not. filegroup%files(i)%must_exist_on_read) then
      call MPI_Allreduce(current_iorank_snd, current_iorank_rcv, 1, MPI_INTEGER, MPI_SUM, mpicomm, mpierr)
      call MPI_Bcast(skip_file(i), 1, MPI_LOGICAL, current_iorank_rcv, mpicomm, mpierr)
    end if      

    ! ========================================================================!
    !                           _____________________                         !
    !                          /                     \                        !
    !                         /                       \                       !
    !                        |      .-----------.      |                      !
    !                        |     /             \     |                      !
    !                        |    |               |    |                      !
    !                        |    |    R.I.P.     |    |                      !
    !                        |    |               |    |                      !
    !                        |    | MULTITHREADED |    |                      !
    !                        |    |    RESTART    |    |                      !
    !                        |    |    READING    |    |                      !
    !                        |    |               |    |                      !
    !                        |    |  2021 - 2025  |    |                      !
    !                        |    |               |    |                      !
    !                        |    | "You had to   |    |                      !
    !                        |    |  be worked    |    |                      !
    !                        |    |  around on    |    |                      !
    !                        |    |  every        |    |                      !
    !                        |    |  machine"     |    |                      !
    !                        |    |               |    |                      !
    !                        |     \             /     |                      !
    !                        |      '-----------'      |                      !
    !                        |                         |                      !
    !                        |_________________________|                      !
    !                                                                         !
    !  Cause of death: MPI threading issues, deadlocks, and performance       !
    !                  degradation on Levante, Albedo, Aleph, Juwels          !
    !                                                                         !
    !  Survived by: Sequential file reading (simple, reliable, boring)        !
    !                                                                         !
    !  For more info and a possible resurrection see:                         !
    !  https://github.com/FESOM/fesom2/pull/801                               !
    ! ========================================================================!

    if(.not. skip_file(i)) then
      call filegroup%files(i)%read_and_scatter_variables()
    end if

    if(skip_file(i)) cycle

    if(filegroup%files(i)%is_iorank()) then
      write(*,*) 'restart from record ', filegroup%files(i)%rec_count(), ' of ', filegroup%files(i)%rec_count(), filegroup%files(i)%path

      ! read the last entry from the iter variable
      call filegroup%files(i)%read_var1(filegroup%files(i)%iter_varindex, [filegroup%files(i)%rec_count()], globalstep)

      ! read the last entry from the time variable
      call filegroup%files(i)%read_var1(filegroup%files(i)%time_varindex(), [filegroup%files(i)%rec_count()], rtime)
      call filegroup%files(i)%close_file()

     if (int(ctime)/=int(rtime)) then
        print *, achar(27)//'[33m'    //'____________________________________________________________'//achar(27)//'[0m'
        print *, achar(27)//'[5;33m'  //' --> WARNING: RESTART TIMESTAMP MISMATCH !!!               '//achar(27)//'[0m'
        write(*,*) 'WARNING: timestamps in restart and in clock files do not match for ', filegroup%files(i)%varname
        write(*,*) '         at path: ', trim(filegroup%files(i)%path)
        write(*,*) '         clock time  =', ctime
        write(*,*) '         restart time=', rtime
        write(*,*) 'WARNING: This mismatch will be ignored and the model will continue running.'
        write(*,*) 'WARNING: Please verify that this is the intended behavior for your simulation.'
        print *, achar(27)//'[33m'    //'____________________________________________________________'//achar(27)//'[0m'
      end if
    end if
  end do

  ! sync globalstep with processes which may have skipped a restart upon reading and thus need to know the globalstep when writing their restart 
  if( any(skip_file .eqv. .true.) ) then
    call MPI_Allreduce(globalstep, max_globalstep, 1, MPI_INTEGER, MPI_MAX, mpicomm, mpierr)
    globalstep = max_globalstep
  end if

  ! sync globalstep with the process responsible for raw restart metadata
  if(filegroup%nfiles >= 1) then
    ! use the first restart I/O process to send the globalstep
    if( filegroup%files(1)%is_iorank() .and. (mype .ne. RAW_RESTART_METADATA_RANK)) then
      call MPI_Send(globalstep, 1, MPI_INTEGER, RAW_RESTART_METADATA_RANK, 42, mpicomm, mpierr)
    else if((mype == RAW_RESTART_METADATA_RANK) .and. (.not. filegroup%files(1)%is_iorank())) then
      call MPI_Recv(globalstep, 1, MPI_INTEGER, MPI_ANY_SOURCE, 42, mpicomm, mpistatus, mpierr)
    end if
  end if
end subroutine read_netcdf_restarts
!
!
!_______________________________________________________________________________
  function is_due(unit, frequency, istep) result(d)
    character(len=*), intent(in) :: unit
    integer, intent(in) :: frequency
    integer, intent(in) :: istep
    logical d
    ! EO parameters
    d = .false.
    
    if(unit.eq.'y') then
      call annual_event(d, frequency)
    else if(unit.eq.'m') then 
      call monthly_event(d, frequency) 
    else if(unit.eq.'d') then
      call daily_event(d, frequency)
    else if(unit.eq.'h') then
      call hourly_event(d, frequency)
    else if(unit.eq.'s') then
      call step_event(d, istep, frequency)
    else if(unit.eq.'off') then
      d = .false.
    else
      write(*,*) 'You did not specify a supported outputflag.'
      write(*,*) 'The program will stop to give you opportunity to do it.'
      stop 1
    stop
    end if
  end function
! !
! !
! !_______________________________________________________________________________
!   function mpirank_to_txt(mpicomm) result(txt)
!     use fortran_utils
!     integer, intent(in) :: mpicomm
!     character(:), allocatable :: txt
!     ! EO parameters
!     integer mype
!     integer npes
!     integer mpierr
!   
!     call MPI_Comm_Rank(mpicomm, mype, mpierr)
!     call MPI_Comm_Size(mpicomm, npes, mpierr)
!     txt = int_to_txt_pad(mype,int(log10(real(npes)))+1) ! pad to the width of the number of processes
!   end function
!!PS --> move this function also to fortran_utils.F90 

end module