!------------------------ ! 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 ! 1: Modules needed for reading and writting startfi: use phyetat0_mod, only: phyetat0 use phyredem, only: physdem0, physdem1 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 ! netcdf is needed to read info like lat and lon in the physiq file use turb_mod, only: q2, wstar ! 1a: Modules specific from the Marsian physiq #ifndef CPP_STD use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil use surfdat_h, only: tsurf, emis,& qsurf,watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & hmons, summit, base,albedo_h2o_frost, & frost_albedo_threshold,emissiv,watercaptag,perenial_co2ice use dimradmars_mod, only: totcloudfrac, albedo use dust_param_mod, only: tauscaling use tracer_mod, only: noms,igcm_h2o_ice,igcm_co2,mmol,igcm_h2o_vap ! tracer names and molar masses #else ! 1b: Modules specific from the Generic physiq 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 ! tracer names use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, & ALBEDO_BAREGROUND,ALBEDO_CO2_ICE_SPECTV, phys_state_var_init use aerosol_mod, only : iniaerosol #endif ! 1c: Modules specific from the 1d marsian physiq #ifndef CPP_1D USE iniphysiq_mod, ONLY: iniphysiq USE control_mod, ONLY: iphysiq, day_step,nsplit_phys #else use time_phylmdz_mod, only: daysec, iphysiq,day_step #endif #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 planete_mod, only: apoastr, periastr, year_day, peri_day, & obliquit #endif USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & subslope_dist,iflat, & major_slope,ini_comslope_h #ifndef CPP_STD USE comcstfi_h, only: r, mugaz #else USE comcstfi_mod, only: r, mugaz #endif USE logic_mod, ONLY: iflag_phys USE mod_const_mpi, ONLY: COMM_LMDZ use time_phylmdz_mod, only: daysec,dtphys USE comconst_mod, ONLY: rad,g,cpp,pi USE infotrac USE geometry_mod, only: latitude_deg use conf_pem_mod, only: conf_pem use pemredem, only: pemdem0,pemdem1 use glaciers_mod, only: co2glaciers_evol,h2oglaciers_evol,co2glaciersflow,h2oglaciersflow use criterion_pem_stop_mod, only: criterion_waterice_stop,criterion_co2_stop use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, & m_co2,m_noco2,threshold_water_frost2perenial,threshold_co2_frost2perenial use evol_co2_ice_s_mod, only: evol_co2_ice_s use evol_h2o_ice_s_mod, only: evol_h2o_ice_s use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, & TI_PEM,inertiedat_PEM, & ! soil thermal inertia tsoil_PEM, mlayer_PEM,layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths fluxgeo, & ! geothermal flux for the PEM and GCM water_reservoir ! Water ressources use adsorption_mod, only : regolith_adsorption,adsorption_pem, & ! bool to check if adsorption, main subroutine ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays co2_adsorbded_phys, h2o_adsorbded_phys ! mass of co2 and h2O adsorbded 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 use ice_table_mod, only: porefillingice_depth,porefillingice_thickness,& end_ice_table_porefilling,ini_ice_table_porefilling, & computeice_table_equilibrium use soil_thermalproperties_mod, only: update_soil_thermalproperties #ifdef CPP_1D use regular_lonlat_mod, only: init_regular_lonlat use physics_distribution_mod, only: init_physics_distribution use mod_grid_phy_lmdz, only : regular_lonlat use init_phys_1d_mod, only : init_phys_1d #endif IMPLICIT NONE include "dimensions.h" include "paramet.h" INTEGER ngridmx PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) 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_start_GCM !(ngrid) pression au sol REAL, dimension(:,:),allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées 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 character*2 str2 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 :: ini_surf_h2o ! Initial surface of sublimating h2o ice REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice 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_pressure! Logical : is the criterion (% of change in the surface pressure) reached? INTEGER :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param real,save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] real ,allocatable :: vmr_co2_gcm(:,:) ! 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_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg] REAL ,allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] 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 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE REAL, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] REAL, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] REAL, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] REAL, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] REAL, dimension(:,:,:),allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] REAL, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice REAL, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state REAL, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state REAL, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year REAL, dimension(:,:), allocatable :: tendencies_co2_ice_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 ! physical pointx slope field : Tendency of evolution of perenial h2o ice REAL, dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope): Flag where there is a CO2 glacier flow REAL, dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) : Flag where there is a CO2 glacier flow REAL, dimension(:,:), allocatable :: flag_h2oflow(:,:) !(ngrid,nslope): Flag where there is a H2O glacier flow REAL, dimension(:), allocatable :: flag_h2oflow_mesh(:) !(ngrid) : Flag where there is a H2O glacier flow !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL REAL, ALLOCATABLE :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] REAL, ALLOCATABLE :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [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(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] 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 :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] REAL,ALLOCATABLE :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] REAL,ALLOCATABLE :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] REAL,ALLOCATABLE :: watersoil_density_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, ALLOCATABLE :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] REAL :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] REAL :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 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] REAL :: watercap_sum ! total mass of water cap [kg/m^2] REAL :: water_sum ! total mass of water in the mesh [kg/m^2] #ifdef CPP_STD 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 :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic REAL,allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic REAL,allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic REAL,allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic REAL,allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic REAL,allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case REAL,allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case REAL,allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case REAL,allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case REAL,allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model LOGICAL, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model REAL,allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case REAL,allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case #endif #ifdef CPP_1D INTEGER :: nsplit_phys integer,parameter:: jjm_value=jjm-1 integer :: ierr #else integer,parameter:: jjm_value=jjm #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 A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 year_day=669 daysec=88775. timestep=year_day*daysec/year_step !------------------------ ! I Initialisation ! I_a READ run.def !------------------------ !----------------------------READ run.def --------------------- #ifndef CPP_1D dtphys=0 CALL conf_gcm( 99, .TRUE. ) call infotrac_init nq=nqtot allocate(q(ip1jmp1,llm,nqtot)) #else ! load tracer names from file 'traceur.def' open(90,file='traceur.def',status='old',form='formatted',& iostat=ierr) if (ierr.ne.0) then write(*,*) 'Cannot find required file "traceur.def"' write(*,*) ' If you want to run with tracers, I need it' write(*,*) ' ... might as well stop here ...' stop else write(*,*) "pem1d: Reading file traceur.def" ! read number of tracers: read(90,*,iostat=ierr) nq print *, "nq",nq nqtot=nq ! set value of nqtot (in infotrac module) as nq if (ierr.ne.0) then write(*,*) "pem1d: error reading number of tracers" write(*,*) " (first line of traceur.def) " stop endif if (nq<1) then write(*,*) "pem1d: error number of tracers" write(*,*) "is nq=",nq," but must be >=1!" stop endif endif nq=nqtot allocate(q(ip1jmp1,llm,nqtot)) allocate(ap(nlayer+1)) allocate(bp(nlayer+1)) call init_phys_1d(llm,nqtot,vcov,ucov, & teta,q,ps, time_0,ap,bp) pi=2.E+0*asin(1.E+0) g=3.72 nsplit_phys=1 #endif CALL conf_pem !------------------------ ! 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 --------------------- !----------------------------READ start.nc --------------------- #ifndef CPP_1D CALL dynetat0(FILE_NAME_start,vcov,ucov, & teta,q,masse,ps,phis, time_0) #endif allocate(ps_start_GCM(ngrid)) #ifndef CPP_1D call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) 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) #else ps_start_GCM(1)=ps(1) #endif #ifndef CPP_1D 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) #endif ! 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) status = nf90_inq_varid(ncid, "soildepth", sdvarid) status = nf90_get_var(ncid, sdvarid, mlayer) status =nf90_close(ncid) !----------------------------READ startfi.nc --------------------- ! First we read the initial state (starfi.nc) #ifndef CPP_STD CALL phyetat0 (FILE_NAME,0,0, & nsoilmx,ngrid,nlayer,nq, & day_ini,time_phys, & tsurf,tsoil,albedo,emis, & q2,qsurf,tauscaling,totcloudfrac,wstar, & watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist) ! Remove unphysical values of surface tracer DO i=1,ngrid DO nnq=1,nqtot DO islope=1,nslope if(qsurf(i,nnq,islope).LT.0) then qsurf(i,nnq,islope)=0. endif enddo enddo enddo call surfini(ngrid,qsurf) #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() allocate(tsurf_read_generic(ngrid)) allocate(qsurf_read_generic(ngrid,nq)) allocate(tsoil_read_generic(ngrid,nsoilmx)) allocate(emis_read_generic(ngrid)) allocate(tsurf(ngrid,1)) allocate(qsurf(ngrid,nq,1)) allocate(tsoil(ngrid,nsoilmx,1)) allocate(emis(ngrid,1)) allocate(watercap(ngrid,1)) allocate(watercaptag(ngrid)) allocate(albedo_read_generic(ngrid,2)) allocate(albedo(ngrid,2,1)) allocate(inertiesoil(ngrid,nsoilmx,1)) call phyetat0(.true., & ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, & day_ini,time_phys,tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,qsurf_read_generic, & cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab, tsea_ice,sea_ice) call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) nslope=1 call ini_comslope_h(ngrid,1) qsurf(:,:,1)=qsurf_read_generic(:,:) tsurf(:,1)=tsurf_read_generic(:) tsoil(:,:,1)=tsoil_read_generic(:,:) emis(:,1)=emis_read_generic(:) watercap(:,1)=0. watercaptag(:)=.false. albedo(:,1,1)=albedo_read_generic(:,1) albedo(:,2,1)=albedo_read_generic(:,2) inertiesoil(:,:,1)=inertiedat(:,:) 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(i,nnq,1)=qsurf_read_generic(i,nnq) if(qsurf(i,nnq,1).LT.0) then qsurf(i,nnq,1)=0. endif enddo enddo #endif DO nnq=1,nqtot ! Why not using ini_tracer ? if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq if(noms(nnq).eq."h2o_vap") then igcm_h2o_vap = nnq mmol(igcm_h2o_vap)=18. endif if(noms(nnq).eq."co2") igcm_co2 = nnq ENDDO r= 8.314511E+0 *1000.E+0/mugaz !------------------------ ! 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)) allocate(flag_h2oflow(ngrid,nslope)) allocate(flag_h2oflow_mesh(ngrid)) flag_co2flow(:,:) = 0 flag_co2flow_mesh(:) = 0 flag_h2oflow(:,:) = 0 flag_h2oflow_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 !------------------------ ! 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(tsoil_ave(ngrid,nsoilmx,nslope)) allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope)) allocate(vmr_co2_gcm(ngrid,timelen)) allocate(ps_timeseries(ngrid,timelen)) allocate(min_co2_ice_1(ngrid,nslope)) allocate(min_h2o_ice_1(ngrid,nslope)) allocate(min_co2_ice_2(ngrid,nslope)) allocate(min_h2o_ice_2(ngrid,nslope)) allocate(tsurf_ave_yr1(ngrid,nslope)) allocate(tsurf_ave(ngrid,nslope)) allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen)) allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) allocate(co2_ice_GCM(ngrid,nslope,timelen)) allocate(watersurf_density_ave(ngrid,nslope)) allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(Tsurfave_before_saved(ngrid,nslope)) allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(delta_co2_adsorbded(ngrid)) allocate(delta_h2o_adsorbded(ngrid)) allocate(vmr_co2_pem_phys(ngrid,timelen)) print *, "Downloading data Y1..." call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,& tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & watersurf_density_ave,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" print *, "Downloading data Y2" call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, & tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & watersurf_density_ave,watersoil_density_timeseries) print *, "Downloading data Y2 done" !------------------------ ! 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) call end_adsorption_h_PEM call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) call end_ice_table_porefilling call ini_ice_table_porefilling(ngrid,nslope) if(soil_pem) then call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) DO l=1,nsoilmx tsoil_PEM(:,l,:)=tsoil_ave(:,l,:) tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:) watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:) ENDDO DO l=nsoilmx+1,nsoilmx_PEM tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:) watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:) ENDDO watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen endif !soil_pem deallocate(tsoil_ave) deallocate(tsoil_GCM_timeseries) !------------------------ ! 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_co2_ice(ngrid,nslope)) allocate(tendencies_co2_ice_ini(ngrid,nslope)) allocate(tendencies_h2o_ice(ngrid,nslope)) ! Compute the tendencies of the evolution of ice over the years call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,& min_co2_ice_2,tendencies_co2_ice) tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:) call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,& min_h2o_ice_2,tendencies_h2o_ice) deallocate(min_co2_ice_1) deallocate(min_co2_ice_2) deallocate(min_h2o_ice_1) deallocate(min_h2o_ice_2) !------------------------ ! 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_co2_ice_sublim(ngrid,nslope)) allocate(initial_co2_ice(ngrid,nslope)) allocate(initial_h2o_ice(ngrid,nslope)) ! We save the places where water ice is sublimating ! We compute the surface of water ice sublimating ini_surf_co2=0. ini_surf_h2o=0. Total_surface=0. do i=1,ngrid Total_surface=Total_surface+cell_area(i) do islope=1,nslope if (tendencies_co2_ice(i,islope).LT.0) then initial_co2_ice_sublim(i,islope)=1. ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) else initial_co2_ice_sublim(i,islope)=0. endif if (qsurf(i,igcm_co2,islope).GT.0) then initial_co2_ice(i,islope)=1. else initial_co2_ice(i,islope)=0. endif if (tendencies_h2o_ice(i,islope).LT.0) then initial_h2o_ice(i,islope)=1. ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) else initial_h2o_ice(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 (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_start_GCM(ig) ENDDO ENDDO global_ave_press_old=0. do i=1,ngrid global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(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("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth,porefillingice_thickness, & tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,& tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,& watersurf_density_ave,watersoil_density_PEM_ave, & co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) do ig = 1,ngrid do islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) enddo enddo if(adsorption_pem) then totmassco2_adsorbded = 0. totmassh2o_adsorbded = 0. do ig = 1,ngrid do islope =1, nslope do l = 1,nsoilmx_PEM - 1 totmassco2_adsorbded = totmassco2_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) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_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 of CO2 in the regolith=", totmassco2_adsorbded write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded endif ! adsorption deallocate(tsurf_ave_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(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 call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) if(adsorption_pem) then do i=1,ngrid global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface enddo print *, 'Global average pressure old time step',global_ave_press_old print *, 'Global average pressure new time step',global_ave_press_new endif ! 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_timeseries(ig,:) ENDDO ENDDO ! II.a.3. Surface pressure timeseries print *, "Recomputing the surface pressure timeserie adapted to the new pressure..." do ig = 1,ngrid ps_timeseries(ig,:) = ps_timeseries(ig,:)*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_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(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope) print *, "Evolution of co2 ice" call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) DO islope=1, nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) ENDDO !------------------------ ! II Run ! II_a update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 & H2O glaciers flows !------------------------ print *, "CO2 glacier flows" if (co2glaciersflow) then call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,& global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) endif print *, "H2O glacier flows" if (h2oglaciersflow) then call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh) endif DO islope=1, nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) ENDDO !------------------------ ! 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=.false. Tsurfave_before_saved(:,:) = tsurf_ave(:,:) DO ig = 1,ngrid DO islope = 1,nslope if(initial_co2_ice(ig,islope).gt.0.5 .and. qsurf(ig,igcm_co2,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(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) bool_sublim=.true. exit endif enddo if (bool_sublim.eqv. .true.) then exit endif enddo else do ig_loop=ig,1,-1 DO islope_loop=islope,iflat if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop) bool_sublim=.true. exit endif enddo if (bool_sublim.eqv. .true.) then exit endif enddo endif initial_co2_ice(ig,islope)=0 if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then albedo(ig,1,islope) = albedo_h2o_frost albedo(ig,2,islope) = albedo_h2o_frost else albedo(ig,1,islope) = albedodat(ig) albedo(ig,2,islope) = albedodat(ig) emis(ig,islope) = emissiv endif elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2 ave=0. do t=1,timelen if(co2_ice_GCM(ig,islope,t).gt.1e-3) then ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) else ave = ave + tsurf_GCM_timeseries(ig,islope,t) endif enddo tsurf_ave(ig,islope)=ave/timelen endif enddo enddo do t = 1,timelen tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:)) enddo ! for the start do ig = 1,ngrid do islope = 1,nslope tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope)) enddo enddo DO islope=1, nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) 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)) print *,"Updating soil temperature" ! Soil averaged do islope = 1,nslope TI_locslope(:,:) = TI_PEM(:,:,islope) do t = 1,timelen Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t) call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) do ig = 1,ngrid do isoil = 1,nsoilmx_PEM watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r) if(isnan(Tsoil_locslope(ig,isoil))) then call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1) endif enddo enddo enddo enddo tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen print *, "Update of soil temperature done" deallocate(TI_locslope) deallocate(Tsoil_locslope) deallocate(Tsurf_locslope) write(*,*) "Compute ice table" ! II_d.3 Update the ice table call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) print *, "Update soil propreties" ! II_d.4 Update the soil thermal properties call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, & porefillingice_depth,porefillingice_thickness,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbded if(adsorption_pem) then call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), & tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) totmassco2_adsorbded = 0. totmassh2o_adsorbded = 0. do ig = 1,ngrid do islope =1, nslope do l = 1,nsoilmx_PEM - 1 totmassco2_adsorbded = totmassco2_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) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_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 of CO2 in the regolith=", totmassco2_adsorbded write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded endif 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,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_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_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & initial_co2_ice_sublim,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 criterion_stop=1 endif if (STOPPING_1_water) then print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water criterion_stop=1 endif if (STOPPING_co2) then print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 criterion_stop=2 endif if (STOPPING_pressure) then print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure criterion_stop=3 endif if (year_iter.ge.year_iter_max) then print *, "STOPPING because maximum number of iterations reached" criterion_stop=4 endif if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) 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) ! H2O ice DO ig=1,ngrid if(watercaptag(ig)) then watercap_sum=0. DO islope=1,nslope if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost else ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap watercap(ig,islope)=watercap(ig,islope)+qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) qsurf(ig,igcm_h2o_ice,islope)=0. endif watercap_sum=watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) watercap(ig,islope)=0. enddo water_reservoir(ig)=water_reservoir(ig)+watercap_sum endif enddo DO ig=1,ngrid water_sum = 0. DO islope=1,nslope water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) ENDDO if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming. if(water_sum.gt.threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 watercaptag(ig)=.true. water_reservoir(ig)=water_reservoir(ig)+threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost DO islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.) ENDDO endif else ! let's check that the infinite source of water disapear if((water_sum + water_reservoir(ig)).lt.threshold_water_frost2perenial) then watercaptag(ig)=.false. DO islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) ENDDO water_reservoir(ig) = 0. endif endif enddo ! CO2 ice DO ig=1,ngrid DO islope=1,nslope if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) endif ENDDO ENDDO ! III_a.2 Tsoil update (for startfi) if(soil_pem) then call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) endif !soil_pem ! 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_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM enddo ! 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_start_GCM(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") .and. (noms(nnq).NE."stormdust_number") .and. (noms(nnq).NE."topdust_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)) #ifndef CPP_1D 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" #else DO nnq = 1, nqtot call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) ENDDO #endif ! 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,def_slope, & subslope_dist) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & tsurf,tsoil,inertiesoil,albedo, & emis,q2,qsurf,tauscaling,totcloudfrac,wstar,& watercap,perenial_co2ice) #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 pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & float(day_ini),0.,nslope,def_slope,subslope_dist) call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , & tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, & co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) call info_run_PEM(year_iter, criterion_stop) print *, "restartfi_PEM.nc has been written" print *, "The PEM had run for ", year_iter, " years." print *, "LL & RV : So far so good" deallocate(vmr_co2_gcm) deallocate(ps_timeseries) deallocate(tsurf_GCM_timeseries) deallocate(q_co2_PEM_phys) deallocate(q_h2o_PEM_phys) deallocate(co2_ice_GCM) deallocate(watersurf_density_ave) deallocate(watersoil_density_timeseries) deallocate(Tsurfave_before_saved) deallocate(tsoil_phys_PEM_timeseries) deallocate(watersoil_density_PEM_timeseries) deallocate(watersoil_density_PEM_ave) deallocate(delta_co2_adsorbded) deallocate(delta_h2o_adsorbded) deallocate(vmr_co2_pem_phys) END PROGRAM pem