!------------------------ ! 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 use phyetat0_mod, only: phyetat0 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 turb_mod, only: q2, wstar use dust_param_mod, only: tauscaling 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 use phyredem, only: physdem0, physdem1 use tracer_mod, only: noms,igcm_h2o_ice ! tracer names ! For phyredem : USE control_mod, ONLY: iphysiq, day_step,nsplit_phys use time_phylmdz_mod, only: daysec use mod_phys_lmdz_para, only: is_parallel, is_sequential, & is_mpi_root, is_omp_root, & is_master use time_phylmdz_mod, only: dtphys USE mod_const_mpi, ONLY: COMM_LMDZ USE comconst_mod, ONLY: rad,g,r,cpp USE logic_mod, ONLY: iflag_phys USE iniphysiq_mod, ONLY: iniphysiq USE infotrac USE comcstfi_h, only: pi 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 USE geometry_mod, only: latitude_deg use pemredem, only: pemdem1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL USE soil_evolution_mod, ONLY: soil_pem, soil_pem_CN use comsoil_h_PEM, only: 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 planete_h, only: aphelie, periheli, year_day, peri_day, & obliquit 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 REAL :: global_ave_press_old ! physical point field : Global average pressure of initial/previous time step REAL :: global_ave_press_new ! physical point field : Global average pressure of current time step REAL , dimension(:,:), allocatable :: zplev_new REAL , dimension(:,:), allocatable :: zplev_gcm REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! same but with the time series REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series 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 ! Initial amount of co2 in the first layer real,save :: m_co2, m_noco2, A , B, mmean real ,allocatable :: vmr_co2_gcm_phys(:,:) !(ngrid) ! co2 volume mixing ratio real ,allocatable :: vmr_co2_pem_phys(:,:) !(ngrid) ! co2 volume mixing ratio real ,allocatable :: q_h2o_GCM_phys(:,:) real ,allocatable :: q_co2_GCM_phys(:,:) real ,allocatable :: q_co2_PEM_phys(:,:) REAL, ALLOCATABLE :: ps_GCM(:,:,:) REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:) REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:) REAL, ALLOCATABLE :: q_h2o_GCM(:,:,:) REAL ,allocatable :: q_h2o_PEM_phys(:,:) integer :: timelen REAL :: ave REAL, ALLOCATABLE :: p(:,:) !(ngrid,llmp1) REAL :: extra_mass ! Extra mass of a tracer if it is greater than 1 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE character*2 :: str2 REAL ,allocatable :: watercap_slope(:,:) !(ngrid,nslope) REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM 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 sublimating point of h2o ice REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating sublimating point of co2 ice REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope ! physical point field : Tendency of evolution of perenial co2 ice over a year REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope_ini ! physical point 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 point field : Tendency of evolution of perenial co2 ice over a year REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: co2_hmax ! Maximum height for CO2 deposit on slopes (m) REAL, PARAMETER :: rho_co2 = 1600 ! CO2 ice density (kg/m^3) INTEGER :: iaval ! Index of the neighboord slope () 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(:,:) ! IG 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 [K] REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K] REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia [SI] REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K] REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG 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(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] REAL, PARAMETER :: year_step = 1 INTEGER :: year_iter !Counter for the number of PEM iteration INTEGER :: year_iter_max REAL :: timestep REAL, ALLOCATABLE :: inertiesoil(:,:) ! Thermal inertia of the mesh for restart REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Averaged GCM Thermal Inertia [SI] REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia [SI] REAL,ALLOCATABLE :: ice_depth(:,:) REAL,ALLOCATABLE :: TI_locslope(:,:) REAL,ALLOCATABLE :: Tsoil_locslope(:,:) REAL,ALLOCATABLE :: Tsurf_locslope(:) REAL,ALLOCATABLE :: alph_locslope(:,:) REAL,ALLOCATABLE :: beta_locslope(:,:) REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) REAL, ALLOCATABLE :: delta_co2_adsorbded(:) REAL :: totmass_adsorbded !!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Loop variable LOGICAL :: bool_sublim INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop REAL :: beta = 3182.48 REAL :: alpha = 23.3494 ! Parallel variables is_sequential=.true. is_parallel=.false. is_mpi_root=.true. is_omp_root=.true. is_master=.true. 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) DO nnq=1,nqtot if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq ENDDO !----------------------------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)) 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,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 !------------------------ ! 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) ! CO2 max thickness (for glaciers flows) allocate(co2_hmax(nslope)) if(nslope.eq.7) then ! ugly way to implement that ... ! CF documentation that explain how the values are computed co2_hmax(1) = 1.5 co2_hmax(7) = co2_hmax(1) co2_hmax(2) = 2.4 co2_hmax(6) = co2_hmax(2) co2_hmax(3) = 5.6 co2_hmax(5) = co2_hmax(3) co2_hmax(4) = 1000000. endif 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)) print *, "Downloading data Y1..." call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,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) ! 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", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,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) 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)) 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) if(soil_pem) then deallocate(tsurf_GCM_timeseries) call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM) DO islope = 1,nslope 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)) 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)) ENDDO ENDDO 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,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,& tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded) 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 CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) 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 print *, "Max_iter_pem", Max_iter_pem 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 if(soil_pem) then global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface endif enddo ! 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 flow" DO ig = 1,ngrid DO islope = 1,nslope IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground ! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate on flat ground IF(co2ice_slope(ig,islope).ge.rho_co2*co2_hmax(islope) * & cos(pi*def_slope_mean(islope)/180.)) THEN ! Second: determine the flatest slopes possible: IF(islope.gt.iflat) THEN iaval=islope-1 ELSE iaval=islope+1 ENDIF do while ((iaval.ne.iflat).and. & (subslope_dist(ig,iaval).eq.0)) IF(iaval.gt.iflat) THEN iaval=iaval-1 ELSE iaval=iaval+1 ENDIF enddo co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + & (co2ice_slope(ig,islope) - rho_co2* co2_hmax(islope) * & cos(pi*def_slope_mean(islope)/180.)) * & subslope_dist(ig,islope)/subslope_dist(ig,iaval) * & cos(pi*def_slope_mean(iaval)/180.) / & cos(pi*def_slope_mean(islope)/180.) co2ice_slope(ig,islope)=rho_co2*co2_hmax(islope) * & cos(pi*def_slope_mean(islope)/180.) flag_co2flow(ig,islope) = 1. flag_co2flow_mesh(ig) = 1. ENDIF ! co2ice > hmax ENDIF ! iflattsoil_lope ENDDO !islope ENDDO !ig !------------------------ ! 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-5) 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-5) 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-5) 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-5).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-5) THEN !Put tsurf as tcond co2 ave=0. do t=1,timelen if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then ave = ave + beta/(alpha-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) Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope) Tsurf_locslope(:) = tsurf_ave_phys(:,islope) alph_locslope(:,:) = alph_PEM(:,:,islope) beta_locslope(:,:) = beta_PEM(:,:,islope) call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) do t = 1,timelen tsoil_phys_PEM_timeseries(:,:,islope,t) = tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope)) enddo TI_PEM(:,:,islope) = TI_locslope(:,:) tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:) tsurf_ave_phys(:,islope) = Tsurf_locslope(:) alph_PEM(:,:,islope) = alph_locslope(:,:) beta_PEM(:,:,islope) = beta_locslope(:,:) enddo print *, "Update of soil temperature done" ! Check Nan in the time series do ig = 1,ngrid do islope = 1,nslope do iloop = 1,nsoilmx_PEM if(isnan(tsoil_PEM(ig,iloop,islope))) then write(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope) stop endif enddo enddo enddo write(*,*) "Compute ice table" ! II_d.3 Update the ice table call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & ps_phys_timeseries, 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 (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 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 Tsurf and 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 DO ig = 1,ngrid tsurf(ig) = 0. DO islope = 1,nslope tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) & * subslope_dist(ig,islope) ENDDO ENDDO ! of DO ig=1,ngrid 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 ! 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 ! 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" .or. 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)) 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 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) 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