MODULE pemetat0_mod implicit none !======================================================================= contains !======================================================================= SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, & tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & m_h2o_regolith_phys,deltam_h2o_regolith_phys) use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var use comsoil_h_PEM, only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock use comsoil_h, only: volcapa, inertiedat use adsorption_mod, only: regolith_adsorption, adsorption_pem use ice_table_mod, only: computeice_table_equilibrium, icetable_equilibrium use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock use soil_thermalproperties_mod, only: update_soil_thermalproperties use tracer_mod, only: mmol, igcm_h2o_vap ! tracer names and molar masses use abort_pem_mod, only: abort_pem use compute_soiltemp_mod, only: ini_tsoil_pem, compute_tsoil_pem use comslope_mod, only: def_slope_mean, subslope_dist #ifndef CPP_STD use comcstfi_h, only: r, mugaz, pi use surfdat_h, only: watercaptag, watercap, perennial_co2ice #else use comcstfi_mod, only: r, mugaz, pi #endif implicit none include "callkeys.h" character(len=*), intent(in) :: filename ! name of the startfi_PEM.nc integer, intent(in) :: ngrid ! # of physical grid points integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM 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) :: timestep ! time step [s] real, intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa] real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K] real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K] real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg] real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg] real, dimension(ngrid,timelen), intent(in) :: ps_inst ! surface pressure [Pa] real, dimension(ngrid,nslope), intent(in) :: tend_h2o_ice ! tendencies on h2o ice real, dimension(ngrid,nslope), intent(in) :: tend_co2_ice ! tendencies on co2 ice real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3) ! outputs real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2] real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2] real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! h2o ice amount [kg/m^2] real, dimension(ngrid,nslope), intent(out) :: co2_ice ! co2 ice amount [kg/m^2] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K] real, dimension(ngrid,nslope), intent(inout) :: ice_table ! Ice table depth [m] real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m] real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! surface water ice density, yearly averaged (kg/m^3) ! local real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM ! soil temperature saved in the start [K] real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! 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 :: 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, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1 ! intermediate soil temperature during yr 1 [K] real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2 ! intermediate soil temperature during yr 2 [K] logical :: startpem_file ! boolean to check if we read the startfile or not #ifdef CPP_STD logical, dimension(ngrid) :: watercaptag watercaptag = .false. #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: read start_pem. Need a specific iostart_PEM !!! !!! Order: 0. Previous year of the PEM run !!! Ice initialization !!! 1. Thermal Inertia !!! 2. Soil Temperature !!! 3. Ice table !!! 4. Mass of CO2 & H2O adsorbed !!! !!! /!\ This order must be respected ! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !0.1 Check if the start_PEM exist. inquire(file = filename,exist = startpem_file) write(*,*)'Is start PEM?',startpem_file !0.2 Set to default values ice_table = -1. ! by default, no ice table ice_table_thickness = -1. !1. Run if (startpem_file) then write(*,*)'There is a "startpem.nc"!' ! open pem initial state file: call open_startphy(filename) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifndef CPP_STD ! h2o ice h2o_ice = 0. call get_field("h2o_ice",h2o_ice,found) if (.not. found) then write(*,*)'Pemetat0: failed loading ' write(*,*)'will reconstruct the values from watercaptag' write(*,*)'with default value ''ini_huge_h2oice''' do ig = 1,ngrid if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.) enddo else ! The variations of infinite reservoirs during the PCM years are taken into account h2o_ice = h2o_ice + watercap endif ! co2 ice co2_ice = perennial_co2ice #endif if (soil_pem) then !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !1. Thermal Inertia ! a. General case do islope = 1,nslope write(num,'(i2.2)') islope call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found) if (.not. found) then write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>' write(*,*)'will reconstruct the values of TI_PEM' do ig = 1,ngrid if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then !!! transition delta = depth_breccia TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + & ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) do iloop=index_breccia+2,index_bedrock TI_PEM(ig,iloop,islope) = TI_breccia enddo else ! we keep the high TI values do iloop = index_breccia + 1,index_bedrock TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope) enddo endif ! TI PEM and breccia comparison !! transition delta = depth_bedrock TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + & ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2)))) do iloop = index_bedrock + 2,nsoil_PEM TI_PEM(ig,iloop,islope) = TI_bedrock enddo enddo else ! found do iloop = nsoil_PCM + 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 write(*,*) '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_PCM inertiedat_PEM(:,iloop) = inertiedat(:,iloop) enddo !!! zone de transition delta = depth_breccia do ig = 1,ngrid if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ & ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) do iloop = index_breccia+2,index_bedrock inertiedat_PEM(ig,iloop) = TI_breccia enddo else do iloop=index_breccia+1,index_bedrock inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM) enddo endif ! comparison ti breccia enddo ! ig !!! zone de transition delta = depth_bedrock do ig = 1,ngrid inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ & (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ & ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2)))) enddo do iloop = index_bedrock + 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,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa ! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1)) ! do iloop=index_breccia+2,index_bedrock ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia)) ! enddo ! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa ! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) ! ! do iloop=index_bedrock+2,nsoil_PEM ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock)) ! enddo ! enddo call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) else ! predictor corrector: restart from year 1 of the PCM and build the evolution of ! tsoil at depth tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope) call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope)) tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope) call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope)) do iloop = nsoil_PCM+1,nsoil_PEM tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope) enddo endif !found do it = 1,timelen do isoil = nsoil_PCM+1,nsoil_PEM tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope) enddo enddo do isoil = nsoil_PCM+1,nsoil_PEM do ig = 1,ngrid watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r) enddo enddo enddo ! islope write(*,*) 'PEMETAT0: SOIL TEMP done' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !3. Ice Table if(icetable_equilibrium) then call get_field("ice_table",ice_table,found) if (.not. found) then write(*,*)'PEM settings: failed loading ' write(*,*)'will reconstruct the values of the ice table given the current state' call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness) call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) enddo endif write(*,*) 'PEMETAT0: ICE TABLE done' endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !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. exit endif enddo do islope=1,nslope write(num,fmt='(i2.2)') islope call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2) if (.not. found2) then m_h2o_regolith_phys = 0. exit endif enddo call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,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) deltam_co2_regolith_phys = 0. if (.not.found2) deltam_h2o_regolith_phys = 0. write(*,*) 'PEMETAT0: CO2 & H2O adsorption done' endif ! adsorption_pem endif ! soil_pem call close_startphy else !No startfi, let's build all by hand ! h2o ice h2o_ice = 0. write(*,*)'There is no "startpem.nc"!' write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.' do ig = 1,ngrid if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.) enddo ! co2 ice write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.' co2_ice = perennial_co2ice if (soil_pem) then !a) Thermal inertia do islope = 1,nslope do ig = 1,ngrid if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then !!! transition delta = depth_breccia TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ & ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) do iloop=index_breccia + 2,index_bedrock TI_PEM(ig,iloop,islope) = TI_breccia enddo else ! we keep the high ti values do iloop=index_breccia + 1,index_bedrock TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope) enddo endif !! transition delta = depth_bedrock TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ & (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ & ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2)))) do iloop = index_bedrock + 2,nsoil_PEM TI_PEM(ig,iloop,islope) = TI_bedrock enddo enddo enddo do iloop = 1,nsoil_PCM inertiedat_PEM(:,iloop) = inertiedat(:,iloop) enddo !!! zone de transition delta = depth_breccia do ig = 1,ngrid if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2)))) do iloop = index_breccia + 2,index_bedrock inertiedat_PEM(ig,iloop) = TI_breccia enddo else do iloop = index_breccia + 1,index_bedrock inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia) enddo endif enddo !!! zone de transition delta = depth_bedrock do ig = 1,ngrid inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + & ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2)))) enddo do iloop = index_bedrock + 2,nsoil_PEM do ig = 1,ngrid inertiedat_PEM(ig,iloop) = TI_bedrock enddo enddo write(*,*) 'PEMETAT0: TI done' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !b) Soil temperature do islope = 1,nslope ! do ig = 1,ngrid ! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa ! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1)) ! ! do iloop=index_breccia+2,index_bedrock ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia)) ! enddo ! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa ! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) ! ! do iloop=index_bedrock+2,nsoil_PEM ! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa ! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock)) ! enddo ! enddo call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) ! First raw initialization do it = 1,timelen do isoil = nsoil_PCM + 1,nsoil_PEM tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope) enddo enddo do it = 1,timelen do isoil = nsoil_PCM + 1,nsoil_PEM call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it)) enddo enddo do isoil = nsoil_PCM + 1,nsoil_PEM do ig = 1,ngrid watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r) enddo enddo enddo !islope write(*,*) 'PEMETAT0: TSOIL done' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !c) Ice table if(icetable_equilibrium) then call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness) call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) enddo write(*,*) 'PEMETAT0: Ice table done' endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !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_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,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 write(*,*) 'PEMETAT0: CO2 adsorption done' endif !soil_pem 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))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1) enddo enddo enddo endif !soil_pem END SUBROUTINE pemetat0 END MODULE pemetat0_mod