MODULE clim_state_init !----------------------------------------------------------------------- ! NAME ! clim_state_init ! ! DESCRIPTION ! Read the start files to initialize the climate state. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4 use io_netcdf, only: open_nc, close_nc, get_var_nc, get_dim_nc, put_var_nc, start_name, start1D_name, startfi_name, startpem_name ! DECLARATION ! ----------- implicit none contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE read_start() !----------------------------------------------------------------------- ! NAME ! read_start ! ! DESCRIPTION ! Read the file "start.nc". ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: nlayer, nlon, nlat, ngrid use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM use tracers, only: nq, qnames, set_q_PCM use stoppage, only: stop_clean use display, only: print_msg ! DECLARATION ! ----------- implicit none ! LOCAL VARIABLES ! --------------- real(dp), dimension(nlayer + 1) :: tmp1d real(dp), dimension(nlon + 1,nlat) :: tmp2d real(dp), dimension(nlon + 1,nlat,nlayer) :: tmp3d logical(k4) :: here integer(di) :: i ! CODE ! ---- ! In case of 1D ! ~~~~~~~~~~~~~ if (ngrid == 1) then call read_start1D() return end if ! In case of 3D ! ~~~~~~~~~~~~~ ! Open inquire(file = start_name,exist = here) if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start_name//'"!',1) call print_msg('> Reading "'//start_name//'"') call open_nc(start_name,'read') ! Get the variables call get_var_nc('ap',tmp1d) call set_ap(tmp1d) call get_var_nc('bp',tmp1d) call set_bp(tmp1d) call get_var_nc('ps',tmp2d) call set_ps_PCM(tmp2d) call get_var_nc('teta',tmp3d) call set_teta_PCM(tmp3d) do i = 1,nq call get_var_nc(trim(qnames(i)),tmp3d) call set_q_PCM(tmp3d,i) end do ! Close call close_nc(start_name) END SUBROUTINE read_start !======================================================================= !======================================================================= SUBROUTINE read_start1D() !----------------------------------------------------------------------- ! NAME ! read_start1D ! ! DESCRIPTION ! Read the file "start1D.txt". ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: nlayer, nslope use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM, set_u_PCM, set_v_PCM, compute_alt_coord, set_ap, set_bp use tracers, only: nq, set_q_PCM use stoppage, only: stop_clean use config, only: read_callphys use display, only: print_msg ! DECLARATION ! ----------- implicit none ! LOCAL VARIABLES ! --------------- integer(di) :: i, j, k, ierr, funit character(30) :: header real(dp), dimension(1,1) :: tmp real(dp), dimension(nslope) :: tmp_1d real(dp), dimension(1,1,nlayer) :: q_tmp, teta_tmp, wind_tmp real(dp), dimension(nlayer + 1) :: ap_tmp, bp_tmp real(dp) :: pa, preff logical(k4) :: here ! CODE ! ---- ! Open "start1D.txt" inquire(file = start1D_name,exist = here) if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start1D_name//'"!',1) call print_msg('> Reading "'//start1D_name//'"') open(newunit = funit,file = start1D_name,status = "old",action = "read",iostat = ierr) if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//start1D_name//'"!',ierr) ! Get the variables read(funit,*) header, tmp, pa, preff call set_ps_PCM(tmp) do i = 1,nq read(funit,*,iostat = ierr) header, (tmp_1d(j),j = 1,nslope), (q_tmp(1,1,k),k = 1,nlayer) if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'not enough atmospheric layers defined in the file "'//start1D_name//'" for the tracer "'//trim(header)//'"!',1) call set_q_PCM(q_tmp,i) end do read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer) call set_u_PCM(wind_tmp) read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer) call set_v_PCM(wind_tmp) read(funit,*,iostat = ierr) header, (tmp_1d(j),j = 1,nslope), (teta_tmp(1,1,k),k = 1,nlayer) call set_teta_PCM(teta_tmp) ! Close close(funit) ! Construct altitude coordinates (not stored in "start1D.txt") call compute_alt_coord(pa,preff,ap_tmp,bp_tmp) call set_ap(ap_tmp) call set_bp(bp_tmp) ! Read the "callphys.def" call read_callphys() ! To get 'CO2cond_ps' END SUBROUTINE read_start1D !======================================================================= !======================================================================= SUBROUTINE read_startfi() !----------------------------------------------------------------------- ! NAME ! read_startfi ! ! DESCRIPTION ! Read the file "startfi.nc". ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil_PCM use stoppage, only: stop_clean use config, only: read_controldata use slopes, only: set_def_slope_mean, set_subslope_dist, set_iflat use surface, only: set_albedodat_PCM, set_albedo_PCM, set_emissivity_PCM use surf_temp, only: set_tsurf_PCM use soil_temp, only: set_tsoil_PCM, set_flux_geo_PCM use frost, only: set_h2ofrost_PCM, set_co2frost_PCM use surf_ice, only: set_is_h2o_perice_PCM, set_co2_perice_PCM use soil, only: set_TI_PCM, set_inertiedat_PCM use display, only: print_msg ! DECLARATION ! ----------- implicit none ! LOCAL VARIABLES ! --------------- real(dp), dimension(:), allocatable :: tmp1d real(dp), dimension(:,:), allocatable :: tmp2d real(dp), dimension(ngrid,nsoil_PCM,nslope) :: tmp3d logical(k4) :: here ! CODE ! ---- inquire(file = startfi_name,exist = here) if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//startfi_name//'"!',1) ! Allocate the array to store the variables call print_msg('> Reading "'//startfi_name//'"') allocate(tmp1d(nslope + 1),tmp2d(ngrid,nslope)) ! Get control data call read_controldata() ! Open call open_nc(startfi_name,'read') ! Get the variables if (nslope > 1) then call get_var_nc('def_slope_mean',tmp1d) call set_def_slope_mean(tmp1d) call get_var_nc('subslope_dist',tmp2d) call set_subslope_dist(tmp2d) end if call get_var_nc('flux_geo',tmp2d) call set_flux_geo_PCM(tmp2d) call get_var_nc('h2o_ice',tmp2d) where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer call set_h2ofrost_PCM(tmp2d) call get_var_nc('co2',tmp2d) where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer call set_co2frost_PCM(tmp2d) call get_var_nc('perennial_co2ice',tmp2d) where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer call set_co2_perice_PCM(tmp2d) deallocate(tmp1d); allocate(tmp1d(ngrid)) call get_var_nc('watercaptag',tmp1d) call set_is_h2o_perice_PCM(tmp1d) call get_var_nc('albedodat',tmp1d) call set_albedodat_PCM(tmp1d) call get_var_nc('albedo',tmp2d) call set_albedo_PCM(tmp2d) call get_var_nc('emis',tmp2d) call set_emissivity_PCM(tmp2d) call get_var_nc('tsurf',tmp2d) call set_tsurf_PCM(tmp2d) call get_var_nc('tsoil',tmp3d) call set_tsoil_PCM(tmp3d) deallocate(tmp2d); allocate(tmp2d(ngrid,nsoil_PCM)) call get_var_nc('inertiedat',tmp2d) call set_inertiedat_PCM(tmp2d) call get_var_nc('inertiesoil',tmp3d) call set_TI_PCM(tmp3d) ! To do? ! h2oice_depth ! d_coef ! lag_co2_ice ! Close/Deallocate call close_nc(startfi_name) deallocate(tmp1d,tmp2d) END SUBROUTINE read_startfi !======================================================================= !======================================================================= SUBROUTINE read_startpem(tsurf_avg_yr1,tsurf_avg_yr2,ps_avg_glob,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,d_h2oice,d_co2ice,h2o_ice,co2_ice, & tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map, & h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) !----------------------------------------------------------------------- ! NAME ! read_startpem ! ! DESCRIPTION ! Read the file "startpem.nc" which stores the PEM state. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use stoppage, only: stop_clean use geometry, only: ngrid, nslope, nsoil, nsoil_PCM use evolution, only: dt use physics, only: mugaz, r use planet, only: TI_breccia, TI_bedrock, alpha_clap_h2o, beta_clap_h2o use frost, only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM use surf_ice, only: is_h2o_perice_PCM, h2oice_huge_ini, co2_perice_PCM use soil, only: do_soil, TI, depth_breccia, depth_bedrock, index_breccia, index_bedrock, inertiedat, layer, inertiedat_PCM use soil_temp, only: ini_tsoil_pem, compute_tsoil use soil_therm_inertia, only: update_soil_TI use ice_table, only: icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium use sorption, only: do_sorption, evolve_regolith_adsorption use tracers, only: mmol, iPCM_qh2o use layered_deposits, only: layering, do_layering, array2map, ini_layering, add_stratum use surf_ice, only: rho_co2ice, rho_h2oice use display, only: print_msg use utility, only: int2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: tsurf_avg_yr1, tsurf_avg_yr2 ! Surface temperature at the last but one PCM run and at the last PCM run real(dp), intent(in) :: ps_avg_glob ! Mean average pressure on the planet real(dp), dimension(:,:), intent(in) :: ps_ts ! Surface pressure timeseries real(dp), dimension(:,:), intent(in) :: q_co2_ts, q_h2o_ts ! MMR of CO2 and H2O real(dp), dimension(:,:), intent(in) :: h2o_surfdensity_avg ! Average of surface water density real(dp), dimension(:,:), intent(in) :: d_h2oice, d_co2ice ! Tendencies on ice real(dp), dimension(:,:), intent(inout) :: h2o_ice, co2_ice ! Surface ice real(dp), dimension(:,:,:), intent(inout) :: tsoil_avg ! Soil temperature real(dp), dimension(:,:,:), intent(inout) :: h2o_soildensity_avg ! Average of soil water soil density real(dp), dimension(:,:), intent(inout) :: icetable_depth ! Ice table depth real(dp), dimension(:,:), intent(inout) :: icetable_thickness ! Ice table thickness real(dp), dimension(:,:,:), intent(inout) :: ice_porefilling ! Subsurface ice pore filling real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg ! Mass of H2O and CO2 adsorbed type(layering), dimension(:,:), intent(inout) :: layerings_map ! Layerings real(dp), dimension(:), intent(inout) :: delta_h2o_ads, delta_co2_ads ! Mass of H2O and CO2 exchanged due to adsorption/desorption ! LOCAL VARIABLES ! --------------- logical(k4) :: here integer(di) :: i, islope, k, nb_str_max, nsoil_startpem real(dp) :: delta ! Depth of the interface regolith/breccia, breccia/bedrock [m] real(dp), dimension(ngrid,nsoil,nslope) :: TI_startpem ! Soil thermal inertia saved in the startpem [SI] real(dp), dimension(ngrid,nsoil,nslope) :: tsoil_startpem ! Soil temperature saved in the startpem [K] real(dp), dimension(:,:,:,:), allocatable :: layerings_array ! Array for layerings ! CODE ! ---- ! Check if the file exists inquire(file = startpem_name,exist = here) ! If the file is not here ! ~~~~~~~~~~~~~~~~~~~~~~~ if (.not. here) then ! It is possibly because it is the very first PEM run so everything needs to be initalized call print_msg('> The file "'//startpem_name//'" was not found (possibly because this is the first PEM run)') ! H2O ice call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true.") h2o_ice(:,:) = 0._dp do i = 1,ngrid if (is_h2o_perice_PCM(i)) h2o_ice(i,:) = h2oice_huge_ini + h2ofrost_PCM(i,:) - h2o_frost4PCM(i,:) end do ! CO2 ice call print_msg("'co2_ice' is initialized with 'perennial_co2ice' found in the PCM.") co2_ice(:,:) = co2_perice_PCM(:,:) + co2frost_PCM(:,:) - co2_frost4PCM(:,:) if (do_soil) then ! Thermal inertia do islope = 1,nslope do i = 1,ngrid if (TI(i,index_breccia,islope) < TI_breccia) then !!! transition delta = depth_breccia TI(i,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/ & (((delta - layer(index_breccia))/(TI(i,index_breccia,islope)**2)) + & ((layer(index_breccia + 1) - delta)/(TI_breccia**2)))) do k = index_breccia + 2,index_bedrock TI(i,k,islope) = TI_breccia end do else ! we keep the high ti values do k = index_breccia + 1,index_bedrock TI(i,k,islope) = TI(i,index_breccia,islope) end do end if !! transition delta = depth_bedrock TI(i,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/ & (((delta - layer(index_bedrock))/(TI(i,index_bedrock,islope)**2)) + & ((layer(index_bedrock + 1) - delta)/(TI_breccia**2)))) do k = index_bedrock + 2,nsoil TI(i,k,islope) = TI_bedrock end do end do end do do k = 1,nsoil_PCM inertiedat(:,k) = inertiedat_PCM(:,k) end do !!! zone de transition delta = depth_breccia do i = 1,ngrid if (inertiedat(i,index_breccia) < TI_breccia) then inertiedat(i,index_breccia + 1) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/ & (((delta - layer(index_breccia))/(inertiedat(i,index_breccia)**2)) + & ((layer(index_breccia + 1) - delta)/(TI_breccia**2)))) do k = index_breccia + 2,index_bedrock inertiedat(i,k) = TI_breccia end do else do k = index_breccia + 1,index_bedrock inertiedat(i,k) = inertiedat(i,index_breccia) end do end if end do !!! zone de transition delta = depth_bedrock do i = 1,ngrid inertiedat(i,index_bedrock + 1) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/ & (((delta - layer(index_bedrock))/(inertiedat(i,index_bedrock)**2)) + & ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2)))) end do do k = index_bedrock + 2,nsoil do i = 1,ngrid inertiedat(i,k) = TI_bedrock end do end do ! Soil temperature do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope)) call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope)) ! First raw initialization h2o_soildensity_avg(:,nsoil_PCM + 1:nsoil,islope) = exp(beta_clap_h2o/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope) + alpha_clap_h2o)/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope)**mmol(iPCM_qh2o)/(mugaz*r) end do ! Ice table if (icetable_equilibrium) then call computeice_table_equilibrium(is_h2o_perice_PCM,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness) call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope)) end do else if (icetable_dynamic) then call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope)) end do end if ! Absorption in regolith if (do_sorption) then h2o_ads_reg = 0._dp co2_ads_reg = 0._dp call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) delta_co2_ads = 0._dp delta_h2o_ads = 0._dp end if ! do_sorption end if ! do_soil ! Layering if (do_layering) then call print_msg('layerings_map is initialized with sub-surface strata.') call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.") do i = 1,ngrid if (is_h2o_perice_PCM(i)) then do islope = 1,nslope call ini_layering(layerings_map(i,islope)) call add_stratum(layerings_map(i,islope),1.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,h2oice_huge_ini/rho_h2oice,0.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,0._dp) end do else do islope = 1,nslope call ini_layering(layerings_map(i,islope)) if (co2_perice_PCM(i,islope) > 0.) call add_stratum(layerings_map(i,islope),1.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0._dp) end do end if end do end if ! do_layering return end if ! If the file is present ! ~~~~~~~~~~~~~~~~~~~~~~ call print_msg('> Reading "'//startpem_name//'"') call open_nc(startpem_name,'read') ! H2O ice h2o_ice(:,:) = 0._dp call get_var_nc('h2o_ice',h2o_ice(:,:)) ! CO2 ice co2_ice(:,:) = 0._dp call get_var_nc('co2_ice',co2_ice(:,:)) if (do_soil) then ! Check if the number of soil layers is compatible call get_dim_nc('subsurface_layers',nsoil_startpem) if (nsoil_startpem /= nsoil) then call print_msg('nsoil (PEM) = '//int2str(nsoil)//' | nsoil ("'//startpem_name//'") = '//int2str(nsoil_startpem)) call stop_clean(__FILE__,__LINE__,'nsoil defined in the PEM is different from the one read in "'//startpem_name//'"!',1) end if do islope = 1,nslope ! Soil temperature call get_var_nc('tsoil',tsoil_startpem(:,:,islope)) ! Predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly ! for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope)) call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope)) call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope)) end do tsoil_avg(:,nsoil_PCM + 1:nsoil,:) = tsoil_startpem(:,nsoil_PCM + 1:nsoil,:) if (any(isnan(tsoil_avg(:,:,:)))) call stop_clean(__FILE__,__LINE__,"NaN detected in 'tsoil_avg'",1) ! Thermal inertia call get_var_nc('TI',TI_startpem(:,:,:)) TI(:,nsoil_PCM + 1:nsoil,:) = TI_startpem(:,nsoil_PCM + 1:nsoil,:) ! 1st layers can change because of the presence of ice at the surface, so we don't change it here call get_var_nc('inertiedat',inertiedat(:,:)) ! Ice table call get_var_nc('icetable_depth',icetable_depth(:,:)) if (icetable_dynamic) call get_var_nc('ice_porefilling',ice_porefilling(:,:,:)) if (icetable_equilibrium) then call get_var_nc('icetable_thickness',icetable_thickness(:,:)) else if (icetable_dynamic) then call get_var_nc('ice_porefilling',ice_porefilling(:,:,:)) end if ! Absorption in regolith if (do_sorption) then call get_var_nc('co2_ads_reg',co2_ads_reg(:,:,:)) call get_var_nc('h2o_ads_reg',h2o_ads_reg(:,:,:)) call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) end if ! do_sorption end if ! do_soil ! Layering if (do_layering) then call get_dim_nc('nb_str_max',nb_str_max) allocate(layerings_array(ngrid,nslope,nb_str_max,6)) layerings_array = 0. call get_var_nc('stratif_top_elevation',layerings_array(:,:,:,1)) call get_var_nc('stratif_h_co2ice',layerings_array(:,:,:,2)) call get_var_nc('stratif_h_h2oice',layerings_array(:,:,:,3)) call get_var_nc('stratif_h_dust',layerings_array(:,:,:,4)) call get_var_nc('stratif_h_pore',layerings_array(:,:,:,5)) call get_var_nc('stratif_poreice_volfrac',layerings_array(:,:,:,6)) call array2map(layerings_array,layerings_map) deallocate(layerings_array) end if ! Close call close_nc(startpem_name) END SUBROUTINE read_startpem !======================================================================= END MODULE clim_state_init