SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, & tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, & global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys, & m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir) use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM,water_reservoir_nom use comsoil_h, only: volcapa,inertiedat use adsorption_mod, only : regolith_adsorption,adsorption_pem USE temps_mod_evol, ONLY: year_PEM USE ice_table_mod, only: computeice_table_equilibrium #ifndef CPP_STD use surfdat_h, only: watercaptag #endif implicit none character(len=*), intent(in) :: filename ! name of the startfi_PEM.nc integer,intent(in) :: ngrid ! # of physical grid points integer,intent(in) :: nsoil_GCM ! # of vertical grid points in the GCM integer,intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM integer,intent(in) :: nslope ! # of sub-grid slopes integer,intent(in) :: timelen ! # time samples real, intent(in) :: tsurf_ave_yr1(ngrid,nslope) ! surface temperature at the first year of GCM call [K] real,intent(in) :: tsurf_ave_yr2(ngrid,nslope) ! surface temperature at the second year of GCM call [K] real,intent(in) :: q_co2(ngrid,timelen) ! MMR tracer co2 [kg/kg] real,intent(in) :: q_h2o(ngrid,timelen) ! MMR tracer h2o [kg/kg] real,intent(in) :: ps_inst(ngrid,timelen) ! surface pressure [Pa] real,intent(in) :: timestep ! time step [s] real,intent(in) :: tend_h2oglaciers(ngrid,nslope) ! tendencies on h2o glaciers real,intent(in) :: tend_co2glaciers(ngrid,nslope) ! tendencies on co2 glaciers real,intent(in) :: co2ice(ngrid,nslope) ! co2 ice amount [kg/m^2] real,intent(in) :: waterice(ngrid,nslope) ! water ice amount [kg/m^2] real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope) ! soil temperature during 1st year [K] real, intent(in) :: watersurf_ave(ngrid,nslope) ! surface water ice density, yearly averaged (kg/m^3) real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3) real, intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa] ! outputs real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! soil (mid-layer) thermal inertia in the PEM grid [SI] real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) ! soil (mid-layer) temperature [K] real,intent(inout) :: ice_table(ngrid,nslope) ! Ice table depth [m] real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature [K] real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope) ! mass of co2 adsorbed [kg/m^2] real,intent(out) :: deltam_co2_regolith_phys(ngrid) ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2] real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope) ! mass of h2o adsorbed [kg/m^2] real,intent(out) :: deltam_h2o_regolith_phys(ngrid) ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2] real,intent(inout) :: water_reservoir(ngrid) ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2] ! local real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope) ! soil temperature saved in the start [K] real :: TI_startPEM(ngrid,nsoil_PEM,nslope) ! soil thermal inertia saved in the start [SI] LOGICAL :: found ! check if variables are found in the start LOGICAL :: found2 ! check if variables are found in the start integer :: iloop,ig,islope,it,isoil ! index for loops REAL :: TI_breccia = 750. ! Thermal inertia of Breccia following Wood 2009 [SI] REAL :: TI_bedrock = 2300. ! Thermal inertia of Bedrock following Wood 2009 [SI] real :: kcond ! Thermal conductivity, intermediate variable [SI] real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m] CHARACTER*2 :: num ! intermediate string to read PEM start sloped variables real :: tsoil_saved(ngrid,nsoil_PEM) ! saved soil temperature [K] real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope) ! intermediate soil temperature during yr1[K] real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope) ! intermediate soil temperature during yr 2 [K] real :: alph_tmp(ngrid,nsoil_PEM-1) ! Intermediate for tsoil computation [] real :: beta_tmp(ngrid,nsoil_PEM-1) ! Intermediate for tsoil computatio [] real :: year_PEM_read ! Year of the PEM previous run real :: alpha_clap_h2o = -6143.7 ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [K^-1] real :: beta_clap_h2o = 28.9074 ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [1] LOGICAL :: startpem_file ! boolean to check if we read the startfile or not !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: read start_pem. Need a specific iostart_PEM !!! !!! Order: 0. Previous year of the PEM run !!! 1. Thermal Inertia !!! 2. Soil Temperature !!! 3. Ice table !!! 4. Mass of CO2 & H2O adsorbed !!! 5. Water reservoir !!! !!! /!\ This order must be respected ! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !0. Check if the start_PEM exist. inquire(file=filename,exist = startpem_file) write(*,*)'Is start PEM?',startpem_file startpem_file = .true. !1. Run if (startpem_file) then ! open pem initial state file: call open_startphy(filename) if(soil_pem) then call get_var("Time",year_PEM_read,found) year_PEM=INT(year_PEM_read) if(.not.found) then write(*,*)'PEMetat0: failed loading year_PEM; take default=0' else write(*,*)'year_PEM of startpem=', year_PEM endif !1. Thermal Inertia ! a. General case DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found) if(.not.found) then write(*,*)'PEM settings: failed loading ' write(*,*)'will reconstruct the values of TI_PEM' do ig = 1,ngrid if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then !!! transition delta = 50. TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ & (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ & ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2)))) do iloop=nsoil_GCM+2,n_1km TI_PEM(ig,iloop,islope) = TI_breccia enddo else ! we keep the high ti values do iloop=nsoil_GCM+1,n_1km TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope) enddo endif ! TI PEM and breccia comparison !! transition delta = 1000. TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ & (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ & ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2)))) do iloop=n_1km+2,nsoil_PEM TI_PEM(ig,iloop,islope) = TI_bedrock enddo enddo else do iloop = nsoil_GCM+1,nsoil_PEM TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope) ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here. enddo endif ! not found ENDDO ! islope print *,'PEMETAT0: THERMAL INERTIA DONE' ! b. Special case for inertiedat, inertiedat_PEM call get_field("inertiedat_PEM",inertiedat_PEM,found) if(.not.found) then do iloop = 1,nsoil_GCM inertiedat_PEM(:,iloop) = inertiedat(:,iloop) enddo !!! zone de transition delta = 50. do ig = 1,ngrid if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ & (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ & ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2)))) do iloop = nsoil_GCM+2,n_1km inertiedat_PEM(ig,iloop) = TI_breccia enddo else do iloop=nsoil_GCM+1,n_1km inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM) enddo endif ! comparison ti breccia enddo!ig !!! zone de transition delta = 1000. do ig = 1,ngrid inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ & (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ & ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2)))) enddo do iloop = n_1km+2, nsoil_PEM do ig = 1,ngrid inertiedat_PEM(ig,iloop) = TI_bedrock enddo enddo endif ! not found !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !2. Soil Temperature DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found) if(.not.found) then write(*,*)'PEM settings: failed loading ' write(*,*)'will reconstruct the values of Tsoil' do ig = 1,ngrid kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1)) do iloop=nsoil_GCM+2,n_1km kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM)) enddo kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) do iloop=n_1km+2,nsoil_PEM kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km)) enddo do iloop = nsoil_GCM+1,nsoil_PEM tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM,islope) enddo enddo ! call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope)) ! do iloop = 1,2000000 ! write(*,*) 'iloop,islope',islope,iloop ! call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope)) ! enddo else ! predictor corrector: restart from year 1 of the GCM and build the evolution of ! tsoil at depth tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope) call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope) call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope)) do iloop = nsoil_GCM+1,nsoil_PEM tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope) enddo endif do it = 1,timelen do isoil = nsoil_GCM+1,nsoil_PEM tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope) enddo enddo do isoil = nsoil_GCM+1,nsoil_PEM do ig = 1,ngrid watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) enddo enddo ENDDO ! islope print *,'PEMETAT0: SOIL TEMP DONE' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !3. Ice Table call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table) call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM) print *,'PEMETAT0: ICE TABLE DONE' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !4. CO2 & H2O Adsorption if(adsorption_pem) then DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found) if((.not.found)) then m_co2_regolith_phys(:,:,:) = 0. endif exit ENDDO DO islope=1,nslope write(num,fmt='(i2.2)') islope call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found) if((.not.found2)) then m_h2o_regolith_phys(:,:,:) = 0. endif exit ENDDO call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, & m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys) if((.not.found)) then deltam_co2_regolith_phys(:) = 0. endif if((.not.found2)) then deltam_h2o_regolith_phys(:) = 0. endif print *,'PEMETAT0: CO2 & H2O adsorption done ' endif endif ! soil_pem !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !. 5 water reservoir #ifndef CPP_STD call get_field("water_reservoir",water_reservoir,found) if(.not.found) then write(*,*)'Pemetat0: failed loading ' write(*,*)'will reconstruct the values from watercaptag' do ig=1,ngrid if(watercaptag(ig)) then water_reservoir(ig)=water_reservoir_nom else water_reservoir(ig)=0. endif enddo endif #endif call close_startphy else !No startfi, let's build all by hand year_PEM=0 if(soil_pem) then !a) Thermal inertia do islope = 1,nslope do ig = 1,ngrid if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then !!! transition delta = 50. TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ & (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ & ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2)))) do iloop=nsoil_GCM+2,n_1km TI_PEM(ig,iloop,islope) = TI_breccia enddo else ! we keep the high ti values do iloop=nsoil_GCM+1,n_1km TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope) enddo endif !! transition delta = 1000. TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ & (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ & ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2)))) do iloop=n_1km+2,nsoil_PEM TI_PEM(ig,iloop,islope) = TI_bedrock enddo enddo enddo do iloop = 1,nsoil_GCM inertiedat_PEM(:,iloop) = inertiedat(:,iloop) enddo !!! zone de transition delta = 50. do ig = 1,ngrid if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ & (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ & ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2)))) do iloop = nsoil_GCM+2,n_1km inertiedat_PEM(ig,iloop) = TI_breccia enddo else do iloop = nsoil_GCM+1,n_1km inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM) enddo endif enddo !!! zone de transition delta = 1000. do ig = 1,ngrid inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ & (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ & ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2)))) enddo do iloop = n_1km+2, nsoil_PEM do ig = 1,ngrid inertiedat_PEM(ig,iloop) = TI_bedrock enddo enddo print *,'PEMETAT0: TI DONE' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !b) Soil temperature do islope = 1,nslope do ig = 1,ngrid kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1)) do iloop=nsoil_GCM+2,n_1km kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM)) enddo kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) do iloop=n_1km+2,nsoil_PEM kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km)) enddo do iloop = nsoil_GCM+1,nsoil_PEM tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM,islope) enddo enddo ! call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope)) ! do iloop = 1,2000000 ! write(*,*) 'islope,iloop',islope,iloop ! call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope)) ! enddo do it = 1,timelen do isoil = nsoil_GCM+1,nsoil_PEM tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope) enddo enddo do isoil = nsoil_GCM+1,nsoil_PEM do ig = 1,ngrid watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) enddo enddo enddo print *,'PEMETAT0: TSOIL DONE ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !c) Ice table call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table) call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! print *,'PEMETAT0: Ice table DONE ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !d) Regolith adsorbed if(adsorption_pem) then m_co2_regolith_phys(:,:,:) = 0. m_h2o_regolith_phys(:,:,:) = 0. call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, & m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys) deltam_co2_regolith_phys(:) = 0. deltam_h2o_regolith_phys(:) = 0. endif print *,'PEMETAT0: CO2 adsorption done ' endif !soil_pem !. e) water reservoir #ifndef CPP_STD do ig=1,ngrid if(watercaptag(ig)) then water_reservoir=water_reservoir_nom else water_reservoir=0. endif enddo #endif endif ! of if (startphy_file) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(soil_pem) then !! Sanity check DO ig = 1,ngrid DO islope = 1,nslope DO iloop = 1,nsoil_PEM if(isnan(tsoil_PEM(ig,iloop,islope))) then call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1) endif ENDDO ENDDO ENDDO endif!soil_pem END SUBROUTINE