MODULE climate_init !----------------------------------------------------------------------- ! NAME ! climate_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, startevo_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 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, LVL_NFO ! 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 ! ---- ! 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//'"',LVL_NFO) 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_startfi() !----------------------------------------------------------------------- ! NAME ! read_startfi ! ! DESCRIPTION ! Read the file "startfi.nc". ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! C. Metz, 04/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil_PCM use stoppage, only: stop_clean use config, only: read_control_data 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 soil, only: set_TI_PCM, set_inertiedat_PCM, set_soildepth_PCM use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! LOCAL VARIABLES ! --------------- real(dp), dimension(:), allocatable :: tmp1d real(dp), dimension(:,:), allocatable :: tmp2d real(dp), dimension(ngrid,nsoil_PCM,nslope) :: tmp3d integer(di) :: islope 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//'"',LVL_NFO) allocate(tmp1d(nslope + 1),tmp2d(ngrid,nslope)) ! Get control data call read_control_data() ! 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 ! cmetz for now we have to hack the startfi to input a value of 0, because it doesn't exist in the gPCM !call get_var_nc('flux_geo',tmp2d) tmp2d(:,:) = 0._dp 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_ice',tmp2d) where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer call set_co2frost_PCM(tmp2d) deallocate(tmp1d); allocate(tmp1d(ngrid)) 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) ! cmetz we don't touch intertie soil in the gPEM but sure, we just keep it and import/export the same thing call get_var_nc('inertiesoil',tmp2d) do islope = 1,nslope tmp3d(:,:,islope) = tmp2d(:,:) end do call set_TI_PCM(tmp3d) deallocate(tmp1d); allocate(tmp1d(nsoil_PCM)) call get_var_nc('soildepth',tmp1d) call set_soildepth_PCM(tmp1d) ! Close/Deallocate call close_nc(startfi_name) deallocate(tmp1d,tmp2d) END SUBROUTINE read_startfi !======================================================================= !======================================================================= SUBROUTINE read_startevo(tsurf_avg_yr1,tsurf_avg_yr2,ps_avg_glob,ps_ts,q_co2_ts,q_h2o_ts,h2o_ice,co2_ice,tsoil_avg) !----------------------------------------------------------------------- ! NAME ! read_startevo ! ! DESCRIPTION ! Read the file "startevo.nc" which stores the PEM state. ! ! AUTHORS & DATE ! L. Lange, 2023 ! JB Clement, 2023-2026 ! C. Metz, 04/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use stoppage, only: stop_clean use geometry, only: ngrid, nslope, nsoil, nsoil_PCM use evolution, only: dt use frost, only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM 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: TI_breccia, TI_bedrock use slopes, only: subslope_dist, def_slope_mean use maths, only: pi use display, only: print_msg, LVL_WRN, LVL_NFO, LVL_ERR 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(inout) :: h2o_ice, co2_ice ! Surface ice real(dp), dimension(:,:,:), intent(inout) :: tsoil_avg ! Soil temperature ! LOCAL VARIABLES ! --------------- logical(k4) :: here integer(di) :: i, islope, k, nsoil_startevo real(dp) :: delta ! Depth of the interface regolith/breccia, breccia/bedrock [m] real(dp), dimension(ngrid,nsoil,nslope) :: TI_startevo ! Soil thermal inertia saved in the startevo [SI] real(dp), dimension(ngrid,nsoil,nslope) :: tsoil_startevo ! Soil temperature saved in the startevo [K] ! CODE ! ---- ! Check if the file exists inquire(file = startevo_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 "'//startevo_name//'" was not found (possibly because this is the first PEM run)',LVL_WRN) ! H2O ice call print_msg("'h2o_ice' is initialized with the yearly minimum of H2O ice found in the PCM.",LVL_WRN) h2o_ice(:,:) = h2ofrost_PCM(:,:) - h2o_frost4PCM(:,:) ! CO2 ice call print_msg("'co2_ice' is initialized with the yearly minimum of CO2 ice found in the PCM.",LVL_WRN) co2_ice(:,:) = 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)) end do end if ! do_soil return end if ! If the file is present ! ~~~~~~~~~~~~~~~~~~~~~~ call print_msg('> Reading "'//startevo_name//'"',LVL_NFO) call open_nc(startevo_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_startevo) if (nsoil_startevo /= nsoil) then call print_msg('nsoil (PEM) = '//int2str(nsoil)//' | nsoil ("'//startevo_name//'") = '//int2str(nsoil_startevo),LVL_ERR) call stop_clean(__FILE__,__LINE__,'nsoil defined in the PEM is different from the one read in "'//startevo_name//'"!',1) end if do islope = 1,nslope ! Soil temperature call get_var_nc('tsoil',tsoil_startevo(:,:,islope)) ! Predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startevo 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_startevo(:,:,islope)) call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startevo(:,:,islope)) call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_startevo(:,:,islope)) end do tsoil_avg(:,nsoil_PCM + 1:nsoil,:) = tsoil_startevo(:,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_startevo(:,:,:)) TI(:,nsoil_PCM + 1:nsoil,:) = TI_startevo(:,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(:,:)) end if ! do_soil ! Close call close_nc(startevo_name) END SUBROUTINE read_startevo !======================================================================= END MODULE climate_init