!------------------------ ! I Initialisation ! I_a READ run.def , run_pem.def ! I_b READ starfi_0.nc and start_0.nc ! II Run ! II_a Compute tendencies ! II_b Save initial situation ! II_c Time loop ! III Output ! III_a Write startfi.nc !------------------------ PROGRAM pem !module needed for INITIALISATION use phyetat0_mod, only: phyetat0 use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer use surfdat_h, only: tsurf, co2ice, emis,& & qsurf,watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & & hmons, summit, base use dimradmars_mod, only: totcloudfrac, albedo, & ini_dimradmars_mod use turb_mod, only: q2, wstar, ini_turb_mod use dust_param_mod, only: tauscaling, ini_dust_param_mod 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,nqmx ! 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 temps_mod_evol, ONLY: nyear, dt_pem USE vertical_layers_mod, ONLY: ap,bp 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 IMPLICIT NONE include "dimensions.h" include "paramet.h" 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 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 ! 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 :: 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 ! Variable for initial situation of water ice 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 ! Variable for initial situation of co2 ice 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 ! Variable for needed to compute tendencies for water 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 :: h2o_ice_first_last_day ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice ! Variable for needed to compute tendencies for co2 ice 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 , dimension(:,:,:), allocatable :: co2_ice_first_last_day ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice REAL :: global_ave_press_GCM ! physical point field : Global average pressure in the GCM output 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 ! ngrid x nlayer+1 : Pressure level at the end of the iteration REAL , dimension(:,:), allocatable :: ZPLEV_OLD ! ngrid x nlayer+1 : Pressure level at the beginning of the iteration INTEGER :: year_iter !Counter for the number of PEM iteration INTEGER :: year ! Year read in the startfi files and written back in it counting the year iterations 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 co2 ice) reached? LOGICAL :: STOPPING_1_co2 ! Logical : is there still co2 ice to sublimate? REAL, dimension(:),allocatable :: q_co2_GCM ! Initial amount of co2 in the first layer REAL ps_GCM(ip1jmp1) ! pression au sol donné par le GCM real ,allocatable :: vmr_co2_gcm_phys(:,:) !(ngrid x Timelen) : co2 volume mixing ratio given by the GCM REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:) ! iim+1 x jjm+1 x Timelen : co2 volume mixing ratio given by the GCM REAL, ALLOCATABLE :: ps_GCM_1(:,:,:) ! iim+1 x jjm+1 x Timelen : Surface pressure in the first GCM year REAL, ALLOCATABLE :: ps_GCM_2(:,:,:) ! iim+1 x jjm+1 x Timelen : Surface pressure in the second GCM year integer :: timelen ! Number of time point in the GCM diagfi output REAL, ALLOCATABLE :: p(:,:) !(ngrid,llmp1) : Pressure !!!!!!!!!!!!!!!!!!!!!!!! SLOPE REAL ,allocatable :: watercap_slope(:,:) ! (ngrid,nslope) REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of co2 ice at each point for the first year REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of co2 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 :: initial_co2_ice_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 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_h2o_ice_phys_slope ! physical point field : Tendency of evolution of perenial co2 ice over a year REAL, PARAMETER :: co2_hmax = 10 ! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Loop variable INTEGER :: i,j,ig0,l,ig,nnq,t,islope ! Constant REAL :: pi,beta,alpha ! 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 !------------------------ !======================================================= ! I Initialisation ! I_a READ run.def !------------------------ !----------------------------READ run.def --------------------- CALL conf_gcm( 99, .TRUE. ) CALL conf_evol( 99, .TRUE. ) !------------------------ ! I Initialisation ! I_a READ run.def ! I_b READ starfi_0.nc !----------------------------READ startfi.nc --------------------- !----------------------------initialisation --------------------- status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) status =nf90_inq_dimid(ncid, "physical_points", phydimid) status =nf90_inquire_dimension(ncid, phydimid, len = ngrid) status =nf90_inq_dimid(ncid, "nlayer", nlayerdimid) status =nf90_inquire_dimension(ncid, nlayerdimid, len = nlayer) status =nf90_inq_dimid(ncid,"number_of_advected_fields",nqdimid) status =nf90_inquire_dimension(ncid, nqdimid, len = nq) 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 --------------------- call infotrac_init allocate(q(ip1jmp1,llm,nqtot)) CALL dynetat0(FILE_NAME_start,vcov,ucov, & teta,q,masse,ps,phis, time_0) CALL iniconst CALL inigeom 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) year=day_ini !----------------------------reading --------------------- ! First we read the initial state (starfi.nc) allocate(watercap_slope(ngrid,nslope)) 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,nslope,tsurf_slope, & tsoil_slope,co2ice_slope,def_slope,def_slope_mean, & subslope_dist,albedo_slope,emiss_slope, & qsurf_slope,watercap_slope) ! Slope : Definition of the flat slope iflat=1 DO islope=2,nslope IF(abs(def_slope_mean(islope)).lt. & abs(def_slope_mean(iflat)))THEN iflat = islope ENDIF ENDDO PRINT*,'Flat slope for islope = ',iflat PRINT*,'corresponding criterium = ',def_slope_mean(iflat) allocate(flag_co2flow(ngrid,nslope)) allocate(flag_co2flow_mesh(ngrid)) flag_co2flow(:,:) = 0. flag_co2flow_mesh(:) = 0. ! Then we read the evolution of water and co2 ice over the first year of the GCM run, saving only the minimum value allocate(min_h2o_ice_s_1(iim+1,jjm+1)) allocate(min_co2_ice_s_1(iim+1,jjm+1)) allocate(ps_GCM_1(iim+1,jjm+1,2676)) allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope)) allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope)) allocate(vmr_co2_gcm(iim+1,jjm+1,2676)) call pemetat0 ("concat_year_one.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,nslope) ! Then we read the evolution of water and co2 ice over the second year of the GCM run, saving only the minimum value 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)) allocate(ps_GCM_2(iim+1,jjm+1,2676)) call pemetat0 ("concat_year_two.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_2,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2,nslope) ! We read the evolutaion of water ice over the first year of the GCM run, saving the first and last state allocate(vmr_co2_gcm_phys(ngrid,timelen)) vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:) ig0 = 2 DO j=2,jjm DO i = 1, iim vmr_co2_gcm_phys(ig0,:) =vmr_co2_gcm(i,j,:) ig0= ig0 + 1 ENDDO ENDDO vmr_co2_gcm_phys(ig0,:) = vmr_co2_gcm(1,jjm+1,:) !======================================================= ! II Run ! II_a Compute tendencies !------------------------ !---------------------------- RUN --------------------- !----- Compute tendencies 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_h2o_ice_slope(iim+1,jjm+1,nslope)) allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope)) ! Compute the tendencies of the evolution of water 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) ! Compute the tendencies of the evolution of co2 ice over the years call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,& min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys) ! Compute the tendencies of the evolution of water ice over the years with slopes 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) ! Compute the tendencies of the evolution of co2 ice over the years with slopes 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) !------------------------ ! II Run ! II_a Compute tendencies ! II_b Save initial situation !------------------------ !----- Save initial situation allocate(initial_h2o_ice(ngrid)) allocate(initial_co2_ice(ngrid)) allocate(q_co2_GCM(ngrid)) allocate(initial_co2_ice_slope(ngrid,nslope)) allocate(initial_h2o_ice_slope(ngrid,nslope)) ! We save the places where water/co2 ice is sublimating do i=1,ngrid if (tendencies_h2o_ice_phys(i).LT.0) then initial_h2o_ice(i)=1. 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_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. else initial_h2o_ice_slope(i,islope)=0. endif enddo enddo ! 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 if (initial_h2o_ice(i).GT.0.5) then ini_surf=ini_surf+cell_area(i) endif do islope=1,nslope if (initial_co2_ice_slope(i,islope).GT.0.5) then ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) endif if (initial_h2o_ice_slope(i,islope).GT.0.5) then ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) endif enddo Total_surface=Total_surface+cell_area(i) enddo ! We put the pressure on the physical grid allocate(ps_phys(ngrid)) ig0 = iim+1 ps_phys(1)=ps(ig0) ig0=ig0+1 DO i = 2, ngrid-1 if(modulo(ig0,iim+1).eq.0 .and. i.gt.3) then ig0=ig0+1 endif ps_phys(i) = ps(ig0) ig0= ig0 + 1 ENDDO ps_phys(ngrid)=ps(ig0) ! We compute the global average pressure of the initial state 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 ! We compute the initial pressure level allocate(zplev_new(ngrid,nlayer+1)) allocate(zplev_old(ngrid,nlayer+1)) DO l=1,nlayer+1 DO ig=1,ip1jmp1 zplev_old(ig,l) = ap(l) + bp(l)*ps(ig) ENDDO ENDDO !======================================================= !------------------------ ! II Run ! II_a Compute tendencies ! II_b Save initial situation ! II_c Time loop !------------------------ !----- Time loop year_iter=0 do while (.true.) ! Compute the new global average pressure given the tendencies of co2 condensation global_ave_press_new=global_ave_press_old do i=1,ngrid do islope=1,nslope global_ave_press_new=global_ave_press_new-dt_pem*g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface enddo enddo ! Compute the new surface pressure at each physical point do i=1,ngrid ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old enddo ! Compute the new surface pressure at each dynamical point do i=1,ip1jmp1 ps(i)=ps(i)*global_ave_press_new/global_ave_press_old enddo ! Compute the new pressure level at each physical point DO l=1,nlayer+1 DO ig=1,ip1jmp1 zplev_new(ig,l) = ap(l) + bp(l)*ps(ig) ENDDO ENDDO ! Compute the new values of tracer given the change of pressure and pressure level DO nnq=1,nqtot ! For each tracer if (nnq.NE.1) then ! If it is not co2 DO l=1,llm-1 ! For each level DO i=1,ip1jmp1 ! Each point q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1)) ENDDO ! i loop q(:,llm,nnq)=q(:,llm-1,nnq) ENDDO ! l loop else ! If it is co2 DO l=1,llm-1 ! For each level DO i=1,ip1jmp1 ! Each point q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1)) & + ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - & (zplev_old(ig,l)-zplev_old(ig,l+1)) ) / & (zplev_new(ig,l)-zplev_new(ig,l+1)) ENDDO ! i loop q(:,llm,nnq)=q(:,llm-1,nnq) ENDDO ! l loop endif ENDDO ! nnq loop !We apply the tendency on the ice call evol_h2o_ice_s_slope(qsurf_slope(:,6,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope) call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope) ! FLOW : Make the co2 ice 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 * & 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 * & 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 * & cos(pi*def_slope_mean(islope)/180.) flag_co2flow(ig,islope) = 1. flag_co2flow_mesh(ig) = 1. ENDIF ! co2ice > hmax ENDIF ! iflat ENDDO !islope ENDDO !ig ! End of FLOW ! We recompute the tendencies of co2 given the new surface pressure call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,vmr_co2_gcm_phys,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) year_iter=year_iter+dt_pem ! We check if we should stop : the criterion is based on the area of ice that is sublimating call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,6,:),STOPPING_water,ngrid,initial_h2o_ice) call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_slope,global_ave_press_GCM,global_ave_press_new,nslope) print *, "STOPPING_water", STOPPING_water, "STOPPING1_water", STOPPING_1_water print *, "STOPPING_co2", STOPPING_co2, "STOPPING1_co2", STOPPING_1_co2 print *, "year_iter=", year_iter year=year+dt_pem ! We check if the orbits parameters have changed too much call orbit_param(year) if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then exit endif ! Preparing for next iteration global_ave_press_old=global_ave_press_new zplev_old=zplev_new enddo ! We put the correct thickness of co2ice given the slope 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 ! We put the correct thickness of water ice given the slope DO ig = 1,ngrid qsurf(ig,6) = 0. DO islope = 1,nslope qsurf(ig,6) = qsurf(ig,6) + qsurf_slope(ig,6,islope) & * subslope_dist(ig,islope) / & cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO ! of DO ig=1,ngrid !------------------------ ! III Output ! III_a Write startfi.nc !------------------------ !----------------------------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) day_ini=year CALL dynredem0("restart_evol.nc", day_ini, phis) CALL dynredem1("restart_evol.nc", & time_0,vcov,ucov,teta,q,masse,ps) !----------------------------WRITE restartfi.nc --------------------- pday=year 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,nslope,co2ice_slope, & tsurf_slope,tsoil_slope, albedo_slope, & emiss_slope,qsurf_slope,watercap_slope) print *, "RV : So far so good" deallocate(longitude) deallocate(latitude) deallocate(cell_area) END PROGRAM pem