!------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil ! I_f Compute tendencies & Save initial situation ! I_g Save initial PCM situation ! I_h Read the PEMstart ! I_i Compute orbit criterion ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 glaciers flows ! II_d Update surface and soil temperatures ! II_e Update the tendencies ! II_f Checking the stopping criterion ! III Output ! III_a Update surface value for the PCM start files ! III_b Write start and starfi.nc ! III_c Write start_pem !------------------------ PROGRAM pem !module needed for INITIALISATION #ifndef CPP_STD use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa use surfdat_h, only: tsurf, co2ice, emis,& qsurf,watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & hmons, summit, base,albedo_h2o_frost, & frost_albedo_threshold,emissiv use dimradmars_mod, only: totcloudfrac, albedo use dust_param_mod, only: tauscaling use tracer_mod, only: noms,igcm_h2o_ice ! tracer names #else use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, & emissiv use tracer_h, only: noms,igcm_h2o_ice,igcm_co2_ice ! tracer names use phys_state_var_mod #endif use phyetat0_mod, only: phyetat0 use phyredem, only: physdem0, physdem1 use turb_mod, only: q2, wstar use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, & nf90_get_var, nf90_inq_varid, nf90_inq_dimid, & nf90_inquire_dimension,nf90_close ! For phyredem : USE control_mod, ONLY: iphysiq, day_step,nsplit_phys USE iniphysiq_mod, ONLY: iniphysiq USE logic_mod, ONLY: iflag_phys #ifndef CPP_STD use mod_phys_lmdz_para, only: is_parallel, is_sequential, & is_mpi_root, is_omp_root, & is_master use planete_h, only: aphelie, periheli, year_day, peri_day, & obliquit #else ! USE comcstfi_mod, ONLY: rad,g,r,cpp,pi ! USE inifis_mod, ONLY: inifis use planete_mod, only: apoastr, periastr, year_day, peri_day, & obliquit #endif USE mod_const_mpi, ONLY: COMM_LMDZ USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & subslope_dist,co2ice_slope, & tsurf_slope,tsoil_slope,fluxgrd_slope,& fluxrad_sky_slope,sky_slope,callsubslope,& co2iceflow, beta_slope, capcal_slope,& albedo_slope,emiss_slope,qsurf_slope,& iflat,major_slope,ini_comslope_h use time_phylmdz_mod, only: daysec,dtphys USE comconst_mod, ONLY: rad,g,r,cpp,pi USE infotrac USE geometry_mod, only: latitude_deg use pemredem, only: pemdem1 use co2glaciers_mod,only: co2glaciers_evol !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, & TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, & ! soil thermal inertia tsoil_PEM, mlayer_PEM,layer_PEM, & !number of subsurface layers, soil mid layer depths fluxgeo, co2_adsorbded_phys ! geothermal flux, mass of co2 in the regolith use adsorption_mod, only : regolith_co2adsorption !!! For orbit parameters USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem use orbit_param_criterion_mod, only : orbit_param_criterion use recomp_orb_param_mod, only: recomp_orb_param IMPLICIT NONE include "dimensions.h" include "paramet.h" INTEGER ngridmx PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) include "comdissnew.h" include "comgeom.h" include "iniprint.h" ! Same variable's name as in the GCM INTEGER :: ngrid !Number of physical grid points INTEGER :: nlayer !Number of vertical layer INTEGER :: nq !Number of tracer INTEGER :: day_ini !First day of the simulation REAL :: pday !Physical day REAL :: time_phys !Same as GCM REAL :: ptimestep !Same as GCM REAL :: ztime_fin !Same as GCM ! Variable for reading start.nc character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM ! variables dynamiques REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants REAL teta(ip1jmp1,llm) ! temperature potentielle REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes REAL ps(ip1jmp1) ! pression au sol REAL, dimension(:),allocatable :: ps_phys !(ngrid) ! pression au sol REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression au sol instantannées REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression au sol instantannées for the first year of the gcm REAL masse(ip1jmp1,llm) ! masse d'air REAL phis(ip1jmp1) ! geopotentiel au sol REAL time_0 ! Variable for reading starfi.nc character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM integer :: ncid, varid,status !Variable for handling opening of files integer :: phydimid, subdimid, nlayerdimid, nqdimid !Variable ID for Netcdf files integer :: lonvarid, latvarid, areavarid,sdvarid !Variable ID for Netcdf files integer :: apvarid,bpvarid !Variable ID for Netcdf files ! Variable for reading starfi.nc and writting restartfi.nc REAL, dimension(:),allocatable :: longitude !Longitude read in FILE_NAME and written in restartfi REAL, dimension(:),allocatable :: latitude !Latitude read in FILE_NAME and written in restartfi REAL, dimension(:),allocatable :: ap !Coefficient ap read in FILE_NAME_start and written in restart REAL, dimension(:),allocatable :: bp !Coefficient bp read in FILE_NAME_start and written in restart REAL, dimension(:),allocatable :: cell_area !Cell_area read in FILE_NAME and written in restartfi REAL :: Total_surface !Total surface of the planet ! Variable for h2o_ice evolution REAL , dimension(:,:), allocatable :: tendencies_h2o_ice ! LON x LAT field : Tendency of evolution of perenial ice over a year REAL, dimension(:),allocatable :: tendencies_h2o_ice_phys ! physical point field : Tendency of evolution of perenial ice over a year REAL , dimension(:,:), allocatable :: tendencies_co2_ice ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year REAL, dimension(:),allocatable :: tendencies_co2_ice_phys ! physical point field : Tendency of evolution of perenial co2 ice over a year REAL :: ini_surf ! Initial surface of sublimating water ice REAL :: ini_surf_h2o ! Initial surface of sublimating water ice REAL, dimension(:),allocatable :: initial_h2o_ice ! physical point field : Logical array indicating sublimating point REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice REAL, dimension(:),allocatable :: initial_co2_ice ! physical point field : Logical array indicating sublimating point of co2 ice REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year REAL , dimension(:,:), allocatable :: min_co2_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year REAL , dimension(:,:), allocatable :: min_co2_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year REAL :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] REAL :: global_ave_press_old ! constant: Global average pressure of initial/previous time step REAL :: global_ave_press_new ! constant: Global average pressure of current time step REAL , dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] REAL , dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step LOGICAL :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate? LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate? REAL, dimension(:,:,:),allocatable :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg] real,save :: m_co2, m_noco2, A , B, mmean ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg] real ,allocatable :: vmr_co2_gcm_phys(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM real ,allocatable :: q_h2o_GCM_phys(:,:) ! Physics x Times h2o mass mixing ratio in the first layer from the GCM [kg/kg] real ,allocatable :: q_co2_GCM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer from the GCM [kg/kg] real ,allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM [kg/kg] REAL, ALLOCATABLE :: ps_GCM(:,:,:) ! Lon x Lat x Times: surface pressure from the GCM [Pa] REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:) ! Lon x Lat x Times: surface pressure from the 1st year of the GCM [Pa] REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:) ! Lon x Lat x Times: co2 volumemixing ratio retrieve from the gcm [m^3/m^3] REAL, ALLOCATABLE :: q_h2o_GCM(:,:,:) ! Lon x Lat x Times: h2o volume mixing ratio retrieved from the GCM REAL ,allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM integer :: timelen ! # time samples REAL :: ave ! intermediate varibale to compute average REAL, ALLOCATABLE :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) REAL :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 REAL :: beta_clap_co2 = 3182.48 ! clapeyron's law for CO2 REAL :: alpha_clap_co2 = 23.3494 ! Clapeyron's law for CO2 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE REAL ,allocatable :: watercap_slope(:,:) ! Physics x Nslope: watercap per slope REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2] REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2] REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year [kg/m^2] REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year [kg/m^2] REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim_slope ! physical point field : Logical array indicating sublimating point of co2 ice REAL, dimension(:,:),allocatable :: initial_h2o_ice_slope ! physical point field : Logical array indicating if there is water ice at initial state REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating if there is co2 ice at initial state REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope ! LON x LAT x nslope field : Tendency of evolution of perenial co2 ice over a year REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope ! LON x LAT x slope field : Tendency of evolution of perenial water ice over a year REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM REAL, dimension(:,:),allocatable :: tendencies_h2o_ice_phys_slope ! physical pointx slope field : Tendency of evolution of perenial co2 ice REAL , dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope) ! To flag where there is a glacier flow REAL , dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) ! To flag where there is a glacier flow !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K] REAL, ALLOCATABLE :: tsurf_ave_phys(:,:) ! Physic x LAT x SLOPE field : Averaged Surface Temperature [K] REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K] REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature during 1st year of the GCM [K] REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature during 1st year [K] REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia [SI] REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) ! LON X LAT x SLOPE XTULES field : Surface Temperature in timeseries [K] REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) ! Physic x SLOPE XTULES field : NOn averaged Surf Temperature [K] REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K] REAL, ALLOCATABLE :: tsurf_ave_phys_yr1(:,:) ! Physic SLOPE field : Averaged Surface Temperature of first call of the gcm [K] REAL, ALLOCATABLE :: inertiesoil(:,:) !Physic x Depth Thermal inertia of the mesh for restart [SI] REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope [SI] REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Same but for the start REAL,ALLOCATABLE :: ice_depth(:,:) ! Physic x SLope: Ice table depth [m] REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] REAL,ALLOCATABLE :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] REAL,ALLOCATABLE :: alph_locslope(:,:) ! Physic x Soil: Intermediate to compute Tsoil [1] REAL,ALLOCATABLE :: beta_locslope(:,:) ! Physic x Soil : Intermediate tocompute Tsoil [K] REAL,ALLOCATABLE :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3] REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3] REAL,ALLOCATABLE :: watersurf_density_phys_timeseries(:,:,:) ! Physic x Slope x Times, water surface density, time series [kg/m^3] REAL,ALLOCATABLE :: watersurf_density_phys_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] REAL,ALLOCATABLE :: watersoil_density_phys_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] REAL,ALLOCATABLE :: watersoil_density_phys_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] REAL, ALLOCATABLE :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] REAL :: totmass_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K] REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1] LOGICAL :: bool_sublim ! logical to check if there is sublimation or not !! Some parameters for the PEM run REAL, PARAMETER :: year_step = 1 ! timestep for the pem INTEGER :: year_iter ! number of iteration INTEGER :: year_iter_max ! maximum number of iterations before stopping REAL :: timestep ! timestep [s] #ifdef CPP_STD ! INTEGER :: nsplit_phys=1 ! LOGICAL :: iflag_phys=.true. REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice REAL :: albedo_h2o_frost ! albedo of h2o frost REAL,ALLOCATABLE :: co2ice(:) ! Physics: co2 ice mesh averaged [kg/m^2] #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Loop variable INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil #ifndef CPP_STD ! Parallel variables is_sequential=.true. is_parallel=.false. is_mpi_root=.true. is_omp_root=.true. is_master=.true. #endif day_ini=0 !test time_phys=0. !test ! Some constants ngrid=ngridmx nlayer=llm m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 year_day=669 daysec=88775. dtphys=0 timestep=year_day*daysec/year_step !------------------------ ! I Initialisation ! I_a READ run.def !------------------------ !----------------------------READ run.def --------------------- CALL conf_gcm( 99, .TRUE. ) call infotrac_init nq=nqtot !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc !------------------------ !----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- ! In the gcm, these values are given to the physic by the dynamic. ! Here we simply read them in the startfi_evol.nc file status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) allocate(longitude(ngrid)) allocate(latitude(ngrid)) allocate(cell_area(ngrid)) status = nf90_inq_varid(ncid, "longitude", lonvarid) status = nf90_get_var(ncid, lonvarid, longitude) status = nf90_inq_varid(ncid, "latitude", latvarid) status = nf90_get_var(ncid, latvarid, latitude) status = nf90_inq_varid(ncid, "area", areavarid) status = nf90_get_var(ncid, areavarid, cell_area) call ini_comsoil_h(ngrid) status = nf90_inq_varid(ncid, "soildepth", sdvarid) status = nf90_get_var(ncid, sdvarid, mlayer) status =nf90_close(ncid) !----------------------------READ start.nc --------------------- allocate(q(ip1jmp1,llm,nqtot)) CALL dynetat0(FILE_NAME_start,vcov,ucov, & teta,q,masse,ps,phis, time_0) CALL iniconst !new CALL inigeom allocate(ap(nlayer+1)) allocate(bp(nlayer+1)) status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid) status = nf90_inq_varid(ncid, "ap", apvarid) status = nf90_get_var(ncid, apvarid, ap) status = nf90_inq_varid(ncid, "bp", bpvarid) status = nf90_get_var(ncid, bpvarid, bp) status =nf90_close(ncid) CALL iniphysiq(iim,jjm,llm, & (jjm-1)*iim+2,comm_lmdz, & daysec,day_ini,dtphys/nsplit_phys, & rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, & iflag_phys) !----------------------------READ startfi.nc --------------------- ! First we read the initial state (starfi.nc) allocate(watercap_slope(ngrid,nslope)) allocate(TI_GCM_start(ngrid,nsoilmx,nslope)) allocate(inertiesoil(ngrid,nsoilmx)) #ifndef CPP_STD CALL phyetat0 (FILE_NAME,0,0, & nsoilmx,ngrid,nlayer,nq, & day_ini,time_phys, & tsurf,tsoil,albedo,emis, & q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, & watercap,inertiesoil,nslope,tsurf_slope, & tsoil_slope,co2ice_slope,def_slope,def_slope_mean, & subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start, & qsurf_slope,watercap_slope) if(soil_pem) then deallocate(TI_GCM_start) !not used then endif ! Remove unphysical values of surface tracer DO i=1,ngrid DO nnq=1,nqtot DO islope=1,nslope if(qsurf_slope(i,nnq,islope).LT.0) then qsurf_slope(i,nnq,islope)=0. endif enddo enddo enddo #else call phys_state_var_init(nq) IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) call initracer(ngrid,nq) call iniaerosol() call phyetat0(.true., & ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, & day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, & cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab, tsea_ice,sea_ice) call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) call ini_comslope_h(ngrid,nsoilmx,nq) allocate(co2ice(ngrid)) co2ice(:)=qsurf(:,igcm_co2_ice) co2ice_slope(:,1)=co2ice(:) tsurf_slope(:,1)=tsurf(:) if (nslope.eq.1) then def_slope(1) = 0 def_slope(2) = 0 def_slope_mean=0 subslope_dist(:,1) = 1. endif ! Remove unphysical values of surface tracer DO i=1,ngrid DO nnq=1,nqtot qsurf_slope(i,nnq,1)=qsurf(i,nnq) if(qsurf(i,nnq).LT.0) then qsurf(i,nnq)=0. endif enddo enddo #endif DO nnq=1,nqtot if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq ENDDO !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation !------------------------ !----------------------------Subslope parametrisation definition --------------------- ! Define some slope statistics iflat=1 DO islope=2,nslope IF(abs(def_slope_mean(islope)).lt. & abs(def_slope_mean(iflat))) THEN iflat = islope ENDIF ENDDO PRINT*,'Flat slope for islope = ',iflat PRINT*,'corresponding criterium = ',def_slope_mean(iflat) allocate(flag_co2flow(ngrid,nslope)) allocate(flag_co2flow_mesh(ngrid)) flag_co2flow(:,:) = 0. flag_co2flow_mesh(:) = 0. !---------------------------- READ GCM data --------------------- ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid !------------------------ ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value call nb_time_step_GCM("data_GCM_Y1.nc",timelen) allocate(min_h2o_ice_s_1(iim+1,jjm+1)) allocate(min_co2_ice_s_1(iim+1,jjm+1)) allocate(vmr_co2_gcm(iim+1,jjm+1,timelen)) allocate(q_h2o_GCM(iim+1,jjm+1,timelen)) allocate(q_co2_GCM(iim+1,jjm+1,timelen)) allocate(ps_GCM(iim+1,jjm+1,timelen)) allocate(ps_GCM_yr1(iim+1,jjm+1,timelen)) allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope)) allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope)) allocate(tsurf_ave(iim+1,jjm+1,nslope)) allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope)) allocate(tsurf_ave_phys(ngrid,nslope)) allocate(tsurf_ave_phys_yr1(ngrid,nslope)) allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen)) allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen)) allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen)) allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen)) allocate(Tsurfave_before_saved(ngrid,nslope)) allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope)) allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope)) allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope)) allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope)) allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen)) allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(delta_co2_adsorbded(ngrid)) allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen)) allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen)) allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen)) allocate(watersurf_density_phys_ave(ngrid,nslope)) allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope)) print *, "Downloading data Y1..." call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm, min_h2o_ice_s_1,min_co2_ice_s_1,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,& nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, & watersurf_density_timeseries,watersoil_density_timeseries) ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value print *, "Downloading data Y1 done" allocate(min_h2o_ice_s_2(iim+1,jjm+1)) allocate(min_co2_ice_s_2(iim+1,jjm+1)) allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope)) allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope)) print *, "Downloading data Y2" call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm ,min_h2o_ice_s_2,min_co2_ice_s_2,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, & nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, & watersurf_density_timeseries,watersoil_density_timeseries) print *, "Downloading data Y2 done" ! The variables in the dynamic grid are transfered to the physical grid allocate(vmr_co2_gcm_phys(ngrid,timelen)) allocate(vmr_co2_pem_phys(ngrid,timelen)) allocate(q_h2o_GCM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) allocate(q_co2_GCM_phys(ngrid,timelen)) allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(ps_phys(ngrid)) allocate(ps_phys_timeseries(ngrid,timelen)) allocate(ps_phys_timeseries_yr1(ngrid,timelen)) CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys) CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys) CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys) call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys) call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries) call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1) CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys) CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1) deallocate(vmr_co2_gcm) deallocate(q_h2o_GCM) deallocate(q_co2_GCM) deallocate(ps_GCM) deallocate(ps_GCM_yr1) deallocate(tsurf_ave) deallocate(tsurf_ave_yr1) q_co2_PEM_phys(:,:)= q_co2_GCM_phys(:,:) q_h2o_PEM_phys(:,:)= q_h2o_GCM_phys(:,:) !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil !------------------------ !---------------------------- Initialisation of the PEM soil and values --------------------- call end_comsoil_h_PEM call ini_comsoil_h_PEM(ngrid,nslope) allocate(ice_depth(ngrid,nslope)) ice_depth(:,:) = 0. allocate(TI_GCM_phys(ngrid,nsoilmx,nslope)) DO islope = 1,nslope if(soil_pem) then CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope)) endif !soil_pem DO t=1,timelen CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t)) enddo ENDDO deallocate(co2_ice_GCM_slope) deallocate(TI_GCM) deallocate(tsurf_GCM_timeseries) if(soil_pem) then call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM) DO islope = 1,nslope DO t=1,timelen CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t)) ENDDO DO l=1,nsoilmx CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope)) DO t=1,timelen CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t)) ENDDO ENDDO DO l=nsoilmx+1,nsoilmx_PEM CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope)) DO t=1,timelen watersoil_density_phys_PEM_timeseries(:,l,islope,t) = watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t) ENDDO ENDDO ENDDO watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen deallocate(watersurf_density_timeseries) deallocate(watersurf_density_phys_timeseries) deallocate(watersoil_density_timeseries) deallocate(tsoil_ave_yr1) deallocate(tsoil_ave) deallocate(tsoil_GCM_timeseries) endif !soil_pem !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil ! I_f Compute tendencies & Save initial situation !----- Compute tendencies from the PCM run allocate(tendencies_h2o_ice(iim+1,jjm+1)) allocate(tendencies_h2o_ice_phys(ngrid)) allocate(tendencies_co2_ice(iim+1,jjm+1)) allocate(tendencies_co2_ice_phys(ngrid)) allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope)) allocate(tendencies_co2_ice_phys_slope(ngrid,nslope)) allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope)) allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope)) allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope)) ! Compute the tendencies of the evolution of ice over the years call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,& min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys) call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,& min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys) call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,& min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope) tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:) call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,& min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope) !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil ! I_f Compute tendencies & Save initial situation ! I_g Save initial PCM situation !---------------------------- Save initial PCM situation --------------------- allocate(initial_h2o_ice(ngrid)) allocate(initial_co2_ice(ngrid)) allocate(initial_co2_ice_sublim_slope(ngrid,nslope)) allocate(initial_co2_ice_slope(ngrid,nslope)) allocate(initial_h2o_ice_slope(ngrid,nslope)) ! We save the places where water ice is sublimating ! We compute the surface of water ice sublimating ini_surf=0. ini_surf_co2=0. ini_surf_h2o=0. Total_surface=0. do i=1,ngrid Total_surface=Total_surface+cell_area(i) if (tendencies_h2o_ice_phys(i).LT.0) then initial_h2o_ice(i)=1. ini_surf=ini_surf+cell_area(i) else initial_h2o_ice(i)=0. endif do islope=1,nslope if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then initial_co2_ice_sublim_slope(i,islope)=1. ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) else initial_co2_ice_sublim_slope(i,islope)=0. endif if (co2ice_slope(i,islope).GT.0) then initial_co2_ice_slope(i,islope)=1. else initial_co2_ice_slope(i,islope)=0. endif if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then initial_h2o_ice_slope(i,islope)=1. ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) else initial_h2o_ice_slope(i,islope)=0. endif enddo enddo print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2 print *, "Total initial surface of h2o ice sublimating=", ini_surf print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o print *, "Total surface of the planet=", Total_surface allocate(zplev_gcm(ngrid,nlayer+1)) DO l=1,nlayer+1 DO ig=1,ngrid zplev_gcm(ig,l) = ap(l) + bp(l)*ps_phys(ig) ENDDO ENDDO global_ave_press_old=0. do i=1,ngrid global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface enddo global_ave_press_GCM=global_ave_press_old global_ave_press_new=global_ave_press_old print *, "Initial global average pressure=", global_ave_press_GCM !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil ! I_f Compute tendencies & Save initial situation ! I_g Save initial PCM situation ! I_h Read the PEMstart !---------------------------- Read the PEMstart --------------------- call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, & tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,& tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,& watersurf_density_phys_ave,watersoil_density_phys_PEM_ave) if(soil_pem) then totmass_adsorbded = 0. do ig = 1,ngrid do islope =1, nslope do l = 1,nsoilmx_PEM - 1 totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & cell_area(ig) enddo enddo enddo write(*,*) "Tot mass in the regolith=", totmass_adsorbded deallocate(tsoil_ave_phys_yr1) endif !soil_pem deallocate(tsurf_ave_phys_yr1) deallocate(ps_phys_timeseries_yr1) !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialisation of the PEM variable and soil ! I_f Compute tendencies & Save initial situation ! I_g Save initial PCM situation ! I_h Read the PEMstar ! I_i Compute orbit criterion #ifndef CPP_STD CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) #else CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit) #endif if(evol_orbit_pem) then call orbit_param_criterion(year_iter_max) else year_iter_max=Max_iter_pem endif !--------------------------- END INITIALISATION --------------------- !---------------------------- RUN --------------------- !------------------------ ! II Run ! II_a update pressure,ice and tracers !------------------------ year_iter=0 do while (year_iter.LT.year_iter_max) ! II.a.1. Compute updated global pressure print *, "Recomputing the new pressure..." do i=1,ngrid do islope=1,nslope global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface enddo enddo print *, 'Global average pressure old time step',global_ave_press_old print *, 'Global average pressure new time step',global_ave_press_new do i=1,ngrid if(soil_pem) then global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface endif enddo print *, 'Global average pressure old time step',global_ave_press_old print *, 'Global average pressure new time step',global_ave_press_new ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen)) print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..." DO l=1,nlayer+1 DO ig=1,ngrid zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:) ENDDO ENDDO ! II.a.3. Surface pressure timeseries print *, "Recomputing the surface pressure timeserie adapted to the new pressure..." do i = 1,ngrid ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old enddo ! II.a.4. New pressure levels timeseries allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen)) print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..." do l=1,nlayer+1 do ig=1,ngrid zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:) enddo enddo ! II.a.5. Tracers timeseries print *, "Recomputing of tracer VMR timeseries for the new pressure..." l=1 DO ig=1,ngrid DO t = 1, timelen q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) if(q_h2o_PEM_phys(ig,t).LT.0) then q_h2o_PEM_phys(ig,t)=1E-30 endif if(q_h2o_PEM_phys(ig,t).GT.1) then q_h2o_PEM_phys(ig,t)=1. endif enddo enddo DO ig=1,ngrid DO t = 1, timelen q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) & + ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - & (zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / & (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) if (q_co2_PEM_phys(ig,t).LT.0) then q_co2_PEM_phys(ig,t)=1E-30 elseif (q_co2_PEM_phys(ig,t).GT.1) then q_co2_PEM_phys(ig,t)=1. endif mmean=1/(A*q_co2_PEM_phys(ig,t) +B) vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 ENDDO ENDDO deallocate(zplev_new_timeseries) deallocate(zplev_old_timeseries) ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II.b. Evolution of the ice print *, "Evolution of h2o ice" call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope) print *, "Evolution of co2 ice" call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope) !------------------------ ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 glaciers flows !------------------------ print *, "Co2 glacier flows" call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,& global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh) !------------------------ ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 glaciers flows ! II_d Update surface and soil temperatures !------------------------ ! II_d.1 Update Tsurf print *, "Updating the new Tsurf" bool_sublim=0 Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:) DO ig = 1,ngrid DO islope = 1,nslope if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice if(latitude_deg(ig).gt.0) then do ig_loop=ig,ngrid DO islope_loop=islope,iflat,-1 if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop) bool_sublim=1 exit endif enddo if (bool_sublim.eq.1) then exit endif enddo else do ig_loop=ig,1,-1 DO islope_loop=islope,iflat if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop) bool_sublim=1 exit endif enddo if (bool_sublim.eq.1) then exit endif enddo endif initial_co2_ice_slope(ig,islope)=0 if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then albedo_slope(ig,1,islope) = albedo_h2o_frost albedo_slope(ig,2,islope) = albedo_h2o_frost else albedo_slope(ig,1,islope) = albedodat(ig) albedo_slope(ig,2,islope) = albedodat(ig) emiss_slope(ig,islope) = emissiv endif elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2 ave=0. do t=1,timelen if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.)) else ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t) endif enddo tsurf_ave_phys(ig,islope)=ave/timelen endif enddo enddo do t = 1,timelen tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:)) enddo ! for the start do ig = 1,ngrid do islope = 1,nslope tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope)) enddo enddo if(soil_pem) then ! II_d.2 Update soil temperature allocate(TI_locslope(ngrid,nsoilmx_PEM)) allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) allocate(Tsurf_locslope(ngrid)) allocate(alph_locslope(ngrid,nsoilmx_PEM-1)) allocate(beta_locslope(ngrid,nsoilmx_PEM-1)) print *,"Updating soil temperature" ! Soil averaged do islope = 1,nslope TI_locslope(:,:) = TI_PEM(:,:,islope) alph_locslope(:,:) = alph_PEM(:,:,islope) beta_locslope(:,:) = beta_PEM(:,:,islope) do t = 1,timelen Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) Tsurf_locslope(:) = tsurf_phys_GCM_timeseries(:,islope,t) call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) do ig = 1,ngrid do isoil = 1,nsoilmx_PEM watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil) if(isnan(Tsoil_locslope(ig,isoil))) then write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil) stop endif enddo enddo enddo alph_PEM(:,:,islope) = alph_locslope(:,:) beta_PEM(:,:,islope) = beta_locslope(:,:) enddo tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen print *, "Update of soil temperature done" deallocate(TI_locslope) deallocate(Tsoil_locslope) deallocate(Tsurf_locslope) deallocate(alph_locslope) deallocate(beta_locslope) write(*,*) "Compute ice table" ! II_d.3 Update the ice table call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth) print *, "Update soil propreties" ! II_d.4 Update the soil thermal properties call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, & ice_depth,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbded call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, & co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded) endif !soil_pem !------------------------ ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 glaciers flows ! II_d Update surface and soil temperatures ! II_e Update the tendencies !------------------------ print *, "Adaptation of the new co2 tendencies to the current pressure" call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,& global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) !------------------------ ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 glaciers flows ! II_d Update surface and soil temperatures ! II_e Update the tendencies ! II_f Checking the stopping criterion !------------------------ call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope) call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope) year_iter=year_iter+dt_pem print *, "Checking all the stopping criterion." if (STOPPING_water) then print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water endif if (year_iter.ge.year_iter_max) then print *, "STOPPING because maximum number of iterations reached" endif if (STOPPING_1_water) then print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water endif if (STOPPING_co2) then print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2 endif if (STOPPING_1_co2) then print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2 endif if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then exit else print *, "We continue!" print *, "Number of iteration of the PEM : year_iter=", year_iter endif global_ave_press_old=global_ave_press_new enddo ! big time iteration loop of the pem !---------------------------- END RUN PEM --------------------- !---------------------------- OUTPUT --------------------- !------------------------ ! III Output ! III_a Update surface value for the PCM start files !------------------------ ! III_a.1 Ice update (for startfi) ! Co2 ice DO ig = 1,ngrid co2ice(ig) = 0. DO islope = 1,nslope co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) & * subslope_dist(ig,islope) / & cos(pi*def_slope_mean(islope)/180.) ENDDO #ifdef CPP_STD qsurf(ig,igcm_co2_ice)=co2ice(ig) #endif ENDDO ! of DO ig=1,ngrid ! H2o ice DO ig = 1,ngrid qsurf(ig,igcm_h2o_ice) = 0. DO islope = 1,nslope qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) & * subslope_dist(ig,islope) / & cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO ! of DO ig=1,ngrid ! III_a.2 Tsoil update (for startfi) if(soil_pem) then call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys) tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen) else TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:) endif !soil_pem #ifndef CPP_STD DO ig = 1,ngrid DO iloop = 1,nsoilmx tsoil(ig,iloop) = 0. inertiesoil(ig,iloop) = 0. DO islope = 1,nslope tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) & * subslope_dist(ig,islope) inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) & * subslope_dist(ig,islope) ENDDO ENDDO ENDDO ! of DO ig=1,ngrid ! III_a.3 Surface optical properties (for startfi) DO ig = 1,ngrid DO l = 1,2 albedo(ig,l) =0. DO islope = 1,nslope albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) & *subslope_dist(ig,islope) ENDDO ENDDO ENDDO DO ig = 1,ngrid emis(ig) =0. DO islope = 1,nslope emis(ig)= emis(ig)+emiss_slope(ig,islope) & *subslope_dist(ig,islope) ENDDO ENDDO DO ig = 1,ngrid tsurf(ig) = 0. DO islope = 1,nslope tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 & * subslope_dist(ig,islope) ENDDO tsurf(ig) = tsurf(ig)**(0.25)/emis(ig) ENDDO ! of DO ig=1,ngrid #endif ! III_a.4 Pressure (for start) do i=1,ip1jmp1 ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM enddo do i = 1,ngrid ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM enddo write(*,*) 'rapport ps',global_ave_press_new/global_ave_press_GCM ! III_a.5 Tracer (for start) allocate(zplev_new(ngrid,nlayer+1)) do l=1,nlayer+1 do ig=1,ngrid zplev_new(ig,l) = ap(l) + bp(l)*ps_phys(ig) enddo enddo DO nnq=1,nqtot if (noms(nnq).NE."co2") then DO l=1,llm-1 DO ig=1,ngrid q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1)) ENDDO q(:,llm,nnq)=q(:,llm-1,nnq) ENDDO else DO l=1,llm-1 DO ig=1,ngrid q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1)) & + ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - & (zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / & (zplev_new(ig,l)-zplev_new(ig,l+1)) ENDDO q(:,llm,nnq)=q(:,llm-1,nnq) ENDDO endif ENDDO ! Conserving the tracers's mass for GCM start files DO nnq=1,nqtot DO ig=1,ngrid DO l=1,llm-1 if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) q(ig,l,nnq)=1. q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2)) write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number" endif if(q(ig,l,nnq).LT.0) then q(ig,l,nnq)=1E-30 endif ENDDO enddo enddo !------------------------ if(evol_orbit_pem) then call recomp_orb_param(year_iter) endif ! III Output ! III_a Update surface value for the PCM start files ! III_b Write start and starfi.nc !------------------------ ! III_b.1 WRITE restart.nc ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys pday=day_ini ztime_fin=0. allocate(p(ip1jmp1,nlayer+1)) CALL pression (ip1jmp1,ap,bp,ps,p) CALL massdair(p,masse) CALL dynredem0("restart_evol.nc", day_ini, phis) CALL dynredem1("restart_evol.nc", & time_0,vcov,ucov,teta,q,masse,ps) print *, "restart_evol.nc has been written" ! III_b.2 WRITE restartfi.nc #ifndef CPP_STD call physdem0("restartfi_evol.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,0.,cell_area, & albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, & hmons,summit,base,nslope,def_slope, & subslope_dist) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & tsurf,tsoil,co2ice,albedo,emis, & q2,qsurf,tauscaling,totcloudfrac,wstar, & watercap,inertiesoil,nslope,co2ice_slope, & tsurf_slope,tsoil_slope, albedo_slope, & emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys) #else call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,time_phys,cell_area, & albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & tsurf,tsoil,emis,q2,qsurf, & cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) #endif print *, "restartfi_evol.nc has been written" !------------------------ ! III Output ! III_a Update surface value for the PCM start files ! III_b Write start and starfi.nc ! III_c Write start_pem !------------------------ call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , & tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys) print *, "restartfi_PEM.nc has been written" print *, "The PEM had run for ", year_iter, " years." print *, "LL & RV : So far so good" END PROGRAM pem