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. diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 55a259930..07b1dba78 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -357,23 +357,35 @@ 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, & - dhi_dt_out=thdgrsn, & - dhs_dt_out=thdgr, & - evap_ocn_out=evaporation ) - - heat_flux(:) = - net_heat_flux(:) - water_flux(:) = - (fresh_wa_flux(:)/1000.0_WP) - runoff(:) - - ! Evaporation - evaporation(:) = - evaporation(:) / 1000.0_WP - ice_sublimation(:) = 0.0_WP + 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 + ! Heat flux + heat_flux(n) = - net_heat_flux(n) + + ! 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(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(n) = - ice_sublimation(n) * inv_rhowat + end do +!$OMP END PARALLEL DO call init_flux_atm_ocn() @@ -515,12 +527,31 @@ 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 + ! -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) & + flux(n) = evaporation(n) & -ice_sublimation(n) & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean - +prec_rain(n) & + +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 @@ -614,6 +645,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 @@ -622,7 +655,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/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 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_io.F90 b/src/icepack_drivers/icedrv_io.F90 index 7d8bb7c4d..e14732c05 100644 --- a/src/icepack_drivers/icedrv_io.F90 +++ b/src/icepack_drivers/icedrv_io.F90 @@ -18,28 +18,32 @@ 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 + use mod_mesh + use io_meandata, only: def_stream - 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, & @@ -50,251 +54,264 @@ module subroutine init_io_icepack(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 ') - 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 ') - 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 ') - 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 ') - 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 ') + !_______________________________________________________________________ + ! define output streams + do i=1, io_listsize + select case (trim(io_list_icepack(i)%id)) + case ('aice0 ') + 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_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_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_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_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_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_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_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_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_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_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 ') + 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 ') + 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') + 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') + 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 ') + 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 ') + 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 + 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 + 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.) - end do - case ('qsnon ') - do k = 1,nslyr ! Separate variable for each snow layer + 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.) - end do - ! Average over categories - case ('iage ') + 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 ') + 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 ') + 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 ') + 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 ') + 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 ') + 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 ') + 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 + 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 + 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) - end do - case ('qsno ') - do k = 1,nslyr ! Separate variable for each snow layer + 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) - 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 + 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_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_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 + end do ! --> do i=1, io_listsize - 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 @@ -312,8 +329,6 @@ module subroutine init_restart_icepack(year, 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, & @@ -341,68 +356,63 @@ 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('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 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,:), mesh, partit, ncat) 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,:) , mesh, partit, ncat) 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,:), 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 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,:), 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 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,:), 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 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,:), 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 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,:), mesh, partit, ncat) + call icepack_files%def_node_var('first_ice', 'distinguishes ice that disappears' , 'logical', first_ice_real(:,:), mesh, partit, ncat) end if !----------------------------------------------------------------- @@ -415,11 +425,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,:), mesh, partit, ncat) 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,:), mesh, partit, ncat) end do ! Snow @@ -428,7 +438,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,:), mesh, partit, ncat) end do ! @@ -438,6 +448,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..6428aad72 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 @@ -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) @@ -820,7 +824,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 +841,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) @@ -883,20 +889,23 @@ 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 use icepack_intfc, only: icepack_compute_tracers 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 5139454bf..df527cd86 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__) @@ -772,12 +776,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 +859,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 +886,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 +904,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) @@ -935,7 +943,6 @@ subroutine ocn_mixed_layer_icepack( & 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) @@ -951,6 +958,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, & @@ -986,17 +997,39 @@ subroutine ocn_mixed_layer_icepack( & 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 !======================================================================= - 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 +1062,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 @@ -1179,8 +1212,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 @@ -1202,18 +1235,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(:) - dhi_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 !----------------------------------------------------------------- @@ -1271,19 +1305,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 !----------------------------------------------------------------- @@ -1294,7 +1348,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/icepack_drivers/icedrv_transfer.F90 b/src/icepack_drivers/icedrv_transfer.F90 index 15d36eb1b..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 @@ -50,12 +44,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 +61,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 +87,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) @@ -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_out) ) evap_out = evap end subroutine icepack_to_fesom diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 39dd8d247..48c8bd395 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 @@ -43,9 +44,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_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 + 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 @@ -54,6 +55,7 @@ module io_fesom_file_module integer, save :: m_nod2d integer, save :: m_elem2d integer, save :: m_nl + integer, save :: m_ncat type fesom_file_type_ptr @@ -103,7 +105,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_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 @@ -111,6 +113,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_ncat type(t_partit), target :: partit ! EO parameters type(fesom_file_type_ptr), allocatable :: tmparr(:) @@ -121,9 +124,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_ncat ... icepack number of ice thickness classes, + if (present(mesh_ncat)) then + m_ncat = mesh_ncat + else + m_ncat = 0 + end if + call this%netcdf_file_type%initialize() allocate(this%used_mesh_dims(0)) @@ -182,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 @@ -224,8 +243,15 @@ 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 + ! 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 + 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 @@ -401,7 +427,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 @@ -441,11 +471,30 @@ 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(:) 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, 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) :: ncat! todo: be able to set precision + ! EO parameters + type(dim_info) level_diminfo, ncat_diminfo + +!PS write(*,*) "--> specify_node_var_2dicepack:", __LINE__, __FILE__, size(local_data) + level_diminfo = obtain_diminfo(this, m_nod2d) + ncat_diminfo = obtain_diminfo(this, size(local_data, dim=2)) + + 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 subroutine specify_node_var_3d(this, name, longname, units, local_data) @@ -456,7 +505,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)) @@ -474,6 +524,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(:) @@ -490,6 +541,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)) @@ -513,14 +565,16 @@ function obtain_diminfo(this, len) result(info) end do ! 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 - info = dim_info( idx=this%add_dim('nz', len), len=len) + else if(len == m_nl ) then + info = dim_info( idx=this%add_dim('nz' , 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 @@ -536,14 +590,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 @@ -558,6 +614,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 3596a9083..f58cc86d1 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 ! !-------------------------------------------------------------------------------------------- ! @@ -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) @@ -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 @@ -1742,7 +1762,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 @@ -1770,13 +1790,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 init_io_icepack(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..a51be4bd6 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,16 @@ 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,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 @@ -262,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 @@ -302,18 +310,31 @@ 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 +#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) 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 +#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) 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 @@ -360,22 +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 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) +#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) + 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) @@ -425,22 +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 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(ice_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 diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index 3c7a327ca..5d5445897 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -44,24 +44,32 @@ 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 - - subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit) + subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit, ncat) 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 - - 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) +!PS write(*,*) "--> def_node_var_3d:", __LINE__, __FILE__ + 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 @@ -74,7 +82,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 @@ -89,18 +97,22 @@ 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 - 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_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_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_ncat ! EO parameters type(restart_file_type), pointer :: f @@ -111,7 +123,16 @@ 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_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 ! this is specific for a restart file f%iter_varindex = f%add_var_int('iter', [f%time_dimindex()]) end subroutine @@ -131,19 +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 - - subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit) + subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mesh, partit, ncat) 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 - - 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