From aff740d72ca3576a352f94c8d4f2a0a275f5bffc Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 7 Oct 2024 16:30:21 +0200 Subject: [PATCH 01/18] try to make icepack to compile. Try to preliminary fix the icepack restart. This still needs to be checked --- src/icepack_drivers/icedrv_io.F90 | 100 +++++++++++++--------------- src/icepack_drivers/icedrv_main.F90 | 17 +++-- src/icepack_drivers/icedrv_step.F90 | 35 ++++++---- src/io_fesom_file.F90 | 39 +++++++++-- src/io_meandata.F90 | 4 +- src/io_restart.F90 | 45 ++++++++++--- src/io_restart_file_group.F90 | 55 +++++++++++++-- 7 files changed, 199 insertions(+), 96 deletions(-) diff --git a/src/icepack_drivers/icedrv_io.F90 b/src/icepack_drivers/icedrv_io.F90 index 7d8bb7c4d..209d15cb7 100644 --- a/src/icepack_drivers/icedrv_io.F90 +++ b/src/icepack_drivers/icedrv_io.F90 @@ -18,7 +18,11 @@ contains - module subroutine init_io_icepack(mesh) + ! + ! + !_________________________________________________________________________ + ! define mean IO output of icepack + module subroutine ini_mean_icepack_io(mesh) use mod_mesh use io_meandata, only: def_stream3D, def_stream2D @@ -269,32 +273,33 @@ module subroutine init_io_icepack(mesh) end select end do - end subroutine init_io_icepack + end subroutine ini_mean_icepack_io ! - !-------------------------------------------------------------------------------------------- ! - - module subroutine init_restart_icepack(year, mesh) + !___________________________________________________________________________ + ! define mean IO output of icepack + module subroutine ini_icepack_io(year, partit, mesh) use mod_mesh + use mod_partit + use mod_parsup use g_config, only: runid, ResultPath - use io_restart, only: ip_id, def_variable_2d, def_dim + use io_restart, only: icepack_files, icepack_path implicit none - type(t_mesh), target, intent(in) :: mesh - - integer, intent(in) :: year - integer (kind=int_kind) :: ncid + type(t_mesh) , intent(in) , target :: mesh + type(t_partit), intent(inout), target :: partit + integer , intent(in) :: year + logical , save :: has_been_called = .false. + integer (kind=int_kind) :: i, j, k, iblk, & ! counting indices nt_Tsfc, nt_sice, nt_qice, nt_qsno, & nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, & nt_vlvl, nt_iage, nt_FY, nt_aero, & ktherm, nt_fbri - integer (kind=int_kind) :: varid character(500) :: longname - character(500) :: filename character(500) :: trname, units character(4) :: cyear @@ -341,68 +346,57 @@ module subroutine init_restart_icepack(year, mesh) write(cyear,'(i4)') year ! Create an icepack restart file ! Only serial output implemented so far - ip_id%filename=trim(ResultPath)//trim(runid)//'.'//cyear//'.icepack.restart.nc' - if (ip_id%is_in_use) return - ip_id%is_in_use=.true. - - ! Define the dimensions of the netCDF file - ! - ! Note that at the moment FESOM2 supports only 3D output and restart (very - ! suboptimal). The different ice layers are thus splitted in different arrays - ! and - ! reconstructed after the restart. Multidimensional variables would solve - ! this. - - call def_dim(ip_id, 'node', nod2D) ! Number of nodes - call def_dim(ip_id, 'ncat', ncat) ! Number of thickness classes - + icepack_path=trim(ResultPath)//trim(runid)//'.'//cyear//'.icepack.restart.nc' + + if(has_been_called) return + has_been_called = .true. + ! Define the netCDF variables for surface ! and vertically constant fields !----------------------------------------------------------------- ! 3D restart fields (ncat) !----------------------------------------------------------------- - - call def_variable_2d(ip_id, 'aicen', (/nod2D, ncat/), 'sea ice concentration', 'none', aicen(:,:)); - call def_variable_2d(ip_id, 'vicen', (/nod2D, ncat/), 'volum per unit area of ice', 'm', vicen(:,:)); - call def_variable_2d(ip_id, 'vsnon', (/nod2D, ncat/), 'volum per unit area of snow', 'm', vsnon(:,:)); - call def_variable_2d(ip_id, 'Tsfc', (/nod2D, ncat/), 'sea ice surf. temperature', 'degC', trcrn(:,nt_Tsfc,:)); + call icepack_files%def_node_var('aicen' , 'sea ice concentration' , 'none', aicen(:,:) , ncat, mesh, partit) + call icepack_files%def_node_var('vicen' , 'volum per unit area of ice' , 'm' , vicen(:,:) , ncat, mesh, partit) + call icepack_files%def_node_var('vsnon' , 'volum per unit area of snow' , 'm' , vsnon(:,:) , ncat, mesh, partit) + call icepack_files%def_node_var('Tsfc' , 'sea ice surf. temperature' , 'degC', trcrn(:,nt_Tsfc,:), ncat, mesh, partit) if (tr_iage) then - call def_variable_2d(ip_id, 'iage', (/nod2D, ncat/), 'sea ice age', 's', trcrn(:,nt_iage,:)); + call icepack_files%def_node_var('iage' , 'sea ice age' , 's' , trcrn(:,nt_iage,:), ncat, mesh, partit) end if if (tr_FY) then - call def_variable_2d(ip_id, 'FY', (/nod2D, ncat/), 'first year ice', 'none', trcrn(:,nt_FY,:)); + call icepack_files%def_node_var('FY' , 'first year ice' , 'none', trcrn(:,nt_FY,:) , ncat, mesh, partit) end if if (tr_lvl) then - call def_variable_2d(ip_id, 'alvl', (/nod2D, ncat/), 'ridged sea ice area', 'none', trcrn(:,nt_alvl,:)); - call def_variable_2d(ip_id, 'vlvl', (/nod2D, ncat/), 'ridged sea ice volume', 'm', trcrn(:,nt_vlvl,:)); + call icepack_files%def_node_var('alvl' , 'ridged sea ice area' , 'none', trcrn(:,nt_alvl,:), ncat, mesh, partit) + call icepack_files%def_node_var('vlvl' , 'ridged sea ice volume' , 'm' , trcrn(:,nt_vlvl,:), ncat, mesh, partit) end if if (tr_pond_cesm) then - call def_variable_2d(ip_id, 'apnd', (/nod2D, ncat/), 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:)); - call def_variable_2d(ip_id, 'hpnd', (/nod2D, ncat/), 'melt pond depth', 'm', trcrn(:,nt_hpnd,:)); + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) end if if (tr_pond_topo) then - call def_variable_2d(ip_id, 'apnd', (/nod2D, ncat/), 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:)); - call def_variable_2d(ip_id, 'hpnd', (/nod2D, ncat/), 'melt pond depth', 'm', trcrn(:,nt_hpnd,:)); - call def_variable_2d(ip_id, 'ipnd', (/nod2D, ncat/), 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:)); + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), ncat, mesh, partit) end if if (tr_pond_lvl) then - call def_variable_2d(ip_id, 'apnd', (/nod2D, ncat/), 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:)); - call def_variable_2d(ip_id, 'hpnd', (/nod2D, ncat/), 'melt pond depth', 'm', trcrn(:,nt_hpnd,:)); - call def_variable_2d(ip_id, 'ipnd', (/nod2D, ncat/), 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:)); - call def_variable_2d(ip_id, 'ffracn', (/nod2D, ncat/), 'fraction of fsurfn over pond used to melt ipond', 'none', ffracn); - call def_variable_2d(ip_id, 'dhsn', (/nod2D, ncat/), 'depth difference for snow on sea ice and pond ice', 'm', dhsn); + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('ffracn', 'fraction of fsurfn over pond used to melt ipond' , 'none', ffracn, mesh, partit); + call icepack_files%def_node_var('dhsn' , 'depth difference for snow on sea ice and pond ice' , 'm' , dhsn , mesh, partit); end if if (tr_brine) then - call def_variable_2d(ip_id, 'fbri', (/nod2D, ncat/), 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:)); - call def_variable_2d(ip_id, 'first_ice', (/nod2D, ncat/), 'distinguishes ice that disappears', 'logical', first_ice_real(:,:)); + call icepack_files%def_node_var('fbri' , 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), ncat, mesh, partit) + call icepack_files%def_node_var('first_ice', 'distinguishes ice that disappears' , 'logical', first_ice_real(:,:), ncat, mesh, partit) end if !----------------------------------------------------------------- @@ -415,11 +409,11 @@ module subroutine init_restart_icepack(year, mesh) write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A21,i1)') 'sea ice salinity lyr:', k units='psu' - call def_variable_2d(ip_id, trim(trname), (/nod2D, ncat/), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:)); + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), ncat, mesh, partit) write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A21,i1)') 'sea ice enthalpy lyr:', k units='J/m3' - call def_variable_2d(ip_id, trim(trname), (/nod2D, ncat/), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:)); + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), ncat, mesh, partit) end do ! Snow @@ -428,7 +422,7 @@ module subroutine init_restart_icepack(year, mesh) write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A18,i1)') 'snow enthalpy lyr:', k units='J/m3' - call def_variable_2d(ip_id, trim(trname), (/nod2D, ncat/), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:)); + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), ncat, mesh, partit) end do ! @@ -438,6 +432,6 @@ module subroutine init_restart_icepack(year, mesh) ! enough to use these options. Lorenzo Zampieri - 16/10/2019. ! - end subroutine init_restart_icepack + end subroutine ini_icepack_io - end submodule icedrv_io +end submodule icedrv_io diff --git a/src/icepack_drivers/icedrv_main.F90 b/src/icepack_drivers/icedrv_main.F90 index 0f5f938ff..aef497d00 100644 --- a/src/icepack_drivers/icedrv_main.F90 +++ b/src/icepack_drivers/icedrv_main.F90 @@ -25,7 +25,7 @@ module icedrv_main ! Subroutines set_icepack, alloc_icepack, init_icepack, step_icepack, & icepack_to_fesom, icepack_to_fesom_single_point, & - init_flux_atm_ocn, init_io_icepack, init_restart_icepack + init_flux_atm_ocn, ini_icepack_io , ini_mean_icepack_io !======================================================================= !--------- Everything else is private @@ -883,19 +883,22 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) end subroutine step_icepack ! Initialize output - module subroutine init_io_icepack(mesh) + module subroutine ini_mean_icepack_io(mesh) use mod_mesh implicit none type(t_mesh), intent(in), target :: mesh - end subroutine init_io_icepack + end subroutine ini_mean_icepack_io ! Initialize restart - module subroutine init_restart_icepack(year, mesh) + module subroutine ini_icepack_io(year, partit, mesh) use mod_mesh + use mod_partit + use mod_parsup implicit none - type(t_mesh), intent(in), target :: mesh - integer(kind=int_kind), intent(in) :: year - end subroutine init_restart_icepack + type(t_mesh), intent(in) , target :: mesh + type(t_partit), intent(inout), target :: partit + integer(kind=int_kind), intent(in) :: year + end subroutine ini_icepack_io ! Cut off Icepack module subroutine cut_off_icepack diff --git a/src/icepack_drivers/icedrv_step.F90 b/src/icepack_drivers/icedrv_step.F90 index 5139454bf..917b7517c 100644 --- a/src/icepack_drivers/icedrv_step.F90 +++ b/src/icepack_drivers/icedrv_step.F90 @@ -772,12 +772,15 @@ end subroutine step_radiation ! Allows heat storage in ocean for uncoupled runs. ! - subroutine ocean_mixed_layer (dt) + subroutine ocean_mixed_layer (ice, dt) use icepack_intfc, only: icepack_ocn_mixed_layer, icepack_atm_boundary - + use mod_ice + implicit none + type(t_ice), target, intent(inout) :: ice + real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -852,7 +855,8 @@ subroutine ocean_mixed_layer (dt) !----------------------------------------------------------------- do i = 1, nx - call ocn_mixed_layer_icepack( alvdr_ocn=alvdr_ocn(i), swvdr=swvdr(i), & + call ocn_mixed_layer_icepack( ice, & + alvdr_ocn=alvdr_ocn(i), swvdr=swvdr(i), & alidr_ocn=alidr_ocn(i), swidr=swidr(i), & alvdf_ocn=alvdf_ocn(i), swvdf=swvdf(i), & alidf_ocn=alidf_ocn(i), swidf=swidf(i), & @@ -878,7 +882,7 @@ end subroutine ocean_mixed_layer !======================================================================= - subroutine ocn_mixed_layer_icepack( & + subroutine ocn_mixed_layer_icepack(ice, & alvdr_ocn, swvdr, & alidr_ocn, swidr, & alvdf_ocn, swvdf, & @@ -896,11 +900,11 @@ subroutine ocn_mixed_layer_icepack( & frzmlt, fsalt, & sss) - use i_therm_param, only: emiss_wat use g_forcing_param, only: use_virt_salt - + use mod_ice + implicit none - + real (kind=dbl_kind), intent(in) :: & alvdr_ocn , & ! visible, direct (fraction) alidr_ocn , & ! near-ir, direct (fraction) @@ -951,6 +955,10 @@ subroutine ocn_mixed_layer_icepack( & ice_ref_salinity character(len=*),parameter :: subname='(icepack_ocn_mixed_layer)' + + type(t_ice), target, intent(inout) :: ice + real(kind=WP), pointer :: emiss_wat + emiss_wat => ice%thermo%emiss_wat call icepack_query_parameters( Tffresh_out=Tffresh, Lfresh_out=Lfresh, & stefan_boltzmann_out=stefan_boltzmann, & @@ -992,11 +1000,14 @@ end subroutine ocn_mixed_layer_icepack !======================================================================= - subroutine coupling_prep(dt) - + subroutine coupling_prep(ice, dt) + + use mod_ice ! local variables implicit none + + type(t_ice), target, intent(inout) :: ice real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -1029,7 +1040,7 @@ subroutine coupling_prep(dt) frzmlt_init (i) = frzmlt(i) enddo - call ocean_mixed_layer (dt) ! ocean surface fluxes and sst + call ocean_mixed_layer (ice, dt) ! ocean surface fluxes and sst !----------------------------------------------------------------- ! Aggregate albedos @@ -1207,7 +1218,7 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) !----------------------------------------------------------------- dhi_dt(:) = ( vice(:) - dhi_dt(:) ) / dt - dhs_dt(:) = ( vsno(:) - dhi_dt(:) ) / dt + dhs_dt(:) = ( vsno(:) - dhs_dt(:) ) / dt ! clean up, update tendency diagnostics @@ -1294,7 +1305,7 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) ! get ready for coupling and the next time step !----------------------------------------------------------------- - call coupling_prep (dt) + call coupling_prep (ice, dt) !----------------------------------------------------------------- ! icepack timing diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 39dd8d247..2770b5345 100644 --- a/src/io_fesom_file.F90 +++ b/src/io_fesom_file.F90 @@ -43,9 +43,9 @@ module io_fesom_file_module procedure, public :: async_read_and_scatter_variables, async_gather_and_write_variables, join, init, is_iorank, rec_count, time_varindex, time_dimindex procedure, public :: read_variables_raw, write_variables_raw procedure, public :: close_file ! inherited procedures we overwrite - generic, public :: specify_node_var => specify_node_var_2d, specify_node_var_3d + generic, public :: specify_node_var => specify_node_var_2d, specify_node_var_3d, specify_node_var_2dicepack generic, public :: specify_elem_var => specify_elem_var_2d, specify_elem_var_3d - procedure, private :: specify_node_var_2d, specify_node_var_3d + procedure, private :: specify_node_var_2d, specify_node_var_3d, specify_node_var_2dicepack procedure, private :: specify_elem_var_2d, specify_elem_var_3d procedure, private :: read_and_scatter_variables, gather_and_write_variables end type @@ -54,6 +54,7 @@ module io_fesom_file_module integer, save :: m_nod2d integer, save :: m_elem2d integer, save :: m_nl + integer, save :: m_nitc type fesom_file_type_ptr @@ -103,7 +104,7 @@ function time_dimindex(this) result(x) end function - subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit) ! todo: would like to call it initialize but Fortran is rather cluncky with overwriting base type procedures + subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) ! todo: would like to call it initialize but Fortran is rather cluncky with overwriting base type procedures use io_netcdf_workaround_module use io_gather_module use MOD_PARTIT @@ -111,6 +112,7 @@ subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit) ! todo: would li integer mesh_nod2d integer mesh_elem2d integer mesh_nl + integer, optional :: mesh_nitc type(t_partit), target :: partit ! EO parameters type(fesom_file_type_ptr), allocatable :: tmparr(:) @@ -121,9 +123,16 @@ subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit) ! todo: would li call init_io_gather(partit) ! get hold of our mesh data for later use (assume the mesh instance will not change) - m_nod2d = mesh_nod2d + m_nod2d = mesh_nod2d m_elem2d = mesh_elem2d - m_nl = mesh_nl + m_nl = mesh_nl + !PS mesh_nitc ... icepack number of ice thickness classes, + if (present(mesh_nitc)) then + m_nitc = mesh_nitc + else + m_nitc = 0 + end if + call this%netcdf_file_type%initialize() allocate(this%used_mesh_dims(0)) @@ -446,6 +455,22 @@ subroutine specify_node_var_2d(this, name, longname, units, local_data) external_local_data_ptr(1:1,1:size(local_data)) => local_data(:) call specify_variable(this, name, [level_diminfo%idx, this%time_dimidx], level_diminfo%len, external_local_data_ptr, .false., longname, units) end subroutine + + subroutine specify_node_var_2dicepack(this, name, longname, units, local_data, nitc) + use, intrinsic :: ISO_C_BINDING + class(fesom_file_type), intent(inout) :: this + character(len=*), intent(in) :: name + character(len=*), intent(in) :: units, longname + real(kind=8), target, intent(inout) :: local_data(:,:) + integer, intent(in) :: nitc! todo: be able to set precision + ! EO parameters + type(dim_info) level_diminfo, nitc_diminfo + + level_diminfo = obtain_diminfo(this, m_nod2d) + nitc_diminfo = obtain_diminfo(this, size(local_data, dim=2)) + + call specify_variable(this, name, [level_diminfo%idx, nitc_diminfo%idx, this%time_dimidx], level_diminfo%len, local_data, .false., longname, units) + end subroutine subroutine specify_node_var_3d(this, name, longname, units, local_data) @@ -520,7 +545,9 @@ function obtain_diminfo(this, len) result(info) else if(len == m_nl-1) then info = dim_info( idx=this%add_dim('nz_1', len), len=len) else if(len == m_nl) then - info = dim_info( idx=this%add_dim('nz', len), len=len) + info = dim_info( idx=this%add_dim('nz' , len), len=len) + else if(len == m_nitc) then + info = dim_info( idx=this%add_dim('nitc', len), len=len) else print *, "error in line ",__LINE__, __FILE__," can not find dimension with size",len stop 1 diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index afc58f87f..cc08d77e4 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -1742,7 +1742,7 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) use iom #endif #if defined (__icepack) - use icedrv_main, only: init_io_icepack + use icedrv_main, only: ini_mean_icepack_io #endif implicit none integer :: istep @@ -1775,7 +1775,7 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) call init_io_gather(partit) #if defined (__icepack) - call init_io_icepack(mesh) !icapack has its copy of p_partit => partit + call ini_mean_icepack_io(mesh) !icapack has its copy of p_partit => partit #endif end if ! --> if (lfirst) then diff --git a/src/io_restart.F90 b/src/io_restart.F90 index fd1121278..482f098a3 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -13,6 +13,9 @@ MODULE io_RESTART use MOD_PARSUP use fortran_utils use mpi +#if defined(__icepack) + use icedrv_main +#endif #if defined(__recom) use recom_glovar use recom_config @@ -27,9 +30,17 @@ MODULE io_RESTART 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 character(:), allocatable, save :: oce_path + + type(restart_file_group) , save :: ice_files character(:), allocatable, save :: ice_path + +#if defined(__icepack) + type(restart_file_group) , save, public :: icepack_files + character(:), allocatable, save, public :: icepack_path + +#endif + #if defined(__recom) type(restart_file_group) , save :: bio_files character(:), allocatable, save :: bio_path @@ -233,9 +244,9 @@ end subroutine ini_bio_io ! subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tracers, partit, mesh) -#if defined(__icepack) - icepack restart not merged here ! produce a compiler error if USE_ICEPACK=ON; todo: merge icepack restart from 68d8b8b -#endif +! #if defined(__icepack) +! icepack restart not merged here ! produce a compiler error if USE_ICEPACK=ON; todo: merge icepack restart from 68d8b8b +! #endif use fortran_utils implicit none @@ -302,18 +313,29 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ! initialise files for netcdf restart if l_read==TRUE --> the restart file ! will be read if (.not. l_read) then - call ini_ocean_io(yearnew, dynamics, tracers, partit, mesh) - if (use_ice) call ini_ice_io (yearnew, ice, partit, mesh) + call ini_ocean_io(yearnew, dynamics, tracers, partit, mesh) + if (use_ice) then + call ini_ice_io(yearnew, ice, partit, mesh) +#if defined(__icepack) + call ini_icepack_io(yearnew, partit, mesh) +#endif + end if #if defined(__recom) if (use_REcoM) call ini_bio_io (yearnew, tracers, partit, mesh) #endif + else - call ini_ocean_io(yearold, dynamics, tracers, partit, mesh) - if (use_ice) call ini_ice_io (yearold, ice, partit, mesh) + call ini_ocean_io(yearold, dynamics, tracers, partit, mesh) + if (use_ice) then + call ini_ice_io (yearold, ice, partit, mesh) +#if defined(__icepack) + call ini_icepack_io(yearnew, partit, mesh) +#endif + end if #if defined(__recom) if (REcoM_restart) call ini_bio_io(yearold, tracers, partit, mesh) #endif - end if + end if ! --> if (.not. l_read) then !___READING OF RESTART________________________________________________________ ! should restart files be readed --> see r_restart in gen_modules_clock.F90 @@ -431,13 +453,16 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr if(use_ice) then if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ice'//achar(27)//'[0m' call write_restart(ice_path, ice_files, istep) +#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_restart(icepack_path, icepack_files, istep) +#endif end if #if defined(__recom) !RECOM restart !write here 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_restart(bio_path, bio_files, istep) end if diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index 3c7a327ca..4431c91eb 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -21,14 +21,14 @@ module restart_file_group_module integer, public :: nfiles = 0 ! todo: allow dynamically allocated size without messing with shallow copied pointers contains - generic, public :: def_node_var => def_node_var_2d, def_node_var_3d + generic, public :: def_node_var => def_node_var_2d, def_node_var_3d, def_node_var_2dicepack generic, public :: def_elem_var => def_elem_var_2d, def_elem_var_3d - procedure, private :: def_node_var_2d, def_node_var_3d + procedure, private :: def_node_var_2d, def_node_var_3d, def_node_var_2dicepack procedure, private :: def_elem_var_2d, def_elem_var_3d ! def_*_optional procedures create a restart variable which does not have to exist when reading the restart file - generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_3d_optional + generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_3d_optional, def_node_var_2dicepack_optional generic, public :: def_elem_var_optional => def_elem_var_2d_optional, def_elem_var_3d_optional - procedure, private :: def_node_var_2d_optional, def_node_var_3d_optional + procedure, private :: def_node_var_2d_optional, def_node_var_3d_optional, def_node_var_2dicepack_optional procedure, private :: def_elem_var_2d_optional, def_elem_var_3d_optional end type @@ -49,6 +49,23 @@ subroutine def_node_var_2d(this, name, longname, units, local_data, mesh, partit call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine + !PS add a seprate case for icepack since here i need to write a 2d vertice file + !PS over the dimension of number of ice thickness classes. Additional input + !PS parameter nitc + subroutine def_node_var_2dicepack(this, name, longname, units, local_data, nitc, mesh, partit) + use mod_mesh + class(restart_file_group), target, intent(inout) :: this + character(len=*), intent(in) :: name + character(len=*), intent(in) :: units, longname + real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + integer :: nitc + type(t_mesh) , intent(in) :: mesh + type(t_partit) , intent(in) :: partit + ! EO parameters + + call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + end subroutine subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit) use mod_mesh @@ -95,12 +112,16 @@ subroutine def_elem_var_3d(this, name, longname, units, local_data, mesh, partit end subroutine - subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_nl, partit) + subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) class(restart_file_group), target, intent(inout) :: g character(len=*), intent(in) :: name logical must_exist_on_read integer mesh_nod2d, mesh_elem2d, mesh_nl type(t_partit), intent(in) :: partit + !PS mesh_nitc ... icepack number of ice thickness classes, do it here as optional + !PS parameter, i assume is the less intrusive option. Therefore it must be + !PS at the end of the argumnent input list + integer, optional :: mesh_nitc ! EO parameters type(restart_file_type), pointer :: f @@ -111,7 +132,12 @@ subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_n f%path = "" allocate(f%varname,source=name) f%must_exist_on_read = must_exist_on_read - call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit) + + if (present(mesh_nitc)) then + call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) + else + call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit) + end if ! this is specific for a restart file f%iter_varindex = f%add_var_int('iter', [f%time_dimindex()]) end subroutine @@ -131,6 +157,23 @@ subroutine def_node_var_2d_optional(this, name, longname, units, local_data, mes call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine + !PS add a seprate case for icepack since here i need to write a 2d vertice file + !PS over the dimension of number of ice thickness classes. Additional input + !PS parameter nitc + subroutine def_node_var_2dicepack_optional(this, name, longname, units, local_data, nitc, mesh, partit) + use mod_mesh + class(restart_file_group), target, intent(inout) :: this + character(len=*), intent(in) :: name + character(len=*), intent(in) :: units, longname + real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + integer :: nitc + type(t_mesh) , intent(in) :: mesh + type(t_partit) , intent(in) :: partit + ! EO parameters + + call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + end subroutine subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit) use mod_mesh From a66454e37094a96b58118879f550a9f569eeb554 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 7 Oct 2024 17:05:05 +0200 Subject: [PATCH 02/18] add old icepack debug changes into ice_oce_coupling.F90 --- src/ice_oce_coupling.F90 | 31 ++++++++++++++++++------- src/icepack_drivers/icedrv_main.F90 | 6 +++-- src/icepack_drivers/icedrv_transfer.F90 | 7 ++++-- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 0d4575668..499502a77 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -362,16 +362,22 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) fhocn_tot_out=net_heat_flux, & fresh_tot_out=fresh_wa_flux, & fsalt_out=real_salt_flux, & - dhi_dt_out=thdgrsn, & - dhs_dt_out=thdgr, & - evap_ocn_out=evaporation ) + dhs_dt_out=thdgrsn, & + dhi_dt_out=thdgr, & + evap_ocn_out=evaporation, & + evap_out=ice_sublimation ) - heat_flux(:) = - net_heat_flux(:) - water_flux(:) = - (fresh_wa_flux(:)/1000.0_WP) - runoff(:) + ! Heat flux + heat_flux = - net_heat_flux - ! Evaporation - evaporation(:) = - evaporation(:) / 1000.0_WP - ice_sublimation(:) = 0.0_WP + ! Freshwater flux (convert units from icepack to fesom) + water_flux = - (fresh_wa_flux * inv_rhowat) - runoff(:) + + ! Evaporation (convert units from icepack to fesom) + evaporation = - evaporation * (1.0_WP - a_ice) * inv_rhowat + + ! Ice-Sublimation is not added to evap in icepack + ice_sublimation = 0.0_WP call init_flux_atm_ocn() @@ -509,12 +515,21 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) ! enforce the total freshwater/salt flux be zero ! 1. water flux ! if (.not. use_virt_salt) can be used! ! we conserve only the fluxes from the database plus evaporation. + ! ICEPACK: adds rain, snow and evap is based on the newly formed ice + ! concentration (a_ice). In our standard ice model rain, snow and evap is + ! added based on the ice concentration of the previous time step (a_ice_old) + ! So for the proper balancing of snow the proper aice has to be choosen !$OMP PARALLEL DO do n=1, myDim_nod2D+eDim_nod2D flux(n) = evaporation(n) & -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean +prec_rain(n) & +#if defined (__icepack) + +prec_snow(n)*(1.0_WP-a_ice(n)) & +#else +prec_snow(n)*(1.0_WP-a_ice_old(n)) & +#endif + #if defined (__oasis) || defined (__ifsinterface) +residualifwflx(n) & ! balance residual ice flux only in coupled case #endif diff --git a/src/icepack_drivers/icedrv_main.F90 b/src/icepack_drivers/icedrv_main.F90 index aef497d00..d749964bf 100644 --- a/src/icepack_drivers/icedrv_main.F90 +++ b/src/icepack_drivers/icedrv_main.F90 @@ -820,7 +820,8 @@ module subroutine icepack_to_fesom( & fhocn_tot_out, fresh_tot_out, & strocnxT_out, strocnyT_out, & dhs_dt_out, dhi_dt_out, & - fsalt_out, evap_ocn_out ) + fsalt_out, evap_ocn_out, & + evap_out ) use mod_mesh implicit none integer (kind=int_kind), intent(in) :: & @@ -836,7 +837,8 @@ module subroutine icepack_to_fesom( & fsalt_out, & dhs_dt_out, & dhi_dt_out, & - evap_ocn_out + evap_ocn_out, & + evap_out end subroutine icepack_to_fesom ! Copy variables from icepack to fesom (single node or element) diff --git a/src/icepack_drivers/icedrv_transfer.F90 b/src/icepack_drivers/icedrv_transfer.F90 index 15d36eb1b..524b9d9ec 100644 --- a/src/icepack_drivers/icedrv_transfer.F90 +++ b/src/icepack_drivers/icedrv_transfer.F90 @@ -188,7 +188,8 @@ module subroutine icepack_to_fesom( nx_in, & fhocn_tot_out, fresh_tot_out, & strocnxT_out, strocnyT_out, & dhs_dt_out, dhi_dt_out, & - fsalt_out, evap_ocn_out ) + fsalt_out, evap_ocn_out, & + evap_out) implicit none @@ -206,7 +207,8 @@ module subroutine icepack_to_fesom( nx_in, & fsalt_out, & dhs_dt_out, & dhi_dt_out, & - evap_ocn_out + evap_ocn_out, & + evap_out character(len=*),parameter :: subname='(icepack_to_fesom)' @@ -222,6 +224,7 @@ module subroutine icepack_to_fesom( nx_in, & if (present(dhs_dt_out) ) dhs_dt_out = dhs_dt if (present(fsalt_out) ) fsalt_out = fsalt if (present(evap_ocn_out) ) evap_ocn_out = evap_ocn + if (present(evap_ocn_out) ) evap_out = evap end subroutine icepack_to_fesom From da877f3523f5c3ab3c28a2d692534d4891b1435e Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 9 Oct 2024 13:30:22 +0200 Subject: [PATCH 03/18] switch defalt icepack update ocean fluxes to true --- config/namelist.icepack | 2 +- config/namelist.icepack.cesm.ponds | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/config/namelist.icepack b/config/namelist.icepack index ed0dd4d4c..3fa487147 100644 --- a/config/namelist.icepack +++ b/config/namelist.icepack @@ -89,7 +89,7 @@ ustar_min = 0.0005 emissivity = 0.95 fbot_xfer_type = 'constant' - update_ocn_f = .false. + update_ocn_f = .true. l_mpond_fresh = .false. tfrz_option = 'linear_salt' oceanmixed_ice = .true. diff --git a/config/namelist.icepack.cesm.ponds b/config/namelist.icepack.cesm.ponds index 51aa33191..9e3bc18df 100644 --- a/config/namelist.icepack.cesm.ponds +++ b/config/namelist.icepack.cesm.ponds @@ -87,7 +87,7 @@ ustar_min = 0.0005 emissivity = 0.95 fbot_xfer_type = 'constant' - update_ocn_f = .false. + update_ocn_f = .true. l_mpond_fresh = .false. tfrz_option = 'linear_salt' oceanmixed_ice = .true. From e62e2c23870d1841e7edbb77d934baf13cd9067f Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 9 Oct 2024 13:31:44 +0200 Subject: [PATCH 04/18] make icepack restart work for the moment --- src/io_fesom_file.F90 | 21 ++++++++++++++------- src/io_restart_file_group.F90 | 22 +++++++++++----------- 2 files changed, 25 insertions(+), 18 deletions(-) diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 2770b5345..0df58c7da 100644 --- a/src/io_fesom_file.F90 +++ b/src/io_fesom_file.F90 @@ -43,9 +43,9 @@ module io_fesom_file_module procedure, public :: async_read_and_scatter_variables, async_gather_and_write_variables, join, init, is_iorank, rec_count, time_varindex, time_dimindex procedure, public :: read_variables_raw, write_variables_raw procedure, public :: close_file ! inherited procedures we overwrite - generic, public :: specify_node_var => specify_node_var_2d, specify_node_var_3d, specify_node_var_2dicepack + generic, public :: specify_node_var => specify_node_var_2d, specify_node_var_2dicepack, specify_node_var_3d generic, public :: specify_elem_var => specify_elem_var_2d, specify_elem_var_3d - procedure, private :: specify_node_var_2d, specify_node_var_3d, specify_node_var_2dicepack + procedure, private :: specify_node_var_2d, specify_node_var_2dicepack, specify_node_var_3d procedure, private :: specify_elem_var_2d, specify_elem_var_3d procedure, private :: read_and_scatter_variables, gather_and_write_variables end type @@ -450,6 +450,7 @@ subroutine specify_node_var_2d(this, name, longname, units, local_data) real(8), pointer :: external_local_data_ptr(:,:) type(dim_info) level_diminfo +!PS write(*,*) "--> specify_node_var_2d:", __LINE__, __FILE__ level_diminfo = obtain_diminfo(this, m_nod2d) external_local_data_ptr(1:1,1:size(local_data)) => local_data(:) @@ -466,6 +467,7 @@ subroutine specify_node_var_2dicepack(this, name, longname, units, local_data, n ! EO parameters type(dim_info) level_diminfo, nitc_diminfo +!PS write(*,*) "--> specify_node_var_2dicepack:", __LINE__, __FILE__ level_diminfo = obtain_diminfo(this, m_nod2d) nitc_diminfo = obtain_diminfo(this, size(local_data, dim=2)) @@ -481,7 +483,8 @@ subroutine specify_node_var_3d(this, name, longname, units, local_data) real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision ! EO parameters type(dim_info) level_diminfo, depth_diminfo - + +!PS write(*,*) "--> specify_node_var_3d:", __LINE__, __FILE__ level_diminfo = obtain_diminfo(this, m_nod2d) depth_diminfo = obtain_diminfo(this, size(local_data, dim=1)) @@ -499,6 +502,7 @@ subroutine specify_elem_var_2d(this, name, longname, units, local_data) real(8), pointer :: external_local_data_ptr(:,:) type(dim_info) level_diminfo +!PS write(*,*) "--> specify_elem_var_2d:", __LINE__, __FILE__ level_diminfo = obtain_diminfo(this, m_elem2d) external_local_data_ptr(1:1,1:size(local_data)) => local_data(:) @@ -515,6 +519,7 @@ subroutine specify_elem_var_3d(this, name, longname, units, local_data) ! EO parameters type(dim_info) level_diminfo, depth_diminfo +!PS write(*,*) "--> specify_elem_var_3d:", __LINE__, __FILE__ level_diminfo = obtain_diminfo(this, m_elem2d) depth_diminfo = obtain_diminfo(this, size(local_data, dim=1)) @@ -537,16 +542,18 @@ function obtain_diminfo(this, len) result(info) end if end do +!PS write(*,*) 'm_n2d=',m_nod2d,', m_e2d=',m_elem2d,', m_nl=', m_nl,', nitc', m_nitc, ', LEN=', len + ! the dim has not been added yet, see if it is one of our allowed mesh related dims - if(len == m_nod2d) then + if (len == m_nod2d) then info = dim_info( idx=this%add_dim('node', len), len=len) else if(len == m_elem2d) then info = dim_info( idx=this%add_dim('elem', len), len=len) - else if(len == m_nl-1) then + else if(len == m_nl-1 ) then info = dim_info( idx=this%add_dim('nz_1', len), len=len) - else if(len == m_nl) then + else if(len == m_nl ) then info = dim_info( idx=this%add_dim('nz' , len), len=len) - else if(len == m_nitc) then + else if(len == m_nitc ) then info = dim_info( idx=this%add_dim('nitc', len), len=len) else print *, "error in line ",__LINE__, __FILE__," can not find dimension with size",len diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index 4431c91eb..bbd139805 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -21,14 +21,14 @@ module restart_file_group_module integer, public :: nfiles = 0 ! todo: allow dynamically allocated size without messing with shallow copied pointers contains - generic, public :: def_node_var => def_node_var_2d, def_node_var_3d, def_node_var_2dicepack + generic, public :: def_node_var => def_node_var_2d, def_node_var_2dicepack, def_node_var_3d generic, public :: def_elem_var => def_elem_var_2d, def_elem_var_3d - procedure, private :: def_node_var_2d, def_node_var_3d, def_node_var_2dicepack + procedure, private :: def_node_var_2d, def_node_var_2dicepack, def_node_var_3d procedure, private :: def_elem_var_2d, def_elem_var_3d ! def_*_optional procedures create a restart variable which does not have to exist when reading the restart file - generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_3d_optional, def_node_var_2dicepack_optional + generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_2dicepack_optional, def_node_var_3d_optional generic, public :: def_elem_var_optional => def_elem_var_2d_optional, def_elem_var_3d_optional - procedure, private :: def_node_var_2d_optional, def_node_var_3d_optional, def_node_var_2dicepack_optional + procedure, private :: def_node_var_2d_optional, def_node_var_2dicepack_optional, def_node_var_3d_optional procedure, private :: def_elem_var_2d_optional, def_elem_var_3d_optional end type @@ -44,7 +44,7 @@ subroutine def_node_var_2d(this, name, longname, units, local_data, mesh, partit type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters - +!PS write(*,*) "--> def_node_var_2d:", __LINE__, __FILE__ call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine @@ -62,9 +62,9 @@ subroutine def_node_var_2dicepack(this, name, longname, units, local_data, nitc, type(t_mesh) , intent(in) :: mesh type(t_partit) , intent(in) :: partit ! EO parameters - +!PS write(*,*) "--> def_node_var_2dicepack:", __LINE__, __FILE__ call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, nitc) end subroutine subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit) @@ -76,7 +76,7 @@ subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters - +!PS write(*,*) "--> def_node_var_3d:", __LINE__, __FILE__ call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine @@ -91,7 +91,7 @@ subroutine def_elem_var_2d(this, name, longname, units, local_data, mesh, partit type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters - +!PS write(*,*) "--> def_elem_var_2d:", __LINE__, __FILE__ call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) call this%files(this%nfiles)%specify_elem_var(name, longname, units, local_data) end subroutine @@ -106,7 +106,7 @@ subroutine def_elem_var_3d(this, name, longname, units, local_data, mesh, partit type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters - +!PS write(*,*) "--> def_elem_var_3d:", __LINE__, __FILE__ call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) call this%files(this%nfiles)%specify_elem_var(name, longname, units, local_data) end subroutine @@ -172,7 +172,7 @@ subroutine def_node_var_2dicepack_optional(this, name, longname, units, local_da ! EO parameters call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, nitc) end subroutine subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit) From 2d0b474f2e1916e7dd248ba0dd1b2d6e6faa4dbb Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 9 Oct 2024 13:40:53 +0200 Subject: [PATCH 05/18] add icepack bugfix from ancient debugging round in 2021, Problem in moment Volume is not conserved --- src/ice_oce_coupling.F90 | 17 ++- src/icepack_drivers/icedrv_allocate.F90 | 4 + src/icepack_drivers/icedrv_main.F90 | 4 + src/icepack_drivers/icedrv_set.F90 | 179 +++++++++++++----------- src/icepack_drivers/icedrv_step.F90 | 54 +++++-- src/icepack_drivers/icedrv_transfer.F90 | 22 +-- 6 files changed, 172 insertions(+), 108 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 499502a77..d9edf1eec 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -376,8 +376,8 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) ! Evaporation (convert units from icepack to fesom) evaporation = - evaporation * (1.0_WP - a_ice) * inv_rhowat - ! Ice-Sublimation is not added to evap in icepack - ice_sublimation = 0.0_WP +!PS ! Ice-Sublimation is not added to evap in icepack +!PS ice_sublimation = 0.0_WP call init_flux_atm_ocn() @@ -519,14 +519,25 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) ! concentration (a_ice). In our standard ice model rain, snow and evap is ! added based on the ice concentration of the previous time step (a_ice_old) ! So for the proper balancing of snow the proper aice has to be choosen + ! -icepack_therm_itd.F90 --> subroutine icepack_step_therm2(...) + ! fresh = fresh + frain*aice + ! -icedrv_step.F90: subroutine ocn_mixed_layer_icepack(... + ! fresh_tot = fresh + (-evap_ocn + frain + fsnow)*(c1-aice) + ! At the end all rain is added to the ocean, only snow needs to be + ! scaled with (1-aice ) + ! -Ice-Sublimation is not added to evap in icepack, therefor we dont need + ! to compensate for it the ice2atmos subplimation does not contribute + ! to the freshwater flux into the ocean + !$OMP PARALLEL DO do n=1, myDim_nod2D+eDim_nod2D flux(n) = evaporation(n) & - -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean +prec_rain(n) & + #if defined (__icepack) +prec_snow(n)*(1.0_WP-a_ice(n)) & #else + -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean +prec_snow(n)*(1.0_WP-a_ice_old(n)) & #endif diff --git a/src/icepack_drivers/icedrv_allocate.F90 b/src/icepack_drivers/icedrv_allocate.F90 index e2d6d6108..1faccbc22 100644 --- a/src/icepack_drivers/icedrv_allocate.F90 +++ b/src/icepack_drivers/icedrv_allocate.F90 @@ -116,6 +116,10 @@ subroutine alloc_flux opening(nx) , & ! rate of opening due to divergence/shear (1/s) dhi_dt(nx) , & ! ice volume tendency due to thermodynamics (m/s) dhs_dt(nx) , & ! snow volume tendency due to thermodynamics (m/s) + dhi_t_dt(nx) , & ! ice volume tendency due to thermodynamics (m/s) + dhs_t_dt(nx) , & ! snow volume tendency due to thermodynamics (m/s) + dhi_r_dt(nx) , & ! ice volume tendency due to ridging (m/s) + dhs_r_dt(nx) , & ! snow volume tendency due to ridging (m/s) dardg1ndt(nx,ncat), & ! rate of area loss by ridging ice (1/s) dardg2ndt(nx,ncat), & ! rate of area gain by new ridges (1/s) dvirdgndt(nx,ncat), & ! rate of ice volume ridged (m/s) diff --git a/src/icepack_drivers/icedrv_main.F90 b/src/icepack_drivers/icedrv_main.F90 index d749964bf..3c609508c 100644 --- a/src/icepack_drivers/icedrv_main.F90 +++ b/src/icepack_drivers/icedrv_main.F90 @@ -153,6 +153,10 @@ module icedrv_main opening(:), & ! rate of opening due to divergence/shear (1/s) dhi_dt(:), & ! ice volume tendency due to thermodynamics (m/s) dhs_dt(:), & ! snow volume tendency due to thermodynamics (m/s) + dhi_t_dt(:), & ! ice volume tendency due to thermodynamics (m/s) + dhs_t_dt(:), & ! snow volume tendency due to thermodynamics (m/s) + dhi_r_dt(:), & ! ice volume tendency due to ridging (m/s) + dhs_r_dt(:), & ! snow volume tendency due to ridging (m/s) ! ridging diagnostics in categories dardg1ndt(:,:), & ! rate of area loss by ridging ice (1/s) dardg2ndt(:,:), & ! rate of area gain by new ridges (1/s) diff --git a/src/icepack_drivers/icedrv_set.F90 b/src/icepack_drivers/icedrv_set.F90 index dbfc55a1d..63bd5d93f 100644 --- a/src/icepack_drivers/icedrv_set.F90 +++ b/src/icepack_drivers/icedrv_set.F90 @@ -24,9 +24,6 @@ module subroutine set_icepack(ice, partit) use MOD_ICE -! use i_param, only: whichEVP -! use i_param, only: cd_oce_ice, Pstar, c_pressure -! use i_therm_param, only: albw implicit none @@ -153,7 +150,10 @@ module subroutine set_icepack(ice, partit) logical (kind=log_kind) :: oceanmixed_ice character (len=char_len) :: wave_spec_type - + ! to snychronize density definition between icepack and fesom2, becomes + ! crucial for waterflux computation and volume conservation + real (kind=dbl_kind) :: rhoice, rhosno, rhowat, rhofwt + !----------------------------------------------------------------- ! Namelist definition !----------------------------------------------------------------- @@ -238,16 +238,16 @@ module subroutine set_icepack(ice, partit) end if do while (nml_error > 0) - if (mype == 0) print*,'Reading namelist file ',nml_filename + if (partit%mype == 0) print*,'Reading namelist file ',nml_filename - if (mype == 0) print*,'Reading env_nml' + if (partit%mype == 0) print*,'Reading env_nml' read(nu_nml, nml=env_nml,iostat=nml_error) if (nml_error /= 0) exit end do if (nml_error == 0) close(nu_nml) if (nml_error /= 0) then - if (mype == 0) write(*,*) 'Error reading env namelist' + if (partit%mype == 0) write(*,*) 'Error reading env namelist' call icedrv_system_abort(file=__FILE__,line=__LINE__) close(nu_nml) end if @@ -339,7 +339,25 @@ module subroutine set_icepack(ice, partit) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) - + + !----------------------------------------------------------------- + ! Synchronize values between Icepack and FESOM2 values + !----------------------------------------------------------------- + ! Make the namelists.ice and namelist.icepack consistent (icepack wins + ! over fesom) + ice%cd_oce_ice = dragio + ice%thermo%albw= albocn + ice%Pstar = P_star + ice%c_pressure = C_star + + ! in terms of density definition fesom wins over icepack, otherwise + ! we can get in trouble with the waterflux computation and thus the + ! volume conservation under zstar + rhoice = ice%thermo%rhoice + rhosno = ice%thermo%rhosno + rhowat = ice%thermo%rhowat + rhofwt = ice%thermo%rhofwt + !----------------------------------------------------------------- ! other default values !----------------------------------------------------------------- @@ -402,7 +420,7 @@ module subroutine set_icepack(ice, partit) if (nml_error == 0) close(nu_nml) if (nml_error /= 0) then - if (mype == 0) write(*,*) 'Error reading iecpack namelists' + if (partit%mype == 0) write(*,*) 'Error reading iecpack namelists' call icedrv_system_abort(file=__FILE__,line=__LINE__) endif close(nu_nml) @@ -411,47 +429,47 @@ module subroutine set_icepack(ice, partit) ! set up diagnostics output and resolve conflicts !----------------------------------------------------------------- - if (mype == 0) write(*,*) 'Diagnostic output will be in file ' - if (mype == 0) write(*,*) ' icepack.diagnostics' + if (partit%mype == 0) write(*,*) 'Diagnostic output will be in file ' + if (partit%mype == 0) write(*,*) ' icepack.diagnostics' diag_filename = 'icepack.errors' open (ice_stderr, file=diag_filename, status='unknown', iostat=diag_error) if (diag_error /= 0) then - if (mype == 0) write(*,*) 'Error while opening error file' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write(*,*) 'Error while opening error file' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif diag_filename = 'icepack.diagnostics' open (nu_diag, file=diag_filename, status='unknown', iostat=diag_error) if (diag_error /= 0) then - if (mype == 0) write(*,*) 'Error while opening diagnostic file' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write(*,*) 'Error while opening diagnostic file' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif - if (mype == 0) write(nu_diag,*) '-----------------------------------' - if (mype == 0) write(nu_diag,*) ' ICEPACK model diagnostic output ' - if (mype == 0) write(nu_diag,*) '-----------------------------------' - if (mype == 0) write(nu_diag,*) ' ' + if (partit%mype == 0) write(nu_diag,*) '-----------------------------------' + if (partit%mype == 0) write(nu_diag,*) ' ICEPACK model diagnostic output ' + if (partit%mype == 0) write(nu_diag,*) '-----------------------------------' + if (partit%mype == 0) write(nu_diag,*) ' ' if (ice%whichEVP == 1 .or. ice%whichEVP == 2) then - if (mype == 0) write (nu_diag,*) 'WARNING: whichEVP = 1 or 2' - if (mype == 0) write (nu_diag,*) 'Adaptive or Modified EVP formulations' - if (mype == 0) write (nu_diag,*) 'are not allowed when using Icepack (yet).' - if (mype == 0) write (nu_diag,*) 'Standard EVP will be used instead' - if (mype == 0) write (nu_diag,*) ' whichEVP = 0' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: whichEVP = 1 or 2' + if (partit%mype == 0) write (nu_diag,*) 'Adaptive or Modified EVP formulations' + if (partit%mype == 0) write (nu_diag,*) 'are not allowed when using Icepack (yet).' + if (partit%mype == 0) write (nu_diag,*) 'Standard EVP will be used instead' + if (partit%mype == 0) write (nu_diag,*) ' whichEVP = 0' ice%whichEVP = 0 endif if (ncat == 1 .and. kitd == 1) then - if (mype == 0) write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1.' - if (mype == 0) write (nu_diag,*) 'Use kitd = 0 (delta function ITD) with kcatbound = 0' - if (mype == 0) write (nu_diag,*) 'or for column configurations use kcatbound = -1' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1.' + if (partit%mype == 0) write (nu_diag,*) 'Use kitd = 0 (delta function ITD) with kcatbound = 0' + if (partit%mype == 0) write (nu_diag,*) 'or for column configurations use kcatbound = -1' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif if (ncat /= 1 .and. kcatbound == -1) then - if (mype == 0) write (nu_diag,*) 'WARNING: ITD required for ncat > 1' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting kitd and kcatbound to default values' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: ITD required for ncat > 1' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting kitd and kcatbound to default values' kitd = 1 kcatbound = 0 endif @@ -467,46 +485,46 @@ module subroutine set_icepack(ice, partit) if (rpcesm + rplvl + rptopo > puny) tr_pond = .true. if (rpcesm + rplvl + rptopo > c1 + puny) then - if (mype == 0) write (nu_diag,*) 'WARNING: Must use only one melt pond scheme' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Must use only one melt pond scheme' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif if (tr_pond_lvl .and. .not. tr_lvl) then - if (mype == 0) write (nu_diag,*) 'WARNING: tr_pond_lvl=T but tr_lvl=F' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting tr_lvl=T' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: tr_pond_lvl=T but tr_lvl=F' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting tr_lvl=T' tr_lvl = .true. endif if (tr_pond_lvl .and. abs(hs0) > puny) then - if (mype == 0) write (nu_diag,*) 'WARNING: tr_pond_lvl=T and hs0/=0' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting hs0=0' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: tr_pond_lvl=T and hs0/=0' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting hs0=0' hs0 = c0 endif if (tr_pond_cesm .and. trim(frzpnd) /= 'cesm') then - if (mype == 0) write (nu_diag,*) 'WARNING: tr_pond_cesm=T' - if (mype == 0) write (nu_diag,*) 'WARNING: frzpnd, dpscale not used' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: tr_pond_cesm=T' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: frzpnd, dpscale not used' frzpnd = 'cesm' endif if (trim(shortwave) /= 'dEdd' .and. tr_pond .and. calc_tsfc) then - if (mype == 0) write (nu_diag,*) 'WARNING: Must use dEdd shortwave' - if (mype == 0) write (nu_diag,*) 'WARNING: with tr_pond and calc_tsfc=T.' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting shortwave = dEdd' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Must use dEdd shortwave' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: with tr_pond and calc_tsfc=T.' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting shortwave = dEdd' shortwave = 'dEdd' endif if (tr_aero .and. n_aero==0) then - if (mype == 0) write (nu_diag,*) 'WARNING: aerosols activated but' - if (mype == 0) write (nu_diag,*) 'WARNING: not allocated in tracer array.' - if (mype == 0) write (nu_diag,*) 'WARNING: Activate in compilation script.' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write (nu_diag,*) 'WARNING: aerosols activated but' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: not allocated in tracer array.' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Activate in compilation script.' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif if (tr_aero .and. trim(shortwave) /= 'dEdd') then - if (mype == 0) write (nu_diag,*) 'WARNING: aerosols activated but dEdd' - if (mype == 0) write (nu_diag,*) 'WARNING: shortwave is not.' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting shortwave = dEdd' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: aerosols activated but dEdd' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: shortwave is not.' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting shortwave = dEdd' shortwave = 'dEdd' endif @@ -514,53 +532,53 @@ module subroutine set_icepack(ice, partit) rfracmax = min(max(rfracmax,c0),c1) if (ktherm == 2 .and. .not. calc_Tsfc) then - if (mype == 0) write (nu_diag,*) 'WARNING: ktherm = 2 and calc_Tsfc = F' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting calc_Tsfc = T' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: ktherm = 2 and calc_Tsfc = F' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting calc_Tsfc = T' calc_Tsfc = .true. endif if (ktherm == 1 .and. trim(tfrz_option) /= 'linear_salt') then - if (mype == 0) write (nu_diag,*) & + if (partit%mype == 0) write (nu_diag,*) & 'WARNING: ktherm = 1 and tfrz_option = ',trim(tfrz_option) - if (mype == 0) write (nu_diag,*) & + if (partit%mype == 0) write (nu_diag,*) & 'WARNING: For consistency, set tfrz_option = linear_salt' endif if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then - if (mype == 0) write (nu_diag,*) & + if (partit%mype == 0) write (nu_diag,*) & 'WARNING: ktherm = 2 and tfrz_option = ',trim(tfrz_option) - if (mype == 0) write (nu_diag,*) & + if (partit%mype == 0) write (nu_diag,*) & 'WARNING: For consistency, set tfrz_option = mushy' endif if (formdrag) then if (trim(atmbndy) == 'constant') then - if (mype == 0) write (nu_diag,*) 'WARNING: atmbndy = constant not allowed with formdrag' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting atmbndy = default' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: atmbndy = constant not allowed with formdrag' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting atmbndy = default' atmbndy = 'default' endif if (.not. calc_strair) then - if (mype == 0) write (nu_diag,*) 'WARNING: formdrag=T but calc_strair=F' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting calc_strair=T' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: formdrag=T but calc_strair=F' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting calc_strair=T' calc_strair = .true. endif if (tr_pond_cesm) then - if (mype == 0) write (ice_stderr,*) 'ERROR: formdrag=T but frzpnd=cesm' - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write (ice_stderr,*) 'ERROR: formdrag=T but frzpnd=cesm' + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif if (.not. tr_lvl) then - if (mype == 0) write (nu_diag,*) 'WARNING: formdrag=T but tr_lvl=F' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting tr_lvl=T' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: formdrag=T but tr_lvl=F' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting tr_lvl=T' tr_lvl = .true. max_ntrcr = max_ntrcr + 2 ! tr_lvl brings two more tracers endif endif if (trim(fbot_xfer_type) == 'Cdn_ocn' .and. .not. formdrag) then - if (mype == 0) write (nu_diag,*) 'WARNING: formdrag=F but fbot_xfer_type=Cdn_ocn' - if (mype == 0) write (nu_diag,*) 'WARNING: Setting fbot_xfer_type = constant' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: formdrag=F but fbot_xfer_type=Cdn_ocn' + if (partit%mype == 0) write (nu_diag,*) 'WARNING: Setting fbot_xfer_type = constant' fbot_xfer_type = 'constant' endif @@ -571,7 +589,7 @@ module subroutine set_icepack(ice, partit) ! Write Icepack configuration !----------------------------------------------------------------- - if (mype == 0) then + if (partit%mype == 0) then write(nu_diag,*) ' Document ice_in namelist parameters:' write(nu_diag,*) ' ==================================== ' @@ -616,6 +634,11 @@ module subroutine set_icepack(ice, partit) write(nu_diag,1000) ' ahmax = ', ahmax endif + write(nu_diag,1000) ' rhos = ', rhosno + write(nu_diag,1000) ' rhoi = ', rhoice + write(nu_diag,1000) ' rhow = ', rhowat + write(nu_diag,1000) ' rhofwt = ', rhofwt + write(nu_diag,1000) ' rfracmin = ', rfracmin write(nu_diag,1000) ' rfracmax = ', rfracmax @@ -656,7 +679,7 @@ module subroutine set_icepack(ice, partit) write(nu_diag,1010) ' update_ocn_f = ', update_ocn_f write(nu_diag,1010) ' wave_spec = ', wave_spec write(nu_diag,1005) ' dragio = ', dragio - write(nu_diag,1005) ' Pstar = ', P_star + write(nu_diag,* ) ' Pstar = ', P_star write(nu_diag,1005) ' Cstar = ', C_star if (kstrength==1) then @@ -760,12 +783,12 @@ module subroutine set_icepack(ice, partit) endif if (ntrcr > max_ntrcr-1) then - if (mype == 0) write(ice_stderr,*) 'max_ntrcr-1 < number of namelist tracers' - if (mype == 0) write(ice_stderr,*) 'max_ntrcr-1 = ',max_ntrcr-1,' ntrcr = ',ntrcr - if (mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) + if (partit%mype == 0) write(ice_stderr,*) 'max_ntrcr-1 < number of namelist tracers' + if (partit%mype == 0) write(ice_stderr,*) 'max_ntrcr-1 = ',max_ntrcr-1,' ntrcr = ',ntrcr + if (partit%mype == 0) call icedrv_system_abort(file=__FILE__,line=__LINE__) endif - if (mype == 0) then + if (partit%mype == 0) then write(nu_diag,*) ' ' write(nu_diag,1020) 'max_ntrcr = ', max_ntrcr @@ -813,17 +836,7 @@ module subroutine set_icepack(ice, partit) !----------------------------------------------------------------- ! set Icepack values - !----------------------------------------------------------------- - - ! Make the namelists.ice and namelist.icepack consistent (icepack wins - ! over fesom) - - ice%cd_oce_ice = dragio - ice%thermo%albw= albocn - ice%Pstar = P_star - ice%c_pressure = C_star - - + !----------------------------------------------------------------- call icepack_init_parameters(ustar_min_in=ustar_min, Cf_in=Cf, & albicev_in=albicev, albicei_in=albicei, & albsnowv_in=albsnowv, albsnowi_in=albsnowi, & @@ -850,7 +863,9 @@ module subroutine set_icepack(ice, partit) fbot_xfer_type_in=fbot_xfer_type, & wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, & ksno_in=ksno, dragio_in=dragio, Cstar_in=C_star, & - Pstar_in=P_star, albocn_in=albocn) + Pstar_in=P_star, albocn_in=albocn, & + update_ocn_f_in=update_ocn_f, & + rhoi_in=rhoice, rhos_in=rhosno, rhow_in=rhowat, rhofresh_in=rhofwt) call icepack_init_tracer_sizes(ntrcr_in=ntrcr, & ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & nfsd_in=nfsd, n_aero_in=n_aero) diff --git a/src/icepack_drivers/icedrv_step.F90 b/src/icepack_drivers/icedrv_step.F90 index 917b7517c..826134f36 100644 --- a/src/icepack_drivers/icedrv_step.F90 +++ b/src/icepack_drivers/icedrv_step.F90 @@ -313,7 +313,10 @@ subroutine step_therm2 (dt) logical (kind=log_kind) :: & tr_fsd ! floe size distribution tracers - + + logical (kind=log_kind) :: & + update_ocn_f ! floe size distribution tracers + character(len=*), parameter :: subname='(step_therm2)' !----------------------------------------------------------------- @@ -322,6 +325,7 @@ subroutine step_therm2 (dt) call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_query_parameters(update_ocn_f_out=update_ocn_f) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -939,7 +943,6 @@ subroutine ocn_mixed_layer_icepack(ice, & frzmlt , & ! freezing/melting potential (W/m^2) fhocn_tot , & ! net total heat flux to ocean (W/m^2) fresh_tot ! fresh total water flux to ocean (kg/m^2/s) - real (kind=dbl_kind), parameter :: & frzmlt_max = c1000 ! max magnitude of frzmlt (W/m^2) @@ -1190,8 +1193,8 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) ! tendencies needed by fesom !----------------------------------------------------------------- - dhi_dt(:) = vice(:) - dhs_dt(:) = vsno(:) + dhi_t_dt(:) = vice(:) + dhs_t_dt(:) = vsno(:) !----------------------------------------------------------------- ! initialize diagnostics @@ -1213,18 +1216,19 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) call step_therm1 (dt) ! vertical thermodynamics call step_therm2 (dt) ! ice thickness distribution thermo + !----------------------------------------------------------------- + ! clean up, update tendency diagnostics + !----------------------------------------------------------------- + offset = dt + call update_state (dt, daidtt, dvidtt, dagedtt, offset) + !----------------------------------------------------------------- ! tendencies needed by fesom !----------------------------------------------------------------- - dhi_dt(:) = ( vice(:) - dhi_dt(:) ) / dt - dhs_dt(:) = ( vsno(:) - dhs_dt(:) ) / dt + dhi_t_dt(:) = ( vice(:) - dhi_t_dt(:) ) / dt + dhs_t_dt(:) = ( vsno(:) - dhs_t_dt(:) ) / dt - ! clean up, update tendency diagnostics - - offset = dt - call update_state (dt, daidtt, dvidtt, dagedtt, offset) - !----------------------------------------------------------------- ! dynamics, transport, ridging !----------------------------------------------------------------- @@ -1282,19 +1286,39 @@ module subroutine step_icepack(ice, mesh, time_evp, time_advec, time_therm) t3 = MPI_Wtime() time_advec = t3 - t2 - + + !----------------------------------------------------------------- + ! initialize tendencies needed by fesom + !----------------------------------------------------------------- + dhi_r_dt(:) = vice(:) + dhs_r_dt(:) = vsno(:) + !----------------------------------------------------------------- ! ridging !----------------------------------------------------------------- - call step_dyn_ridge (dt_dyn, ndtd) ! clean up, update tendency diagnostics offset = c0 call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset) - + + !----------------------------------------------------------------- + ! tendencies needed by fesom + !----------------------------------------------------------------- + ! --> ridging adds fresh water need to know for compensation of thermodynamic + ! growth rate of ice and snow in fesom + dhi_r_dt(:) = ( vice(:) - dhi_r_dt(:) ) / dt + dhs_r_dt(:) = ( vsno(:) - dhs_r_dt(:) ) / dt + enddo - + + !----------------------------------------------------------------- + ! total tendencies of thermodynamic and ridging needed by fesom + !----------------------------------------------------------------- + ! --> needed for total compenstion of fresh of thgr and thgrsn + dhi_dt(:) = dhi_r_dt(:) + dhi_t_dt(:) + dhs_dt(:) = dhs_r_dt(:) + dhs_t_dt(:) + !----------------------------------------------------------------- ! albedo, shortwave radiation !----------------------------------------------------------------- diff --git a/src/icepack_drivers/icedrv_transfer.F90 b/src/icepack_drivers/icedrv_transfer.F90 index 524b9d9ec..174b48bfc 100644 --- a/src/icepack_drivers/icedrv_transfer.F90 +++ b/src/icepack_drivers/icedrv_transfer.F90 @@ -50,12 +50,13 @@ module subroutine fesom_to_icepack(ice, mesh) frcidf = 0.17_dbl_kind, & ! frac of incoming sw in near IR diffuse band R_dry = 287.05_dbl_kind, & ! specific gas constant for dry air (J/K/kg) R_vap = 461.495_dbl_kind, & ! specific gas constant for water vapo (J/K/kg) - rhowat = 1025.0_dbl_kind, & ! Water density - cc = rhowat*4190.0_dbl_kind, & ! Volumetr. heat cap. of water [J/m**3/K](cc = rhowat*cp_water) ex = 0.286_dbl_kind - - integer(kind=dbl_kind) :: i, n, k, elem - real (kind=int_kind) :: tx, ty, tvol + !PS rhowat = 1025.0_dbl_kind, & ! Water density + !PS cc = rhowat*4190.0_dbl_kind, & ! Volumetr. heat cap. of water [J/m**3/K](cc = rhowat*cp_water) + + integer(kind=int_kind) :: i, n, k, elem + real (kind=dbl_kind) :: tx, ty, tvol + real (kind=dbl_kind) :: rhowat real (kind=dbl_kind) :: & aux, & @@ -66,6 +67,11 @@ module subroutine fesom_to_icepack(ice, mesh) real(kind=WP), dimension(:), pointer :: u_ice, v_ice, S_oc_array, T_oc_array, & u_w, v_w, stress_atmice_x, stress_atmice_y #include "associate_mesh.h" + + ! make sure we use the same water density to scale precip and snow in + ! fesom and icepack + call icepack_query_parameters(rhow_out=rhowat) + u_ice => ice%uice(:) v_ice => ice%vice(:) S_oc_array => ice%srfoce_salt(:) @@ -87,8 +93,8 @@ module subroutine fesom_to_icepack(ice, mesh) vatm(:) = v_wind(:) fsw(:) = shortwave(:) flw(:) = longwave(:) - frain(:) = prec_rain(:) * 1000.0_dbl_kind - fsnow(:) = prec_snow(:) * 1000.0_dbl_kind + frain(:) = prec_rain(:) * rhowat + fsnow(:) = prec_snow(:) * rhowat wind(:) = sqrt(uatm(:)**2 + vatm(:)**2) @@ -224,7 +230,7 @@ module subroutine icepack_to_fesom( nx_in, & if (present(dhs_dt_out) ) dhs_dt_out = dhs_dt if (present(fsalt_out) ) fsalt_out = fsalt if (present(evap_ocn_out) ) evap_ocn_out = evap_ocn - if (present(evap_ocn_out) ) evap_out = evap + if (present(evap_out) ) evap_out = evap end subroutine icepack_to_fesom From 15a2ba4b6ebad1eaf2c3af78aae787c01ce50cdf Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 10 Oct 2024 14:01:11 +0200 Subject: [PATCH 06/18] make last changes so that icepack is fully volume conserving under zstar, there are also two changes in icepack itself --- src/ice_oce_coupling.F90 | 18 +++++++++++------- src/icepack_drivers/icedrv_step.F90 | 21 ++++++++++++++++++++- 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index d9edf1eec..ac8ed2e1a 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -376,8 +376,10 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) ! Evaporation (convert units from icepack to fesom) evaporation = - evaporation * (1.0_WP - a_ice) * inv_rhowat -!PS ! Ice-Sublimation is not added to evap in icepack -!PS ice_sublimation = 0.0_WP + ! Ice-Sublimation is added to to the freshwater in icepack --> see + ! icepack_therm_vertical.90 --> subroutine thermo_vertical(...): Line: 453 + ! freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt , evapn==sublimation + ice_sublimation = - ice_sublimation * inv_rhowat call init_flux_atm_ocn() @@ -531,13 +533,12 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) !$OMP PARALLEL DO do n=1, myDim_nod2D+eDim_nod2D - flux(n) = evaporation(n) & - +prec_rain(n) & - + flux(n) = evaporation(n) & + -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean + +prec_rain(n) & #if defined (__icepack) +prec_snow(n)*(1.0_WP-a_ice(n)) & #else - -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean +prec_snow(n)*(1.0_WP-a_ice_old(n)) & #endif @@ -623,6 +624,8 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) !$OMP END PARALLEL DO end if + + !___________________________________________________________________________ !---fwf-code-begin if(use_landice_water) then !$OMP PARALLEL DO @@ -631,7 +634,8 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) end do !$OMP END PARALLEL DO end if - + + !___________________________________________________________________________ if(lwiso .and. use_landice_water) then !$OMP PARALLEL DO do n=1, myDim_nod2D+eDim_nod2D diff --git a/src/icepack_drivers/icedrv_step.F90 b/src/icepack_drivers/icedrv_step.F90 index 826134f36..df527cd86 100644 --- a/src/icepack_drivers/icedrv_step.F90 +++ b/src/icepack_drivers/icedrv_step.F90 @@ -997,7 +997,26 @@ subroutine ocn_mixed_layer_icepack(ice, & fresh = fresh - lfs_corr * ice_ref_salinity / sss endif - fresh_tot = fresh + (-evap_ocn + frain + fsnow)*(c1-aice) + !PS fresh_tot = fresh + (-evap_ocn + frain + fsnow)*(c1-aice) + fresh_tot = fresh + frain + (-evap_ocn + fsnow)*(c1-aice) + ! | | | |--> depending on ice concentration eiter snow adds to the + ! | | | freshwater in the ocean or accumulates on the ice as snoe layer + ! | | | + ! | | |--> evaporation ocean-->atmosphere + ! | | + ! | |--> add here the total rain, at the end all the rain + ! | drains through the ice, therefor comment the line + ! | in ice_pack_therm_itd.F90, subroutine icepack_step_therm2() + ! | Line: 1999 + ! | !!! If i dont do this here im not able to balance + ! | the ocean volume under zstar to maschine precision !!! + ! | + ! |--> at that point fresh contains the freshwater flux contributions + ! from the thermodynamic growth rates of ice and snow but also + ! the contributions from the sublimation of ice-->atmos (see. + ! icepack_therm_vertical.F90-->subroutine thermo_vertical(...) + ! Line: 453 --> freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt + ! evapn == evaporative water flux (kg/m^2/s) from sublimation end subroutine ocn_mixed_layer_icepack From 2efe5b02560331b6d15f77b7244f31de2fde67cf Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 11 Oct 2024 11:24:43 +0200 Subject: [PATCH 07/18] add openmp loop in ice_oce_:coupling.F90, delete uneccessary comment in ../src/icepack_drivers/icedrv_transfer.F90 --- src/ice_oce_coupling.F90 | 24 ++++++++++++++---------- src/icepack_drivers/icedrv_transfer.F90 | 6 ------ 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index ac8ed2e1a..f61f96742 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -367,19 +367,23 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) evap_ocn_out=evaporation, & evap_out=ice_sublimation ) - ! Heat flux - heat_flux = - net_heat_flux +!$OMP PARALLEL DO + do n=1, myDim_nod2d+eDim_nod2d + ! Heat flux + heat_flux(n) = - net_heat_flux(n) - ! Freshwater flux (convert units from icepack to fesom) - water_flux = - (fresh_wa_flux * inv_rhowat) - runoff(:) + ! Freshwater flux (convert units from icepack to fesom) + water_flux(n) = - (fresh_wa_flux(n) * inv_rhowat) - runoff(n) - ! Evaporation (convert units from icepack to fesom) - evaporation = - evaporation * (1.0_WP - a_ice) * inv_rhowat + ! Evaporation (convert units from icepack to fesom) + evaporation(n) = - evaporation(n) * (1.0_WP - a_ice(n)) * inv_rhowat - ! Ice-Sublimation is added to to the freshwater in icepack --> see - ! icepack_therm_vertical.90 --> subroutine thermo_vertical(...): Line: 453 - ! freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt , evapn==sublimation - ice_sublimation = - ice_sublimation * inv_rhowat + ! Ice-Sublimation is added to to the freshwater in icepack --> see + ! icepack_therm_vertical.90 --> subroutine thermo_vertical(...): Line: 453 + ! freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt , evapn==sublimation + ice_sublimation(n) = - ice_sublimation(n) * inv_rhowat + end do +!$OMP END PARALLEL DO call init_flux_atm_ocn() diff --git a/src/icepack_drivers/icedrv_transfer.F90 b/src/icepack_drivers/icedrv_transfer.F90 index 174b48bfc..de317f285 100644 --- a/src/icepack_drivers/icedrv_transfer.F90 +++ b/src/icepack_drivers/icedrv_transfer.F90 @@ -18,12 +18,6 @@ module subroutine fesom_to_icepack(ice, mesh) use g_forcing_param, only: ncar_bulk_z_wind, ncar_bulk_z_tair, & ncar_bulk_z_shum use g_sbf, only: l_mslp -! use i_arrays, only: S_oc_array, T_oc_array, & ! Ocean and sea ice fields -! u_w, v_w, & -! stress_atmice_x, stress_atmice_y -! ! u_ice, v_ice, & - -! use i_param, only: cd_oce_ice ! Sea ice parameters use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_parameters use icepack_intfc, only: icepack_sea_freezing_temperature From fcf69bbffe8ab3623ff824d64e379f0d14583807 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Sun, 13 Oct 2024 10:50:27 +0200 Subject: [PATCH 08/18] fix silly misstake... --- src/icb_coupling.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icb_coupling.F90 b/src/icb_coupling.F90 index 35b29ab9f..7f8f6b7e1 100644 --- a/src/icb_coupling.F90 +++ b/src/icb_coupling.F90 @@ -121,7 +121,7 @@ subroutine prepare_icb2fesom(mesh, partit, ib,i_have_element,localelement,depth_ end if dz = abs( lev_low - lev_up ) if( (abs(lev_low)>=abs(depth_ib)) .and. (abs(lev_up) Date: Sun, 13 Oct 2024 10:50:45 +0200 Subject: [PATCH 09/18] write out 3D icb heat flux --- src/io_meandata.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index afc58f87f..5b2fc8876 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -1049,7 +1049,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwbv', 'basal iceberg melting', 'm/s', ibfwbv(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwl', 'lateral iceberg melting', 'm/s', ibfwl(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwe', 'iceberg erosion', 'm/s', ibfwe(:), 1, 'm', i_real4, partit, mesh) - call def_stream(nod2D, myDim_nod2D, 'ibhf', 'heat flux from iceberg melting', 'W/m2', ibhf(:), 1, 'm', i_real4, partit, mesh) + call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) end if !------------------------------------------ !_______________________________________________________________________________ From cd374e76223ca65bfe560ed80743f521c8219e9e Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 15 Oct 2024 16:23:52 +0200 Subject: [PATCH 10/18] fix icepack io_mean variable writing + fix and improve restart writing and reading --- src/icepack_drivers/icedrv_io.F90 | 316 +++++++++++++++------------- src/icepack_drivers/icedrv_main.F90 | 2 +- src/io_fesom_file.F90 | 82 +++++--- src/io_meandata.F90 | 9 +- src/io_restart.F90 | 15 +- src/io_restart_file_group.F90 | 34 +-- 6 files changed, 249 insertions(+), 209 deletions(-) diff --git a/src/icepack_drivers/icedrv_io.F90 b/src/icepack_drivers/icedrv_io.F90 index 209d15cb7..036b7d07b 100644 --- a/src/icepack_drivers/icedrv_io.F90 +++ b/src/icepack_drivers/icedrv_io.F90 @@ -24,26 +24,26 @@ ! define mean IO output of icepack module subroutine ini_mean_icepack_io(mesh) - use mod_mesh - use io_meandata, only: def_stream3D, def_stream2D + use mod_mesh + use io_meandata, only: def_stream3D, def_stream2D - implicit none + implicit none - type(t_mesh), target, intent(in) :: mesh + type(t_mesh), target, intent(in) :: mesh - integer :: i, j, k, & - nt_Tsfc, nt_sice, nt_qice, nt_qsno, & - nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, & - nt_vlvl, nt_iage, nt_FY, nt_aero, & - ktherm, nt_fbri + integer :: i, j, k, & + nt_Tsfc, nt_sice, nt_qice, nt_qsno, & + nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, & + nt_vlvl, nt_iage, nt_FY, nt_aero, & + ktherm, nt_fbri - integer, save :: nm_io_unit = 102 ! unit to open namelist file - integer, save :: nm_icepack_unit = 103 - integer :: iost - character(len=10) :: id_string - character(500) :: longname, trname, units + integer, save :: nm_io_unit = 102 ! unit to open namelist file + integer, save :: nm_icepack_unit = 103 + integer :: iost + character(len=10) :: id_string + character(500) :: longname, trname, units - logical (kind=log_kind) :: & + logical (kind=log_kind) :: & solve_zsal, skl_bgc, z_tracers, & tr_iage, tr_FY, tr_lvl, tr_aero, tr_pond_cesm, & tr_pond_topo, tr_pond_lvl, tr_brine, & @@ -54,224 +54,236 @@ module subroutine ini_mean_icepack_io(mesh) tr_zaero, tr_bgc_Fe, & tr_bgc_hum - integer, save :: io_listsize=0 - type io_entry - character(len=10) :: id ='unknown ' - integer :: freq =0 - character :: unit ='' - integer :: precision =0 - end type + integer, save :: io_listsize=0 + type io_entry + character(len=10) :: id ='unknown ' + integer :: freq =0 + character :: unit ='' + integer :: precision =0 + end type - type(io_entry), save, allocatable, target :: io_list_icepack(:) + type(io_entry), save, allocatable, target :: io_list_icepack(:) - namelist /nml_listsize / io_listsize - namelist /nml_list_icepack / io_list_icepack + namelist /nml_general / io_listsize + namelist /nml_list_icepack / io_list_icepack #include "associate_mesh.h" - ! Get the tracers information from icepack - call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & - nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) - call icepack_query_tracer_indices( & - nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & - nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & - nt_iage_out=nt_iage, nt_FY_out=nt_FY, & - nt_qice_out=nt_qice, nt_sice_out=nt_sice, nt_fbri_out=nt_fbri, & - nt_aero_out=nt_aero, nt_qsno_out=nt_qsno) - call icepack_query_parameters(solve_zsal_out=solve_zsal, & - skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, ktherm_out=ktherm) - call icepack_query_tracer_flags( & - tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, & - tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & - tr_brine_out=tr_brine, tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, & - tr_bgc_Nit_out=tr_bgc_Nit, tr_bgc_Sil_out=tr_bgc_Sil, & - tr_bgc_DMS_out=tr_bgc_DMS, & - tr_bgc_chl_out=tr_bgc_chl, tr_bgc_Am_out=tr_bgc_Am, & - tr_bgc_PON_out=tr_bgc_PON, tr_bgc_DON_out=tr_bgc_DON, & - tr_zaero_out=tr_zaero, tr_bgc_Fe_out=tr_bgc_Fe, & - tr_bgc_hum_out=tr_bgc_hum) + ! Get the tracers information from icepack + call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & + nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) + cAll icepack_query_tracer_indices( & + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, & + nt_qice_out=nt_qice, nt_sice_out=nt_sice, nt_fbri_out=nt_fbri, & + nt_aero_out=nt_aero, nt_qsno_out=nt_qsno) + call icepack_query_parameters(solve_zsal_out=solve_zsal, & + skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, ktherm_out=ktherm) + call icepack_query_tracer_flags( & + tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, & + tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_brine_out=tr_brine, tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, & + tr_bgc_Nit_out=tr_bgc_Nit, tr_bgc_Sil_out=tr_bgc_Sil, & + tr_bgc_DMS_out=tr_bgc_DMS, & + tr_bgc_chl_out=tr_bgc_chl, tr_bgc_Am_out=tr_bgc_Am, & + tr_bgc_PON_out=tr_bgc_PON, tr_bgc_DON_out=tr_bgc_DON, & + tr_zaero_out=tr_zaero, tr_bgc_Fe_out=tr_bgc_Fe, & + tr_bgc_hum_out=tr_bgc_hum) - ! OPEN and read namelist for icepack I/O - open( unit=nm_io_unit, file='namelist.io', form='formatted', access='sequential', status='old', iostat=iost ) - if (iost == 0) then - if (mype==0) write(*,*) ' file : ', 'namelist.io',' open ok' - else - if (mype==0) write(*,*) 'ERROR: --> bad opening file : ','namelist.io',' ; iostat=',iost - call par_ex(p_partit%MPI_COMM_FESOM, p_partit%mype) - stop - end if - open( unit=nm_icepack_unit, file='namelist.icepack', form='formatted', access='sequential', status='old', iostat=iost ) - if (iost == 0) then - if (mype==0) write(*,*) ' file : ', 'namelist.icepack',' open ok' - else - if (mype==0) write(*,*) 'ERROR: --> bad opening file : ','namelist.icepack',' ; iostat=',iost - call par_ex(p_partit%MPI_COMM_FESOM, p_partit%mype) - stop - end if - - read(nm_io_unit, nml=nml_listsize, iostat=iost ) - allocate(io_list_icepack(io_listsize)) - read(nm_icepack_unit, nml=nml_list_icepack, iostat=iost ) - close(nm_icepack_unit) + !_______________________________________________________________________ + ! OPEN and read namelist.io --> need to extract variable io_listsize + open( unit=nm_io_unit, file='namelist.io', form='formatted', access='sequential', status='old', iostat=iost ) + if (iost == 0) then + if (mype==0) write(*,*) ' file : ', 'namelist.io',' open ok' + else + if (mype==0) write(*,*) 'ERROR: --> bad opening file : ','namelist.io',' ; iostat=',iost + call par_ex(p_partit%MPI_COMM_FESOM, p_partit%mype) + stop + end if + + ! read list_size from namelist.io for allocation + read(nm_io_unit, nml=nml_general, iostat=iost ) + close(nm_io_unit) + allocate(io_list_icepack(io_listsize)) + !_______________________________________________________________________ + ! OPEN and read namelist.icepack --> need to extract io_list_icepack + open( unit=nm_icepack_unit, file='namelist.icepack', form='formatted', access='sequential', status='old', iostat=iost ) + if (iost == 0) then + if (mype==0) write(*,*) ' file : ', 'namelist.icepack',' open ok' + else + if (mype==0) write(*,*) 'ERROR: --> bad opening file : ','namelist.icepack',' ; iostat=',iost + call par_ex(p_partit%MPI_COMM_FESOM, p_partit%mype) + stop + end if + + ! read io_list_icepack from namelist to fill up what has been previously + ! allocated --> allocate(io_list_icepack(io_listsize)) + read(nm_icepack_unit, nml=nml_list_icepack, iostat=iost ) + close(nm_icepack_unit) - do i=1, io_listsize - if (trim(io_list_icepack(i)%id)=='unknown ') then + !_______________________________________________________________________ + ! reduce running index to the number that is actually filt up + do i=1, io_listsize + if (trim(io_list_icepack(i)%id)=='unknown ') then if (mype==0) write(*,*) 'io_listsize will be changed from ', io_listsize, ' to ', i-1, '!' io_listsize=i-1 exit - end if - end do + end if + end do - do i=1, io_listsize - select case (trim(io_list_icepack(i)%id)) - case ('aice0 ') - call def_stream2D(nod2D, nx_nh, 'aice0', 'open water fraction', 'none', aice0(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('aicen ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'aicen', 'sea ice concentration', 'none', aicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - case ('vicen ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vicen', 'volume per unit area of ice', 'm', vicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - case ('vsnon ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vsnon', 'volume per unit area of snow', 'm', vsnon(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - case ('aice ') - call def_stream2D(nod2D, nx_nh, 'aice', 'sea ice concentration', 'none', aice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('vice ') - call def_stream2D(nod2D, nx_nh, 'vice', 'volume per unit area of ice', 'm', vice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('vsno ') - call def_stream2D(nod2D, nx_nh, 'vsno', 'volume per unit area of snow', 'm', vsno(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - ! Sea ice velocity components - case ('uvel ') + !_______________________________________________________________________ + ! define output streams + do i=1, io_listsize + select case (trim(io_list_icepack(i)%id)) + case ('aice0 ') + call def_stream2D(nod2D, nx_nh, 'aice0', 'open water fraction', 'none', aice0(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + case ('aicen ') + call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'aicen', 'sea ice concentration', 'none', aicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + case ('vicen ') + call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vicen', 'volume per unit area of ice', 'm', vicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + case ('vsnon ') + call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vsnon', 'volume per unit area of snow', 'm', vsnon(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + case ('aice ') + call def_stream2D(nod2D, nx_nh, 'aice', 'sea ice concentration', 'none', aice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + case ('vice ') + call def_stream2D(nod2D, nx_nh, 'vice', 'volume per unit area of ice', 'm', vice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + case ('vsno ') + call def_stream2D(nod2D, nx_nh, 'vsno', 'volume per unit area of snow', 'm', vsno(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + ! Sea ice velocity components + case ('uvel ') call def_stream2D(nod2D, nx_nh, 'uvel', 'x-component of ice velocity', 'm/s', uvel(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('vvel ') + case ('vvel ') call def_stream2D(nod2D, nx_nh, 'vvel', 'y-component of ice velocity', 'm/s', vvel(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - ! Sea ice or snow surface temperature - case ('Tsfc ') + ! Sea ice or snow surface temperature + case ('Tsfc ') call def_stream2D(nod2D, nx_nh, 'Tsfc', 'sea ice surf. temperature', 'degC', trcr(:,nt_Tsfc), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('Tsfcn ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'Tsfcn', 'sea ice surf. temperature', 'degC', trcrn(:,nt_Tsfc,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('strength ') + case ('Tsfcn ') + call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'Tsfcn', 'sea ice surf. temperature', 'degC', trcrn(:,nt_Tsfc,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + case ('strength ') call def_stream2D(nod2D, nx_nh, 'strength', 'sea ice strength', 'N', strength(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - ! If the following tracers are not defined they will not be outputed - case ('iagen ') + ! If the following tracers are not defined they will not be outputed + case ('iagen ') if (tr_iage) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'iage', 'sea ice age', 's', trcrn(:,nt_iage,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('FYn ') + case ('FYn ') if (tr_FY) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'FY', 'first year ice', 'none', trcrn(:,nt_FY,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('lvln ') + case ('lvln ') if (tr_lvl) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'alvl', 'ridged sea ice area', 'none', trcrn(:,nt_alvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vlvl', 'ridged sea ice volume', 'm', trcrn(:,nt_vlvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('pond_cesmn') + case ('pond_cesmn') if (tr_pond_cesm) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('pond_topon') + case ('pond_topon') if (tr_pond_topo) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd', 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('pond_lvln ') + case ('pond_lvln ') if (tr_pond_lvl) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd', 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('brinen ') + case ('brinen ') if (tr_brine) then call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'fbri', 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if - case ('qicen ') - do k = 1,nilyr ! Separate variable for each sea ice layer + case ('qicen ') + do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A22,i1)') 'sea ice enthalpy lyr: ', k units='J/m3' call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - end do - case ('sicen ') - do k = 1,nilyr ! Separate variable for each sea ice layer + end do + case ('sicen ') + do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A22,i1)') 'sea ice salinity lyr: ', k units='psu' call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - end do - case ('qsnon ') - do k = 1,nslyr ! Separate variable for each snow layer + end do + case ('qsnon ') + do k = 1,nslyr ! Separate variable for each snow layer write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A19,i1)') 'snow enthalpy lyr: ', k units='J/m3' call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - end do - ! Average over categories - case ('iage ') + end do + ! Average over categories + case ('iage ') if (tr_iage) then - call def_stream2D(nod2D, nx_nh, 'iage', 'sea ice age', 's', trcr(:,nt_iage), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'iage', 'sea ice age', 's', trcr(:,nt_iage), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('FY ') + case ('FY ') if (tr_FY) then - call def_stream2D(nod2D, nx_nh, 'FY', 'first year ice', 'none', trcr(:,nt_FY), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'FY', 'first year ice', 'none', trcr(:,nt_FY), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('lvl ') + case ('lvl ') if (tr_lvl) then - call def_stream2D(nod2D, nx_nh, 'alvl', 'ridged sea ice area', 'none', trcr(:,nt_alvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'vlvl', 'ridged sea ice volume', 'm', trcr(:,nt_vlvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'alvl', 'ridged sea ice area', 'none', trcr(:,nt_alvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'vlvl', 'ridged sea ice volume', 'm', trcr(:,nt_vlvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('pond_cesm ') + case ('pond_cesm ') if (tr_pond_cesm) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('pond_topo ') + case ('pond_topo ') if (tr_pond_topo) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('pond_lvl ') + case ('pond_lvl ') if (tr_pond_lvl) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - !call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + !call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('brine ') + case ('brine ') if (tr_brine) then - call def_stream2D(nod2D, nx_nh, 'fbri', 'volume fraction of ice with dynamic salt', 'none', trcr(:,nt_fbri), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream2D(nod2D, nx_nh, 'fbri', 'volume fraction of ice with dynamic salt', 'none', trcr(:,nt_fbri), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if - case ('qice ') - do k = 1,nilyr ! Separate variable for each sea ice layer + case ('qice ') + do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A22,i1)') 'sea ice enthalpy lyr: ', k units='J/m3' call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_qice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - end do - case ('sice ') - do k = 1,nilyr ! Separate variable for each sea ice layer + end do + case ('sice ') + do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A22,i1)') 'sea ice salinity lyr: ', k units='psu' call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_sice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - end do - case ('qsno ') - do k = 1,nslyr ! Separate variable for each snow layer + end do + case ('qsno ') + do k = 1,nslyr ! Separate variable for each snow layer write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A19,i1)') 'snow enthalpy lyr: ', k units='J/m3' call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_qsno+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - end do - case ('rdg_conv ') - call def_stream2D(nod2D, nx_nh, 'rdg_conv', 'Convergence term for ridging', '1/s', rdg_conv(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case ('rdg_shear ') - call def_stream2D(nod2D, nx_nh, 'rdg_shear', 'Shear term for ridging', '1/s', rdg_shear(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - case default - if (mype==0) write(*,*) 'stream ', io_list_icepack(i)%id, ' is not defined !' - end select - end do + end do + case ('rdg_conv ') + call def_stream2D(nod2D, nx_nh, 'rdg_conv', 'Convergence term for ridging', '1/s', rdg_conv(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + case ('rdg_shear ') + call def_stream2D(nod2D, nx_nh, 'rdg_shear', 'Shear term for ridging', '1/s', rdg_shear(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + case default + if (mype==0) write(*,*) 'stream ', io_list_icepack(i)%id, ' is not defined !' + end select + end do ! --> do i=1, io_listsize end subroutine ini_mean_icepack_io diff --git a/src/icepack_drivers/icedrv_main.F90 b/src/icepack_drivers/icedrv_main.F90 index 3c609508c..6428aad72 100644 --- a/src/icepack_drivers/icedrv_main.F90 +++ b/src/icepack_drivers/icedrv_main.F90 @@ -905,7 +905,7 @@ module subroutine ini_icepack_io(year, partit, mesh) type(t_partit), intent(inout), target :: partit integer(kind=int_kind), intent(in) :: year end subroutine ini_icepack_io - + ! Cut off Icepack module subroutine cut_off_icepack use icepack_intfc, only: icepack_compute_tracers diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 0df58c7da..1347cfe2c 100644 --- a/src/io_fesom_file.F90 +++ b/src/io_fesom_file.F90 @@ -15,6 +15,7 @@ module io_fesom_file_module real(kind=8), allocatable :: global_level_data(:) integer :: global_level_data_size = 0 logical is_elem_based + logical :: is_icepackvar2=.false. character(:), allocatable :: varname ! todo: maybe use a getter in netcdf_file_type to get the name end type @@ -54,7 +55,7 @@ module io_fesom_file_module integer, save :: m_nod2d integer, save :: m_elem2d integer, save :: m_nl - integer, save :: m_nitc + integer, save :: m_ncat type fesom_file_type_ptr @@ -104,7 +105,7 @@ function time_dimindex(this) result(x) end function - subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) ! todo: would like to call it initialize but Fortran is rather cluncky with overwriting base type procedures + subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_ncat) ! todo: would like to call it initialize but Fortran is rather cluncky with overwriting base type procedures use io_netcdf_workaround_module use io_gather_module use MOD_PARTIT @@ -112,7 +113,7 @@ subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) ! tod integer mesh_nod2d integer mesh_elem2d integer mesh_nl - integer, optional :: mesh_nitc + integer, optional :: mesh_ncat type(t_partit), target :: partit ! EO parameters type(fesom_file_type_ptr), allocatable :: tmparr(:) @@ -126,11 +127,11 @@ subroutine init(this, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) ! tod m_nod2d = mesh_nod2d m_elem2d = mesh_elem2d m_nl = mesh_nl - !PS mesh_nitc ... icepack number of ice thickness classes, - if (present(mesh_nitc)) then - m_nitc = mesh_nitc + !PS mesh_ncat ... icepack number of ice thickness classes, + if (present(mesh_ncat)) then + m_ncat = mesh_ncat else - m_nitc = 0 + m_ncat = 0 end if call this%netcdf_file_type%initialize() @@ -191,9 +192,18 @@ subroutine read_and_scatter_variables(this) do i=1, this%nvar_infos var => this%var_infos(i) - nlvl = size(var%external_local_data_ptr,dim=1) + + if (var%is_icepackvar2) then + !PS for icepack external_local_data_ptr has still dimension [nod2, ncat] + !PS but usual fesom output would have [nlev, nod2] so we need to change here + !PS dim= parameter to get the proper dimension + nlvl = size(var%external_local_data_ptr,dim=2) + allocate(laux( size(var%external_local_data_ptr,dim=1) )) ! i.e. myDim_elem2D+eDim_elem2D or myDim_nod2D+eDim_nod2D + else + nlvl = size(var%external_local_data_ptr,dim=1) + allocate(laux( size(var%external_local_data_ptr,dim=2) )) ! i.e. myDim_elem2D+eDim_elem2D or myDim_nod2D+eDim_nod2D + end if is_2d = (nlvl == 1) - allocate(laux( size(var%external_local_data_ptr,dim=2) )) ! i.e. myDim_elem2D+eDim_elem2D or myDim_nod2D+eDim_nod2D if(this%is_iorank()) then ! todo: choose how many levels we read at once @@ -233,8 +243,13 @@ subroutine read_and_scatter_variables(this) dynamics_workaround%uv_rhsAB(2,lvl,:) = laux else #endif - ! the data from our pointer is not contiguous (if it is 3D data), so we can not pass the pointer directly to MPI - var%external_local_data_ptr(lvl,:) = laux ! todo: remove this buffer and pass the data directly to MPI (change order of data layout to be levelwise or do not gather levelwise but by columns) + if (var%is_icepackvar2) then + ! the data from our pointer is not contiguous (if it is 3D data), so we can not pass the pointer directly to MPI + var%external_local_data_ptr(:,lvl) = laux ! todo: remove this buffer and pass the data directly to MPI (change order of data layout to be levelwise or do not gather levelwise but by columns) + else + ! the data from our pointer is not contiguous (if it is 3D data), so we can not pass the pointer directly to MPI + var%external_local_data_ptr(lvl,:) = laux ! todo: remove this buffer and pass the data directly to MPI (change order of data layout to be levelwise or do not gather levelwise but by columns) + end if #ifdef ENABLE_NVHPC_WORKAROUNDS end if #endif @@ -410,7 +425,11 @@ subroutine async_gather_and_write_variables(this) var%local_data_copy = dynamics_workaround%uv_rhsAB(2,:,:) else #endif - var%local_data_copy = var%external_local_data_ptr + if (var%is_icepackvar2) then + var%local_data_copy = transpose(var%external_local_data_ptr) + else + var%local_data_copy = var%external_local_data_ptr + end if #ifdef ENABLE_NVHPC_WORKAROUNDS end if #endif @@ -457,21 +476,22 @@ subroutine specify_node_var_2d(this, name, longname, units, local_data) call specify_variable(this, name, [level_diminfo%idx, this%time_dimidx], level_diminfo%len, external_local_data_ptr, .false., longname, units) end subroutine - subroutine specify_node_var_2dicepack(this, name, longname, units, local_data, nitc) + subroutine specify_node_var_2dicepack(this, name, longname, units, local_data, ncat) use, intrinsic :: ISO_C_BINDING class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname real(kind=8), target, intent(inout) :: local_data(:,:) - integer, intent(in) :: nitc! todo: be able to set precision + integer, intent(in) :: ncat! todo: be able to set precision ! EO parameters - type(dim_info) level_diminfo, nitc_diminfo + type(dim_info) level_diminfo, ncat_diminfo -!PS write(*,*) "--> specify_node_var_2dicepack:", __LINE__, __FILE__ +!PS write(*,*) "--> specify_node_var_2dicepack:", __LINE__, __FILE__, size(local_data) level_diminfo = obtain_diminfo(this, m_nod2d) - nitc_diminfo = obtain_diminfo(this, size(local_data, dim=2)) + ncat_diminfo = obtain_diminfo(this, size(local_data, dim=2)) - call specify_variable(this, name, [level_diminfo%idx, nitc_diminfo%idx, this%time_dimidx], level_diminfo%len, local_data, .false., longname, units) + call specify_variable(this, name, [level_diminfo%idx, ncat_diminfo%idx, this%time_dimidx], level_diminfo%len, local_data, .false., longname, units, ncat) + end subroutine @@ -542,8 +562,6 @@ function obtain_diminfo(this, len) result(info) end if end do -!PS write(*,*) 'm_n2d=',m_nod2d,', m_e2d=',m_elem2d,', m_nl=', m_nl,', nitc', m_nitc, ', LEN=', len - ! the dim has not been added yet, see if it is one of our allowed mesh related dims if (len == m_nod2d) then info = dim_info( idx=this%add_dim('node', len), len=len) @@ -553,8 +571,8 @@ function obtain_diminfo(this, len) result(info) info = dim_info( idx=this%add_dim('nz_1', len), len=len) else if(len == m_nl ) then info = dim_info( idx=this%add_dim('nz' , len), len=len) - else if(len == m_nitc ) then - info = dim_info( idx=this%add_dim('nitc', len), len=len) + else if(len == m_ncat ) then + info = dim_info( idx=this%add_dim('ncat', len), len=len) else print *, "error in line ",__LINE__, __FILE__," can not find dimension with size",len stop 1 @@ -570,14 +588,16 @@ function obtain_diminfo(this, len) result(info) end function - subroutine specify_variable(this, name, dim_indices, global_level_data_size, local_data, is_elem_based, longname, units) - type(fesom_file_type), intent(inout) :: this - character(len=*), intent(in) :: name - integer, intent(in) :: dim_indices(:) - integer global_level_data_size - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision? - logical, intent(in) :: is_elem_based - character(len=*), intent(in) :: units, longname + subroutine specify_variable(this, name, dim_indices, global_level_data_size, local_data, is_elem_based, longname, units, ncat) + type(fesom_file_type), intent(inout) :: this + character(len=*) , intent(in) :: name + integer , intent(in) :: dim_indices(:) + integer , intent(in) :: global_level_data_size + real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision? + logical , intent(in) :: is_elem_based + character(len=*) , intent(in) :: units, longname + integer , intent(in) , optional :: ncat + ! EO parameters integer var_index @@ -592,6 +612,8 @@ subroutine specify_variable(this, name, dim_indices, global_level_data_size, loc this%var_infos(this%nvar_infos)%global_level_data_size = global_level_data_size this%var_infos(this%nvar_infos)%is_elem_based = is_elem_based this%var_infos(this%nvar_infos)%varname = name + if (present(ncat)) this%var_infos(this%nvar_infos)%is_icepackvar2=.true. + end subroutine diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index cc08d77e4..afa5e0a3f 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -400,7 +400,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'dflux', 'density flux', 'kg/(m3*s)', dens_flux(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) CASE ('runoff ') sel_forcvar(10)= 1 - call def_stream(nod2D, myDim_nod2D, 'runoff', 'river runoff', 'm/s', runoff(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(nod2D, myDim_nod2D, 'runoff', 'river runoff', 'm/s', runoff(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) CASE ('evap ') sel_forcvar(7) = 1 call def_stream(nod2D, myDim_nod2D, 'evap', 'evaporation', 'm/s', evaporation(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) @@ -1770,13 +1770,14 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) ! define output streams-->dimension, variable, long_name, units, array, freq, unit, precision !PS if (partit%flag_debug .and. partit%mype==0) print *, achar(27)//'[32m'//' -I/O-> call ini_mean_io'//achar(27)//'[0m' call ini_mean_io(ice, dynamics, tracers, partit, mesh) + +#if defined (__icepack) + call ini_mean_icepack_io(mesh) !icapack has its copy of p_partit => partit +#endif !PS if (partit%flag_debug .and. partit%mype==0) print *, achar(27)//'[33m'//' -I/O-> call init_io_gather'//achar(27)//'[0m' call init_io_gather(partit) -#if defined (__icepack) - call ini_mean_icepack_io(mesh) !icapack has its copy of p_partit => partit -#endif end if ! --> if (lfirst) then !___________________________________________________________________________ diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 482f098a3..f076a6d8e 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -38,7 +38,6 @@ MODULE io_RESTART #if defined(__icepack) type(restart_file_group) , save, public :: icepack_files character(:), allocatable, save, public :: icepack_path - #endif #if defined(__recom) @@ -244,11 +243,8 @@ end subroutine ini_bio_io ! subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tracers, partit, mesh) -! #if defined(__icepack) -! icepack restart not merged here ! produce a compiler error if USE_ICEPACK=ON; todo: merge icepack restart from 68d8b8b -! #endif - use fortran_utils + implicit none ! this is the main restart subroutine ! if l_read is TRUE the restart file will be read @@ -273,6 +269,7 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr integer, intent(out):: which_readr integer :: cstep + !_____________________________________________________________________________ ! initialize directory for core dump restart if(.not. initialized_raw) then @@ -316,7 +313,7 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr call ini_ocean_io(yearnew, dynamics, tracers, partit, mesh) if (use_ice) then call ini_ice_io(yearnew, ice, partit, mesh) -#if defined(__icepack) +#if defined(__icepack) call ini_icepack_io(yearnew, partit, mesh) #endif end if @@ -329,7 +326,7 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr if (use_ice) then call ini_ice_io (yearold, ice, partit, mesh) #if defined(__icepack) - call ini_icepack_io(yearnew, partit, mesh) + call ini_icepack_io(yearold, partit, mesh) #endif end if #if defined(__recom) @@ -387,6 +384,10 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr if (use_ice) then if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ice'//achar(27)//'[0m' call read_restart(ice_path, ice_files, partit%MPI_COMM_FESOM, partit%mype) +#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_restart(icepack_path, icepack_files, partit%MPI_COMM_FESOM, partit%mype) +#endif end if #if defined(__recom) diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index bbd139805..ac4b12e9f 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -51,20 +51,20 @@ subroutine def_node_var_2d(this, name, longname, units, local_data, mesh, partit !PS add a seprate case for icepack since here i need to write a 2d vertice file !PS over the dimension of number of ice thickness classes. Additional input - !PS parameter nitc - subroutine def_node_var_2dicepack(this, name, longname, units, local_data, nitc, mesh, partit) + !PS parameter ncat + subroutine def_node_var_2dicepack(this, name, longname, units, local_data, ncat, mesh, partit) use mod_mesh class(restart_file_group), target, intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - integer :: nitc + integer :: ncat type(t_mesh) , intent(in) :: mesh type(t_partit) , intent(in) :: partit ! EO parameters !PS write(*,*) "--> def_node_var_2dicepack:", __LINE__, __FILE__ - call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, nitc) + call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) end subroutine subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit) @@ -112,16 +112,16 @@ subroutine def_elem_var_3d(this, name, longname, units, local_data, mesh, partit end subroutine - subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) + subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_ncat) class(restart_file_group), target, intent(inout) :: g character(len=*), intent(in) :: name logical must_exist_on_read integer mesh_nod2d, mesh_elem2d, mesh_nl type(t_partit), intent(in) :: partit - !PS mesh_nitc ... icepack number of ice thickness classes, do it here as optional + !PS mesh_ncat ... icepack number of ice thickness classes, do it here as optional !PS parameter, i assume is the less intrusive option. Therefore it must be !PS at the end of the argumnent input list - integer, optional :: mesh_nitc + integer, optional :: mesh_ncat ! EO parameters type(restart_file_type), pointer :: f @@ -133,8 +133,12 @@ subroutine add_file(g, name, must_exist_on_read, mesh_nod2d, mesh_elem2d, mesh_n allocate(f%varname,source=name) f%must_exist_on_read = must_exist_on_read - if (present(mesh_nitc)) then - call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_nitc) + if (present(mesh_ncat)) then + call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit, mesh_ncat) + ! make sure we can identify a 2d icepackvariable, problem here is, We have a new + ! dimension ncat number of ice thickness categories and icepack 2d vars have + ! dimensions (nod2, ncat) but fesom IO expects (ncat, nod2) --> means we + ! need to sneek in a matrix transpose somehow!!! else call f%fesom_file_type%init(mesh_nod2d, mesh_elem2d, mesh_nl, partit) end if @@ -159,20 +163,20 @@ subroutine def_node_var_2d_optional(this, name, longname, units, local_data, mes !PS add a seprate case for icepack since here i need to write a 2d vertice file !PS over the dimension of number of ice thickness classes. Additional input - !PS parameter nitc - subroutine def_node_var_2dicepack_optional(this, name, longname, units, local_data, nitc, mesh, partit) + !PS parameter ncat + subroutine def_node_var_2dicepack_optional(this, name, longname, units, local_data, ncat, mesh, partit) use mod_mesh class(restart_file_group), target, intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - integer :: nitc + integer :: ncat type(t_mesh) , intent(in) :: mesh type(t_partit) , intent(in) :: partit ! EO parameters - call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, nitc) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, nitc) + call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) end subroutine subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit) From bfea5f19023300e0984ded968fea1e4febe97a18 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 15 Oct 2024 16:59:57 +0200 Subject: [PATCH 11/18] add some comments --- src/io_fesom_file.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 1347cfe2c..48c8bd395 100644 --- a/src/io_fesom_file.F90 +++ b/src/io_fesom_file.F90 @@ -245,6 +245,8 @@ subroutine read_and_scatter_variables(this) #endif if (var%is_icepackvar2) then ! the data from our pointer is not contiguous (if it is 3D data), so we can not pass the pointer directly to MPI + ! PS Again icepack variable demension is [nod2, ncat] thats why we need + ! PS to switch here the index where we sort in the 2d slices var%external_local_data_ptr(:,lvl) = laux ! todo: remove this buffer and pass the data directly to MPI (change order of data layout to be levelwise or do not gather levelwise but by columns) else ! the data from our pointer is not contiguous (if it is 3D data), so we can not pass the pointer directly to MPI From 6951b113be158ac3b8902e3301bf8055d325d528 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 16 Oct 2024 16:28:01 +0200 Subject: [PATCH 12/18] sort and structure icepack io meandata output routine --- src/icepack_drivers/icedrv_io.F90 | 138 +++++++++++++++--------------- 1 file changed, 69 insertions(+), 69 deletions(-) diff --git a/src/icepack_drivers/icedrv_io.F90 b/src/icepack_drivers/icedrv_io.F90 index 036b7d07b..0d628b570 100644 --- a/src/icepack_drivers/icedrv_io.F90 +++ b/src/icepack_drivers/icedrv_io.F90 @@ -25,7 +25,7 @@ module subroutine ini_mean_icepack_io(mesh) use mod_mesh - use io_meandata, only: def_stream3D, def_stream2D + use io_meandata, only: def_stream implicit none @@ -139,147 +139,147 @@ module subroutine ini_mean_icepack_io(mesh) do i=1, io_listsize select case (trim(io_list_icepack(i)%id)) case ('aice0 ') - call def_stream2D(nod2D, nx_nh, 'aice0', 'open water fraction', 'none', aice0(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'aice0' , 'open water fraction' , 'none', aice0(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('aicen ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'aicen', 'sea ice concentration', 'none', aicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'aicen' , 'sea ice concentration' , 'none', aicen(:,:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) case ('vicen ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vicen', 'volume per unit area of ice', 'm', vicen(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'vicen' , 'volume per unit area of ice' , 'm' , vicen(:,:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) case ('vsnon ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vsnon', 'volume per unit area of snow', 'm', vsnon(:,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'vsnon' , 'volume per unit area of snow', 'm' , vsnon(:,:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) case ('aice ') - call def_stream2D(nod2D, nx_nh, 'aice', 'sea ice concentration', 'none', aice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'aice' , 'sea ice concentration' , 'none', aice(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('vice ') - call def_stream2D(nod2D, nx_nh, 'vice', 'volume per unit area of ice', 'm', vice(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'vice' , 'volume per unit area of ice' , 'm' , vice(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('vsno ') - call def_stream2D(nod2D, nx_nh, 'vsno', 'volume per unit area of snow', 'm', vsno(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'vsno' , 'volume per unit area of snow', 'm' , vsno(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) ! Sea ice velocity components case ('uvel ') - call def_stream2D(nod2D, nx_nh, 'uvel', 'x-component of ice velocity', 'm/s', uvel(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'uvel' , 'x-component of ice velocity' , 'm/s' , uvel(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('vvel ') - call def_stream2D(nod2D, nx_nh, 'vvel', 'y-component of ice velocity', 'm/s', vvel(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'vvel' , 'y-component of ice velocity' , 'm/s' , vvel(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) ! Sea ice or snow surface temperature case ('Tsfc ') - call def_stream2D(nod2D, nx_nh, 'Tsfc', 'sea ice surf. temperature', 'degC', trcr(:,nt_Tsfc), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'Tsfc' , 'sea ice surf. temperature' , 'degC', trcr(:,nt_Tsfc) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('Tsfcn ') - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'Tsfcn', 'sea ice surf. temperature', 'degC', trcrn(:,nt_Tsfc,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - case ('strength ') - call def_stream2D(nod2D, nx_nh, 'strength', 'sea ice strength', 'N', strength(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'Tsfcn' , 'sea ice surf. temperature' , 'degC', trcrn(:,nt_Tsfc,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + case ('strength ') + call def_stream(nod2D , nx_nh , 'strength', 'sea ice strength' , 'N' , strength(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) ! If the following tracers are not defined they will not be outputed case ('iagen ') if (tr_iage) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'iage', 'sea ice age', 's', trcrn(:,nt_iage,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'iage' , 'sea ice age' , 's' , trcrn(:,nt_iage,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('FYn ') if (tr_FY) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'FY', 'first year ice', 'none', trcrn(:,nt_FY,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'FY' , 'first year ice' , 'none', trcrn(:,nt_FY,:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('lvln ') if (tr_lvl) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'alvl', 'ridged sea ice area', 'none', trcrn(:,nt_alvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'vlvl', 'ridged sea ice volume', 'm', trcrn(:,nt_vlvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'alvl' , 'ridged sea ice area' , 'none', trcrn(:,nt_alvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'vlvl' , 'ridged sea ice volume' , 'm' , trcrn(:,nt_vlvl,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('pond_cesmn') if (tr_pond_cesm) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('pond_topon') if (tr_pond_topo) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd', 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd' , 'melt pond refrozen lid thickness', 'm',trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('pond_lvln ') if (tr_pond_lvl) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd', 'melt pond area fraction', 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd', 'melt pond depth', 'm', trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd', 'melt pond refrozen lid thickness', 'm', trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'ipnd' , 'melt pond refrozen lid thickness','m', trcrn(:,nt_ipnd,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('brinen ') if (tr_brine) then - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), 'fbri', 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), 'fbri' , 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end if case ('qicen ') do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A22,i1)') 'sea ice enthalpy lyr: ', k units='J/m3' - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end do case ('sicen ') do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A22,i1)') 'sea ice salinity lyr: ', k units='psu' - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end do case ('qsnon ') do k = 1,nslyr ! Separate variable for each snow layer write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A19,i1)') 'snow enthalpy lyr: ', k units='J/m3' - call def_stream3D((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) + call def_stream((/ncat, nod2D/), (/ncat, nx_nh/), trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh, .true.) end do ! Average over categories case ('iage ') if (tr_iage) then - call def_stream2D(nod2D, nx_nh, 'iage', 'sea ice age', 's', trcr(:,nt_iage), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'iage' , 'sea ice age' , 's' , trcr(:,nt_iage) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('FY ') if (tr_FY) then - call def_stream2D(nod2D, nx_nh, 'FY', 'first year ice', 'none', trcr(:,nt_FY), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'FY' , 'first year ice' , 'none', trcr(:,nt_FY) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('lvl ') if (tr_lvl) then - call def_stream2D(nod2D, nx_nh, 'alvl', 'ridged sea ice area', 'none', trcr(:,nt_alvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'vlvl', 'ridged sea ice volume', 'm', trcr(:,nt_vlvl), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'alvl' , 'ridged sea ice area' , 'none', trcr(:,nt_alvl) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'vlvl' , 'ridged sea ice volume' , 'm' , trcr(:,nt_vlvl) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('pond_cesm ') if (tr_pond_cesm) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'apnd' , 'melt pond area fraction' , 'none', trcr(:,nt_apnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'hpnd' , 'melt pond depth' , 'm' , trcr(:,nt_hpnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('pond_topo ') if (tr_pond_topo) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'apnd' , 'melt pond area fraction' , 'none', trcr(:,nt_apnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'hpnd' , 'melt pond depth' , 'm' , trcr(:,nt_hpnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'ipnd' , 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('pond_lvl ') if (tr_pond_lvl) then - call def_stream2D(nod2D, nx_nh, 'apnd', 'melt pond area fraction', 'none', trcr(:,nt_apnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - call def_stream2D(nod2D, nx_nh, 'hpnd', 'melt pond depth', 'm', trcr(:,nt_hpnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) - !call def_stream2D(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'apnd' , 'melt pond area fraction' , 'none', trcr(:,nt_apnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'hpnd' , 'melt pond depth' , 'm' , trcr(:,nt_hpnd) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + !call def_stream(nod2D, nx_nh, 'ipnd', 'melt pond refrozen lid thickness', 'm', trcr(:,nt_ipnd), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('brine ') if (tr_brine) then - call def_stream2D(nod2D, nx_nh, 'fbri', 'volume fraction of ice with dynamic salt', 'none', trcr(:,nt_fbri), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'fbri' , 'volume fraction of ice with dynamic salt', 'none', trcr(:,nt_fbri), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end if case ('qice ') do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A22,i1)') 'sea ice enthalpy lyr: ', k units='J/m3' - call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_qice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , trim(trname), trim(longname) , trim(units), trcr(:,nt_qice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end do case ('sice ') do k = 1,nilyr ! Separate variable for each sea ice layer write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A22,i1)') 'sea ice salinity lyr: ', k units='psu' - call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_sice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , trim(trname), trim(longname) , trim(units), trcr(:,nt_sice+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end do case ('qsno ') do k = 1,nslyr ! Separate variable for each snow layer write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A19,i1)') 'snow enthalpy lyr: ', k units='J/m3' - call def_stream2D(nod2D, nx_nh, trim(trname), trim(longname), trim(units), trcr(:,nt_qsno+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , trim(trname), trim(longname) , trim(units), trcr(:,nt_qsno+k-1), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) end do case ('rdg_conv ') - call def_stream2D(nod2D, nx_nh, 'rdg_conv', 'Convergence term for ridging', '1/s', rdg_conv(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'rdg_conv' , 'Convergence term for ridging', '1/s', rdg_conv(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case ('rdg_shear ') - call def_stream2D(nod2D, nx_nh, 'rdg_shear', 'Shear term for ridging', '1/s', rdg_shear(:), io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) + call def_stream(nod2D , nx_nh , 'rdg_shear', 'Shear term for ridging' , '1/s' , rdg_shear(:) , io_list_icepack(i)%freq, io_list_icepack(i)%unit, io_list_icepack(i)%precision, p_partit, mesh) case default if (mype==0) write(*,*) 'stream ', io_list_icepack(i)%id, ' is not defined !' end select @@ -369,46 +369,46 @@ module subroutine ini_icepack_io(year, partit, mesh) !----------------------------------------------------------------- ! 3D restart fields (ncat) !----------------------------------------------------------------- - call icepack_files%def_node_var('aicen' , 'sea ice concentration' , 'none', aicen(:,:) , ncat, mesh, partit) - call icepack_files%def_node_var('vicen' , 'volum per unit area of ice' , 'm' , vicen(:,:) , ncat, mesh, partit) - call icepack_files%def_node_var('vsnon' , 'volum per unit area of snow' , 'm' , vsnon(:,:) , ncat, mesh, partit) - call icepack_files%def_node_var('Tsfc' , 'sea ice surf. temperature' , 'degC', trcrn(:,nt_Tsfc,:), ncat, mesh, partit) + call icepack_files%def_node_var('aicen' , 'sea ice concentration' , 'none', aicen(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('vicen' , 'volum per unit area of ice' , 'm' , vicen(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('vsnon' , 'volum per unit area of snow' , 'm' , vsnon(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('Tsfc' , 'sea ice surf. temperature' , 'degC', trcrn(:,nt_Tsfc,:), mesh, partit, ncat) if (tr_iage) then - call icepack_files%def_node_var('iage' , 'sea ice age' , 's' , trcrn(:,nt_iage,:), ncat, mesh, partit) + call icepack_files%def_node_var('iage' , 'sea ice age' , 's' , trcrn(:,nt_iage,:), mesh, partit, ncat) end if if (tr_FY) then - call icepack_files%def_node_var('FY' , 'first year ice' , 'none', trcrn(:,nt_FY,:) , ncat, mesh, partit) + call icepack_files%def_node_var('FY' , 'first year ice' , 'none', trcrn(:,nt_FY,:) , mesh, partit, ncat) end if if (tr_lvl) then - call icepack_files%def_node_var('alvl' , 'ridged sea ice area' , 'none', trcrn(:,nt_alvl,:), ncat, mesh, partit) - call icepack_files%def_node_var('vlvl' , 'ridged sea ice volume' , 'm' , trcrn(:,nt_vlvl,:), ncat, mesh, partit) + call icepack_files%def_node_var('alvl' , 'ridged sea ice area' , 'none', trcrn(:,nt_alvl,:), mesh, partit, ncat) + call icepack_files%def_node_var('vlvl' , 'ridged sea ice volume' , 'm' , trcrn(:,nt_vlvl,:), mesh, partit, ncat) end if if (tr_pond_cesm) then - call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) - call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), mesh, partit, ncat) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), mesh, partit, ncat) end if if (tr_pond_topo) then - call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) - call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) - call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), mesh, partit, ncat) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), mesh, partit, ncat) + call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), mesh, partit, ncat) end if if (tr_pond_lvl) then - call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), ncat, mesh, partit) - call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), ncat, mesh, partit) - call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), ncat, mesh, partit) + call icepack_files%def_node_var('apnd' , 'melt pond area fraction' , 'none', trcrn(:,nt_apnd,:), mesh, partit, ncat) + call icepack_files%def_node_var('hpnd' , 'melt pond depth' , 'm' , trcrn(:,nt_hpnd,:), mesh, partit, ncat) + call icepack_files%def_node_var('ipnd' , 'melt pond refrozen lid thickness' , 'm' , trcrn(:,nt_ipnd,:), mesh, partit, ncat) call icepack_files%def_node_var('ffracn', 'fraction of fsurfn over pond used to melt ipond' , 'none', ffracn, mesh, partit); call icepack_files%def_node_var('dhsn' , 'depth difference for snow on sea ice and pond ice' , 'm' , dhsn , mesh, partit); end if if (tr_brine) then - call icepack_files%def_node_var('fbri' , 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), ncat, mesh, partit) - call icepack_files%def_node_var('first_ice', 'distinguishes ice that disappears' , 'logical', first_ice_real(:,:), ncat, mesh, partit) + call icepack_files%def_node_var('fbri' , 'volume fraction of ice with dynamic salt', 'none', trcrn(:,nt_fbri,:), mesh, partit, ncat) + call icepack_files%def_node_var('first_ice', 'distinguishes ice that disappears' , 'logical', first_ice_real(:,:), mesh, partit, ncat) end if !----------------------------------------------------------------- @@ -421,11 +421,11 @@ module subroutine ini_icepack_io(year, partit, mesh) write(trname,'(A6,i1)') 'sicen_', k write(longname,'(A21,i1)') 'sea ice salinity lyr:', k units='psu' - call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), ncat, mesh, partit) + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_sice+k-1,:), mesh, partit, ncat) write(trname,'(A6,i1)') 'qicen_', k write(longname,'(A21,i1)') 'sea ice enthalpy lyr:', k units='J/m3' - call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), ncat, mesh, partit) + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qice+k-1,:), mesh, partit, ncat) end do ! Snow @@ -434,7 +434,7 @@ module subroutine ini_icepack_io(year, partit, mesh) write(trname,'(A6,i1)') 'qsnon_', k write(longname,'(A18,i1)') 'snow enthalpy lyr:', k units='J/m3' - call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), ncat, mesh, partit) + call icepack_files%def_node_var(trim(trname), trim(longname), trim(units), trcrn(:,nt_qsno+k-1,:), mesh, partit, ncat) end do ! From 1de84cdb0852efce91f1f03362df519d9b919007 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 16 Oct 2024 16:29:24 +0200 Subject: [PATCH 13/18] make sure the io mean routine is able to write an array for the ncat dimension properly --- src/io_meandata.F90 | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index afa5e0a3f..7e70534af 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -15,7 +15,7 @@ module io_MEANDATA implicit none #include "netcdf.inc" private - public :: def_stream2D, def_stream3D, output, finalize_output + public :: def_stream, def_stream2D, def_stream3D, output, finalize_output ! !-------------------------------------------------------------------------------------------- ! @@ -1426,7 +1426,12 @@ subroutine create_new_file(entry, ice, dynamics, partit, mesh) type(t_ice) , intent(in) :: ice type(Meandata), intent(inout) :: entry - character(len=*), parameter :: global_attributes_prefix = "FESOM_" + character(len=*), parameter :: global_attributes_prefix = "FESOM_" +#if defined(__icepack) + integer, allocatable :: ncat_arr(:) + integer :: ii +#endif + ! Serial output implemented so far if (partit%mype/=entry%root_rank) return ! create an ocean output file @@ -1443,15 +1448,21 @@ subroutine create_new_file(entry, ice, dynamics, partit, mesh) else if (entry%ndim==2) then call assert_nf( nf_def_dim(entry%ncid, entry%dimname(1), entry%glsize(1), entry%dimID(1)), __LINE__) - call assert_nf( nf_def_var(entry%ncid, entry%dimname(1), NF_DOUBLE, 1, entry%dimID(1), entry%dimvarID(1)), __LINE__) if (entry%dimname(1)=='nz') then + call assert_nf( nf_def_var(entry%ncid, entry%dimname(1), NF_DOUBLE, 1, entry%dimID(1), entry%dimvarID(1)), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%dimvarID(1), 'long_name', len_trim('depth at layer interface'),'depth at layer interface'), __LINE__) + elseif (entry%dimname(1)=='nz1') then + call assert_nf( nf_def_var(entry%ncid, entry%dimname(1), NF_DOUBLE, 1, entry%dimID(1), entry%dimvarID(1)), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%dimvarID(1), 'long_name', len_trim('depth at layer midpoint'),'depth at layer midpoint'), __LINE__) + elseif (entry%dimname(1)=='ncat') then + call assert_nf( nf_def_var(entry%ncid, entry%dimname(1), NF_INT , 1, entry%dimID(1), entry%dimvarID(1)), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%dimvarID(1), 'long_name', len_trim('sea-ice thickness class'),'sea-ice thickness class'), __LINE__) + else if (partit%mype==0) write(*,*) 'WARNING: unknown first dimension in 2d mean I/O data' + end if call assert_nf( nf_put_att_text(entry%ncid, entry%dimvarID(1), 'units', len_trim('m'),'m'), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%dimvarID(1), 'positive', len_trim('down'),'down'), __LINE__) @@ -1519,6 +1530,15 @@ subroutine create_new_file(entry, ice, dynamics, partit, mesh) call assert_nf( nf_put_var_double(entry%ncid, entry%dimvarID(1), abs(mesh%zbar)), __LINE__) elseif (entry%dimname(1)=='nz1') then call assert_nf( nf_put_var_double(entry%ncid, entry%dimvarID(1), abs(mesh%Z)), __LINE__) +#if defined(__icepack) + elseif (entry%dimname(1)=='ncat') then + allocate(ncat_arr(entry%glsize(1))) + do ii= 1, entry%glsize(1) + ncat_arr(ii)=ii + end do + call assert_nf( nf_put_var_int(entry%ncid, entry%dimvarID(1), ncat_arr), __LINE__) + deallocate(ncat_arr) +#endif else if (partit%mype==0) write(*,*) 'WARNING: unknown first dimension in 2d mean I/O data' end if From dbb76f81756775ba27946aee03c7c61c44f47489 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 16 Oct 2024 16:31:40 +0200 Subject: [PATCH 14/18] reduce number of def_nod_var subroutines, make ncat index an optional parameter for def_node_var_3d --- src/io_restart_file_group.F90 | 100 ++++++++++++++-------------------- 1 file changed, 41 insertions(+), 59 deletions(-) diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index ac4b12e9f..5d5445897 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -21,14 +21,14 @@ module restart_file_group_module integer, public :: nfiles = 0 ! todo: allow dynamically allocated size without messing with shallow copied pointers contains - generic, public :: def_node_var => def_node_var_2d, def_node_var_2dicepack, def_node_var_3d + generic, public :: def_node_var => def_node_var_2d, def_node_var_3d generic, public :: def_elem_var => def_elem_var_2d, def_elem_var_3d - procedure, private :: def_node_var_2d, def_node_var_2dicepack, def_node_var_3d + procedure, private :: def_node_var_2d, def_node_var_3d procedure, private :: def_elem_var_2d, def_elem_var_3d ! def_*_optional procedures create a restart variable which does not have to exist when reading the restart file - generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_2dicepack_optional, def_node_var_3d_optional + generic, public :: def_node_var_optional => def_node_var_2d_optional, def_node_var_3d_optional generic, public :: def_elem_var_optional => def_elem_var_2d_optional, def_elem_var_3d_optional - procedure, private :: def_node_var_2d_optional, def_node_var_2dicepack_optional, def_node_var_3d_optional + procedure, private :: def_node_var_2d_optional, def_node_var_3d_optional procedure, private :: def_elem_var_2d_optional, def_elem_var_3d_optional end type @@ -49,36 +49,27 @@ subroutine def_node_var_2d(this, name, longname, units, local_data, mesh, partit call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine - !PS add a seprate case for icepack since here i need to write a 2d vertice file - !PS over the dimension of number of ice thickness classes. Additional input - !PS parameter ncat - subroutine def_node_var_2dicepack(this, name, longname, units, local_data, ncat, mesh, partit) + subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit, ncat) use mod_mesh - class(restart_file_group), target, intent(inout) :: this - character(len=*), intent(in) :: name - character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - integer :: ncat - type(t_mesh) , intent(in) :: mesh - type(t_partit) , intent(in) :: partit - ! EO parameters -!PS write(*,*) "--> def_node_var_2dicepack:", __LINE__, __FILE__ - call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) - end subroutine - - subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit) - use mod_mesh - class(restart_file_group), intent(inout) :: this - character(len=*), intent(in) :: name - character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - type(t_mesh), intent(in) :: mesh - type(t_partit), intent(in) :: partit + class(restart_file_group), intent(inout) :: this + character(len=*) , intent(in) :: name + character(len=*) , intent(in) :: units, longname + real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision + type(t_mesh) , intent(in) :: mesh + type(t_partit) , intent(in) :: partit + integer , intent(in) , optional :: ncat ! EO parameters !PS write(*,*) "--> def_node_var_3d:", __LINE__, __FILE__ - call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + if (present(ncat)) then + !PS add a seprate case for icepack since here i need to write a 2d vertice file + !PS over the dimension of number of ice thickness classes. Additional input + !PS parameter ncat + call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) + else + call add_file(this, name, .true., mesh%nod2d, mesh%elem2d, mesh%nl, partit) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + end if end subroutine @@ -161,36 +152,27 @@ subroutine def_node_var_2d_optional(this, name, longname, units, local_data, mes call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) end subroutine - !PS add a seprate case for icepack since here i need to write a 2d vertice file - !PS over the dimension of number of ice thickness classes. Additional input - !PS parameter ncat - subroutine def_node_var_2dicepack_optional(this, name, longname, units, local_data, ncat, mesh, partit) + subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit, ncat) use mod_mesh - class(restart_file_group), target, intent(inout) :: this - character(len=*), intent(in) :: name - character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - integer :: ncat - type(t_mesh) , intent(in) :: mesh - type(t_partit) , intent(in) :: partit + class(restart_file_group), intent(inout) :: this + character(len=*) , intent(in) :: name + character(len=*) , intent(in) :: units, longname + real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision + type(t_mesh) , intent(in) :: mesh + type(t_partit) , intent(in) :: partit + integer , intent(in) , optional :: ncat ! EO parameters - - call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) - end subroutine - - subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit) - use mod_mesh - class(restart_file_group), intent(inout) :: this - character(len=*), intent(in) :: name - character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision - type(t_mesh), intent(in) :: mesh - type(t_partit), intent(in) :: partit - ! EO parameters - - call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit) - call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + if (present(ncat)) then + !PS add a seprate case for icepack since here i need to write a 2d vertice file + !PS over the dimension of number of ice thickness classes. Additional input + !PS parameter ncat + call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit, ncat) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data, ncat) + else + call add_file(this, name, .false., mesh%nod2d, mesh%elem2d, mesh%nl, partit) + call this%files(this%nfiles)%specify_node_var(name, longname, units, local_data) + end if + end subroutine From 4a5d97de00459f2de12f67a6a2328f922af60515 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Oct 2024 11:55:05 +0200 Subject: [PATCH 15/18] fix the icepack netcdf restart, add additional restart variables, now also the restart is volume conserving. Also only write eiter .ice. or .icepack. restart files --- src/ice_oce_coupling.F90 | 22 ++++++------- src/icepack_drivers/icedrv_io.F90 | 14 +++++--- src/io_restart.F90 | 53 ++++++++++++++++++++----------- 3 files changed, 54 insertions(+), 35 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index f61f96742..373740545 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -355,17 +355,17 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh) ! #if defined (__icepack) - call icepack_to_fesom (nx_in=(myDim_nod2D+eDim_nod2D), & - aice_out=a_ice, & - vice_out=m_ice, & - vsno_out=m_snow, & - fhocn_tot_out=net_heat_flux, & - fresh_tot_out=fresh_wa_flux, & - fsalt_out=real_salt_flux, & - dhs_dt_out=thdgrsn, & - dhi_dt_out=thdgr, & - evap_ocn_out=evaporation, & - evap_out=ice_sublimation ) + call icepack_to_fesom (nx_in = (myDim_nod2D+eDim_nod2D), & + aice_out = a_ice, & + vice_out = m_ice, & + vsno_out = m_snow, & + fhocn_tot_out = net_heat_flux, & + fresh_tot_out = fresh_wa_flux, & + fsalt_out = real_salt_flux, & + dhs_dt_out = thdgrsn, & + dhi_dt_out = thdgr, & + evap_ocn_out = evaporation, & + evap_out = ice_sublimation ) !$OMP PARALLEL DO do n=1, myDim_nod2d+eDim_nod2d diff --git a/src/icepack_drivers/icedrv_io.F90 b/src/icepack_drivers/icedrv_io.F90 index 0d628b570..e14732c05 100644 --- a/src/icepack_drivers/icedrv_io.F90 +++ b/src/icepack_drivers/icedrv_io.F90 @@ -329,8 +329,6 @@ module subroutine ini_icepack_io(year, partit, mesh) #include "associate_mesh.h" ! Get the tracers information from icepack - call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & - nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) call icepack_query_tracer_indices( & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & @@ -369,10 +367,16 @@ module subroutine ini_icepack_io(year, partit, mesh) !----------------------------------------------------------------- ! 3D restart fields (ncat) !----------------------------------------------------------------- - call icepack_files%def_node_var('aicen' , 'sea ice concentration' , 'none', aicen(:,:) , mesh, partit, ncat) - call icepack_files%def_node_var('vicen' , 'volum per unit area of ice' , 'm' , vicen(:,:) , mesh, partit, ncat) - call icepack_files%def_node_var('vsnon' , 'volum per unit area of snow' , 'm' , vsnon(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('aice' , 'sea ice concentration' , 'none', aice(:) , mesh, partit) + call icepack_files%def_node_var('vice' , 'volum per unit area of ice' , 'm' , vice(:) , mesh, partit) + call icepack_files%def_node_var('vsno' , 'volum per unit area of snow' , 'm' , vsno(:) , mesh, partit) + + call icepack_files%def_node_var('aicen' , 'sea ice concentration per class' , 'none', aicen(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('vicen' , 'volum per unit area of ice per class' , 'm' , vicen(:,:) , mesh, partit, ncat) + call icepack_files%def_node_var('vsnon' , 'volum per unit area of snow per class' , 'm' , vsnon(:,:) , mesh, partit, ncat) call icepack_files%def_node_var('Tsfc' , 'sea ice surf. temperature' , 'degC', trcrn(:,nt_Tsfc,:), mesh, partit, ncat) + call icepack_files%def_node_var('uvel' , 'zonal component of ice velocity' , 'm/s' , uvel(:) , mesh, partit) + call icepack_files%def_node_var('vvel' , 'meridional component of ice velocity' , 'm/s' , vvel(:) , mesh, partit) if (tr_iage) then call icepack_files%def_node_var('iage' , 'sea ice age' , 's' , trcrn(:,nt_iage,:), mesh, partit, ncat) diff --git a/src/io_restart.F90 b/src/io_restart.F90 index f076a6d8e..b2297bec8 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -312,9 +312,10 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr if (.not. l_read) then call ini_ocean_io(yearnew, dynamics, tracers, partit, mesh) if (use_ice) then - call ini_ice_io(yearnew, ice, partit, mesh) #if defined(__icepack) call ini_icepack_io(yearnew, partit, mesh) +#else + call ini_ice_io(yearnew, ice, partit, mesh) #endif end if #if defined(__recom) @@ -324,9 +325,10 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr else call ini_ocean_io(yearold, dynamics, tracers, partit, mesh) if (use_ice) then - call ini_ice_io (yearold, ice, partit, mesh) #if defined(__icepack) call ini_icepack_io(yearold, partit, mesh) +#else + call ini_ice_io (yearold, ice, partit, mesh) #endif end if #if defined(__recom) @@ -379,26 +381,33 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ! read netcdf file restart else 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_restart(oce_path, oce_files, partit%MPI_COMM_FESOM, partit%mype) + + !_______________________________________________________________________ + ! read ICE/ICEPACK restart if (use_ice) then - if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ice'//achar(27)//'[0m' - call read_restart(ice_path, ice_files, partit%MPI_COMM_FESOM, partit%mype) #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_restart(icepack_path, icepack_files, partit%MPI_COMM_FESOM, partit%mype) + call read_restart(icepack_path, 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_restart(ice_path, ice_files, partit%MPI_COMM_FESOM, partit%mype) #endif end if #if defined(__recom) -!RECOM restart -!read here - 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_restart(bio_path, bio_files, partit%MPI_COMM_FESOM, partit%mype) - end if + !_______________________________________________________________________ + ! 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_restart(bio_path, bio_files, partit%MPI_COMM_FESOM, partit%mype) + end if #endif + !_______________________________________________________________________ ! immediately create a raw core dump restart if(raw_restart_length_unit /= "off") then call write_all_raw_restarts(istep, partit%MPI_COMM_FESOM, partit%mype) @@ -448,25 +457,31 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ! finally write restart for netcdf, core dump and derived type binary ! write netcdf restart if(is_portable_restart_write) then + !___________________________________________________________________________ + ! write OCEAN restart ! if(partit%mype==0) write(*,*)'Do output (netCDF, restart) ...' if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ocean'//achar(27)//'[0m' call write_restart(oce_path, oce_files, istep) + + !___________________________________________________________________________ + ! write ICE/ICEPACK restart if(use_ice) then - if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ice'//achar(27)//'[0m' - call write_restart(ice_path, ice_files, istep) #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_restart(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_restart(ic e_path, ice_files, istep) #endif end if #if defined(__recom) -!RECOM restart -!write here - 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_restart(bio_path, bio_files, istep) - end if + !___________________________________________________________________________ + ! 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_restart(bio_path, bio_files, istep) + end if #endif end if From 5801e2788a9f1ed0c6b6c9ae2d3e62e73ab60c68 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Oct 2024 13:44:58 +0200 Subject: [PATCH 16/18] fix minor issue --- src/io_restart.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io_restart.F90 b/src/io_restart.F90 index b2297bec8..a51be4bd6 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -471,7 +471,7 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr call write_restart(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_restart(ic e_path, ice_files, istep) + call write_restart(ice_path, ice_files, istep) #endif end if From d9a4ae5448db4d1a12626b2ee8056a98fea9923a Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 18 Oct 2024 15:00:10 +0200 Subject: [PATCH 17/18] change default ICEPACK branch in download_icepack.sh from icepack_fesom2 into icepack_fesom2_bugfix --- src/icepack_drivers/download_icepack.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/icepack_drivers/download_icepack.sh b/src/icepack_drivers/download_icepack.sh index a579268eb..59a2c07d9 100755 --- a/src/icepack_drivers/download_icepack.sh +++ b/src/icepack_drivers/download_icepack.sh @@ -6,5 +6,6 @@ DIR="./Icepack" if [ ! -d "$DIR" ]; then git clone https://github.com/lzampier/Icepack.git cd $DIR - git checkout icepack_fesom2 + # git checkout icepack_fesom2 + git checkout icepack_fesom2_bugfix fi From f65da846b098477d5d5c2f20d607aeae75234521 Mon Sep 17 00:00:00 2001 From: ackerlar Date: Sun, 20 Oct 2024 21:50:18 +0200 Subject: [PATCH 18/18] solve icb related segfault --- src/icb_step.F90 | 9 +++++++-- src/io_meandata.F90 | 2 +- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/icb_step.F90 b/src/icb_step.F90 index a7da3baa8..a3cfefe46 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -464,7 +464,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt i_have_element=.false. !if the first node belongs to this processor.. (just one processor enters here!) !if( local_idx_of(iceberg_elem) > 0 .and. elem2D_nodes(1,local_idx_of(iceberg_elem)) <= myDim_nod2D ) then -if((local_idx_of(iceberg_elem)>0) .and. (local_idx_of(iceberg_elem)<=partit%myDim_elem2D+partit%eDim_elem2D) ) then +if((local_idx_of(iceberg_elem)>0) .and. (local_idx_of(iceberg_elem)<=partit%myDim_elem2D) ) then if( elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D ) then i_have_element=.true. @@ -804,14 +804,19 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single write(*,*) " * set iceberg ", ib, " back to elem ", old_element, " from elem ", iceberg_elem write(*,*) " * area_ib = ", length_ib_single * width_ib_single, "; area_ib_tot = ", area_ib_tot, "; elem_area = ", elem_area_tmp end if + left_mype = 0.0 lon_rad = old_lon lat_rad = old_lat lon_deg = lon_rad/rad lat_deg = lat_rad/rad - iceberg_elem = old_element u_ib = 0. v_ib = 0. + + i_have_element= (local_idx_of(iceberg_elem) .ne. 0) + if(i_have_element) then + i_have_element= mesh%elem2D_nodes(1,local_idx_of(iceberg_elem)) <= partit%myDim_nod2D !1 PE still .true. + end if end if else if (mype==0) write(*,*) 'iceberg ',ib, ' changed PE or was very fast' diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 5b2fc8876..3596a9083 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -1049,7 +1049,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwbv', 'basal iceberg melting', 'm/s', ibfwbv(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwl', 'lateral iceberg melting', 'm/s', ibfwl(:), 1, 'm', i_real4, partit, mesh) call def_stream(nod2D, myDim_nod2D, 'ibfwe', 'iceberg erosion', 'm/s', ibfwe(:), 1, 'm', i_real4, partit, mesh) - call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'm/s', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) + call def_stream((/nl,nod2D/), (/nl,myDim_nod2D/), 'ibhf', 'heat flux from iceberg melting', 'W/m2', ibhf_n(:,:), 1, 'm', i_real4, partit, mesh) end if !------------------------------------------ !_______________________________________________________________________________