source: trunk/LMDZ.COMMON/libf/evolution/pem.F90 @ 4071

Last change on this file since 4071 was 4071, checked in by jbclement, 3 weeks ago

PEM:

  • Making the computation of ice tendencies more reliable by doing it after 'read_startpem' to know exactly where perennial ice is (no matter if there is a "startpem.nc" or not). Moreover, when there is no "startpem.nc", location of perennial ice depends now on 'watercaptag' and on huge amount of frost. This prevents negative ice tendency while there is no ice which can happen with weird PCM inputs (i.e. 'watercaptag' = F & 'watercap' /= 0 & no "stratpem.nc").
  • Few small safeguards throughout the code.

JBC
File size: 29.5 KB
Line 
1PROGRAM pem
2!-----------------------------------------------------------------------
3! NAME
4!     pem
5!
6! DESCRIPTION
7!     Main entry point for the Planetary Evolution Model (PEM).
8!
9! AUTHORS & DATE
10!     R. Vandemeulebrouck, 22/07/2022
11!     L. Lange, 22/07/2022
12!     JB Clement, 2023-2025
13!
14! NOTES
15!
16!-----------------------------------------------------------------------
17
18! DEPENDENCIES
19! ------------
20! Common modules
21use job_timelimit_mod, only: timelimit, antetime, timewall
22use parse_args_mod,    only: parse_args
23! PEM modules
24use allocation,         only: ini_allocation, end_allocation
25use atmosphere,         only: ps_PCM, evolve_pressure, CO2cond_ps_PCM, build4PCM_atmosphere
26use clim_state_init,    only: read_start, read_startfi, read_startpem
27use clim_state_rec,     only: write_restart, write_restartfi, write_restartpem
28use config,             only: read_rundef
29use display,            only: print_ini, print_end, print_msg
30use evolution,          only: n_yr_run, n_yr_sim, nmax_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
31use geometry,           only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface, nlayer
32use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
33use ice_table,          only: icetable_equilibrium, icetable_dynamic, icetable_depth, icetable_thickness, ice_porefilling, evolve_ice_table
34use info,               only: read_info, update_info
35use layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering
36use maths,              only: pi
37use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
38use orbit,              only: evol_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
39use output,             only: write_diagevol, dim_ngrid, dim_nsoil
40use physics,            only: g
41use slopes,             only: subslope_dist, def_slope_mean
42use soil,               only: do_soil, set_soil, TI, build4PCM_soil
43use soil_temp,          only: tsoil_PCM, shift_tsoil2surf, evolve_soil_temp
44use soil_therm_inertia, only: update_soil_TI
45use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
46use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice
47use surface,            only: emissivity_PCM, build4PCM_surf_rad_prop
48use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice
49use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf
50use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
51use tracers,            only: adapt_tracers2pressure, build4PCM_tracers, nq
52use utility,            only: real2str
53use xios_data,          only: load_xios_data
54
55! DECLARATION
56! -----------
57implicit none
58
59! LOCAL VARIABLES
60! ---------------
61! Utility-related:
62integer(li)               :: cr     ! Number of clock ticks per second (count rate)
63integer(li)               :: c1, c2 ! Counts of processor clock
64character(:), allocatable :: num    ! To write slope variables
65integer(di)               :: i, islope
66! Pressure-related:
67real(dp), dimension(:),   allocatable :: ps_avg          ! Average surface pressure [Pa]
68real(dp), dimension(:),   allocatable :: ps_dev          ! Deviation of surface pressure [Pa]
69real(dp), dimension(:,:), allocatable :: ps_ts           ! Surface pressure timeseries [Pa]
70real(dp)                              :: ps_avg_glob_ini ! Global average pressure at initialization [Pa]
71real(dp)                              :: ps_avg_glob_old ! Global average pressure of previous time step [Pa]
72real(dp)                              :: ps_avg_glob     ! Global average pressure of current time step [Pa]
73real(dp), dimension(:),   allocatable :: ps4PCM          ! Surface pressure reconstruction to feed back into PCM [Pa]
74real(dp), dimension(:,:), allocatable :: teta4PCM        ! Potential temperature reconstruction to feed back into PCM [K]
75real(dp), dimension(:,:), allocatable :: air_mass4PCM    ! Air mass reconstruction to feed back into PCM [kg]
76real(dp)                              :: pa4PCM          ! Transition pressure (for hybrid coord.) [Pa]
77real(dp)                              :: preff4PCM       ! Reference pressure [Pa]
78! Ice-related:
79real(dp),    dimension(:,:),   allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
80real(dp),    dimension(:,:),   allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
81real(dp)                                   :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
82real(dp)                                   :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
83logical(k4), dimension(:,:),   allocatable :: is_h2oice_ini              ! Initial location of H2O ice
84logical(k4), dimension(:,:),   allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
85logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
86logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
87real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice              ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
88real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice              ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
89real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost               ! Minimum of H2O frost over the last PCM year [kg.m-2]
90real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost               ! Minimum of CO2 frost over the last PCM year [kg.m-2]
91real(dp),    dimension(:,:),   allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
92real(dp),    dimension(:,:),   allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
93logical(k4), dimension(:),     allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
94logical(k4), dimension(:,:),   allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
95! Surface-related:
96real(dp), dimension(:,:), allocatable :: tsurf_avg           ! Average surface temperature [K]
97real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1       ! Average surface temperature of the second to last PCM run [K]
98real(dp), dimension(:,:), allocatable :: tsurf_dev           ! Deviation of surface temperature [K]
99real(dp), dimension(:,:), allocatable :: h2o_surfdensity_avg ! Average water surface density [kg/m^3]
100real(dp), dimension(:,:), allocatable :: zshift_surf         ! Elevation shift for the surface [m]
101real(dp), dimension(:,:), allocatable :: zlag                ! Newly built lag thickness [m]
102real(dp), dimension(:,:), allocatable :: albedo4PCM          ! Albedo reconstruction to feed back into PCM
103real(dp), dimension(:,:), allocatable :: emissivity4PCM      ! Emissivity reconstruction to feed back into PCM
104real(dp), dimension(:,:), allocatable :: tsurf4PCM           ! Surface temperature reconstruction to feed back into PCM [K]
105! Soil-related:
106real(dp), dimension(:,:,:),   allocatable :: tsoil_avg                ! Average soil temperature [K]
107real(dp), dimension(:,:,:),   allocatable :: tsoil_dev                ! Deviation pf soil temperature [K]
108real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts                 ! Soil temperature timeseries [K]
109real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts_old             ! Soil temperature timeseries at the previous time step [K]
110real(dp), dimension(:,:,:),   allocatable :: h2o_soildensity_avg      ! Average of soil water soil density [kg/m^3]
111real(dp), dimension(:),       allocatable :: delta_co2_ads            ! Quantity of CO2 exchanged due to adsorption/desorption [kg/m^2]
112real(dp), dimension(:),       allocatable :: delta_h2o_ads            ! Quantity of H2O exchanged due to adsorption/desorption [kg/m^2]
113real(qp)                                  :: totmass_adsh2o           ! Total mass of H2O exchanged because of adsorption/desorption [kg]
114real(dp), dimension(:,:,:),   allocatable :: inertiesoil4PCM          ! Soil thermal inertia reconstruction to feed back into PCM [SI]
115real(dp), dimension(:,:,:),   allocatable :: tsoil4PCM                ! Average soil temperature reconstruction to feed back into PCM [K]
116real(dp), dimension(:,:),     allocatable :: flux_geo4PCM             ! Geothermal flux reconstruction to feed back into PCM [W/m2]
117real(dp), dimension(:,:,:),   allocatable :: h2o_ads_reg              ! H2O adsorbed in the regolith [kg/m^2]
118real(dp), dimension(:,:,:),   allocatable :: co2_ads_reg              ! CO2 adsorbed in the regolith [kg/m^2]
119real(dp), dimension(:,:),     allocatable :: icetable_depth_old       ! Old depth of the ice table [m]
120real(dp), dimension(:),       allocatable :: delta_icetable ! Total mass of the H2O that has sublimated / condenses from the ice table [kg]
121! Tracer-related:
122real(dp), dimension(:,:),   allocatable :: q_co2_ts     ! CO2 mass mixing ratio in the first layer [kg/kg]
123real(dp), dimension(:,:),   allocatable :: q_co2_ts_ini ! Initial CO2 mass mixing ratio in the first layer [kg/kg]
124real(dp), dimension(:,:),   allocatable :: q_h2o_ts     ! H2O mass mixing ratio in the first layer [kg/kg]
125real(dp), dimension(:,:,:), allocatable :: q4PCM        ! Tracers reconstruction to feed back into PCM [kg/kg]
126! Tendency-related:
127real(dp), dimension(:,:), allocatable :: d_co2ice     ! Tendency of perennial CO2 ice [kg/m2/y]
128real(dp), dimension(:,:), allocatable :: d_co2ice_ini ! Tendency of perennial CO2 ice at the beginning [kg/m2/y]
129real(dp), dimension(:,:), allocatable :: d_h2oice     ! Tendency of perennial H2O ice [kg/m2/y]
130real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y]
131! Layering-related:
132type(layering), dimension(:,:), allocatable :: layerings_map    ! Layering for each grid point and slope
133type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
134real(dp),       dimension(:,:), allocatable :: h2oice_depth     ! Depth of subsurface ice layer
135real(dp),       dimension(:,:), allocatable :: h2oice_depth_old ! Old depth of subsurface ice layer
136logical(k4),    dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm
137! Evolution-related:
138real(dp)        :: nmax_yr_runorb ! Maximum number of years for the run due to orbital parameters
139type(stopFlags) :: stopcrit       ! Stopping criteria
140! Balance-related
141real(qp) :: totmass_co2ice, totmass_atmco2, totmass_adsco2             ! Current total CO2 masses (surface ice|atmospheric|adsorbed)
142real(qp) :: totmass_co2ice_ini, totmass_atmco2_ini, totmass_adsco2_ini ! Initial total CO2 masses (surface ice|atmospheric|adsorbed)
143real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
144real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm   ! Variables to balance H2O ice reservoirs
145
146! CODE
147! ----
148! Pre-processing step
149!~~~~~~~~~~~~~~~~~~~~
150! Elapsed time with system clock
151call system_clock(count_rate = cr)
152call system_clock(c1)
153
154! Parse command-line options
155call parse_args()
156
157! Initialization
158! ~~~~~~~~~~~~~~
159! Header
160call print_ini()
161
162! Allocate module arrays
163call ini_allocation()
164
165! Read the duration information of the simulation
166call read_info()
167
168! Read the PEM parameters
169call read_rundef()
170
171! Read the orbital parameters
172call read_orbitpm(n_yr_sim)
173
174! Read the "start.nc"
175call read_start()
176
177! Read the "startfi.nc"
178call read_startfi()
179
180! Read the PCM data given by XIOS
181allocate(ps_avg(ngrid),ps_ts(ngrid,nday))
182allocate(tsurf_avg(ngrid,nslope),tsurf_avg_yr1(ngrid,nslope),h2o_surfdensity_avg(ngrid,nslope))
183allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope))
184allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday))
185allocate(minPCM_h2operice(ngrid,nslope,2),minPCM_co2perice(ngrid,nslope,2),minPCM_h2ofrost(ngrid,nslope,2),minPCM_co2frost(ngrid,nslope,2))
186
187call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
188                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
189
190! Initiate soil settings and TI
191if (do_soil) call set_soil(TI)
192
193! Compute the deviation from the average
194allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoil_PCM,nslope))
195ps_dev(:) = ps_PCM(:) - ps_avg(:)
196tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:)
197tsoil_dev(:,:,:) = tsoil_PCM(:,:,:) - tsoil_avg(:,1:nsoil_PCM,:)
198
199! Compute global surface pressure
200ps_avg_glob = sum(cell_area*ps_avg)/total_surface
201
202! Compute the accepted maximum number of years due to orbital parameters (if needed)
203call compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)
204
205! Read the "startpem.nc"
206allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,nslope))
207allocate(h2o_ads_reg(ngrid,nsoil,nslope),co2_ads_reg(ngrid,nsoil,nslope),delta_h2o_ads(ngrid),delta_co2_ads(ngrid))
208allocate(layerings_map(ngrid,nslope))
209call read_startpem(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
210                   tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map,                         &
211                   h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
212deallocate(tsurf_avg_yr1)
213
214! Compute ice tendencies from yearly minima
215allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
216call print_msg('> Computing surface ice tendencies')
217call compute_tendice(minPCM_h2operice + minPCM_h2ofrost,h2o_ice(:,:) > 0._dp,d_h2oice)
218call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)))
219call compute_tendice(minPCM_co2perice + minPCM_co2frost,co2_ice(:,:) > 0._dp,d_co2ice)
220call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
221deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
222
223! Save initial set-up useful for the next computations
224call print_msg('> Saving some initial climate state variables')
225allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday))
226allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
227ps_avg_glob_ini = ps_avg_glob
228d_co2ice_ini(:,:) = d_co2ice(:,:)
229q_co2_ts_ini(:,:) = q_co2_ts(:,:)
230h2oice_sublim_coverage_ini = 0._dp
231co2ice_sublim_coverage_ini = 0._dp
232is_h2oice_ini(:,:) = .false.
233is_co2ice_ini(:,:) = .false.
234totmass_co2ice_ini = 0._qp
235totmass_atmco2_ini = 0._qp
236totmass_adsco2_ini = 0._qp
237totmass_adsh2o = 0._qp
238do i = 1,ngrid
239    totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
240    do islope = 1,nslope
241        totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
242        if (co2_ice(i,islope) > 0._dp) then
243            is_co2ice_ini(i,islope) = .true.
244            if (d_co2ice(i,islope) < 0._dp) co2ice_sublim_coverage_ini = co2ice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
245        end if
246        if (h2o_ice(i,islope) > 0._dp) then
247            is_h2oice_ini(i,islope) = .true.
248            if (d_co2ice(i,islope) < 0._dp) h2oice_sublim_coverage_ini = h2oice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
249        end if
250    end do
251end do
252call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini))
253call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini))
254call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini))
255if (do_sorption) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2_ini,totmass_adsh2o)
256
257! Main evolution loop
258! ~~~~~~~~~~~~~~~~~~~
259call print_msg('')
260call print_msg('********* Evolution *********')
261if (do_layering) then
262    allocate(h2oice_depth(ngrid,nslope),h2oice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
263    new_str = .true.
264    new_lag = .true.
265    do islope = 1,nslope
266        do i = 1,ngrid
267            current(i,islope)%p => layerings_map(i,islope)%top
268        end do
269    end do
270    ! Conversion to surface ice
271    call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
272    ! We put the sublimating tendency coming from subsurface ice into the overall tendency
273    !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot
274end if
275if (nslope == 1) then ! No slope
276    allocate(character(0) :: num)
277else ! Using slopes
278    allocate(character(8) :: num)
279end if
280allocate(delta_icetable(ngrid),icetable_depth_old(ngrid,nslope),is_co2ice_disappeared(ngrid,nslope),tsoil_ts_old(ngrid,nsoil,nslope,nday))
281is_co2ice_disappeared(:,:) = .false.
282delta_icetable(:) = 0._dp
283n_yr_run = 0
284idt = 0
285do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < nmax_yr_sim)
286    call print_msg('**** Iteration of the PEM leg (Planetary years): '//real2str(n_yr_run + dt))
287    ! Evolve global surface pressure
288    call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)
289
290    ! Adapt tracers according to global surface pressure
291    call adapt_tracers2pressure(ps_avg_glob_old,ps_avg_glob,ps_ts,q_h2o_ts,q_co2_ts)
292
293    ! Evolve surface ice
294    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
295    if (do_layering) then
296        h2oice_depth_old(:,:) = h2oice_depth(:,:)
297        do islope = 1,nslope
298            do i = 1,ngrid
299                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
300                !call print_layering(layerings_map(i,islope))
301            end do
302        end do
303        ! Conversion to surface ice
304        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
305        ! Balance H2O ice reservoirs
306        allocate(d_h2oice_new(ngrid,nslope))
307        call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
308        call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
309        deallocate(d_h2oice_new)
310    else
311        zlag(:,:) = 0._dp
312        call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
313        call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
314    end if
315
316    ! Flow glaciers according to surface ice
317    if (co2ice_flow) then
318        allocate(is_co2ice_flow(ngrid,nslope))
319        call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
320    end if
321    if (h2oice_flow) then
322        allocate(is_h2oice_flow(ngrid,nslope))
323        call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
324    end if
325    if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map)
326
327    ! Adapt surface temperature if surface ice disappeared
328    call adapt_tsurf2disappearedice(co2_ice,is_co2ice_ini,is_co2ice_disappeared,tsurf_avg)
329    !call adapt_tsurf2disappearedice(h2o_ice,is_h2oice_ini,tsurf_avg)
330
331    if (do_soil) then
332        ! Shift soil temperature to surface
333        call shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf_avg,tsoil_avg)
334
335        ! Evolve soil temperature
336        call evolve_soil_temp(tsoil_avg,tsurf_avg,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)
337
338        ! Evolve ice table
339        call evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable,q_h2o_ts, &
340                              ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)
341
342        ! Update soil thermal properties according to ice table and soil temperature
343        call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
344
345        ! Evolve adsorbed species
346        if (do_sorption) then
347            call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
348            call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
349        else
350            totmass_adsh2o = 0._dp
351            totmass_adsco2 = 0._dp
352        end if
353    end if ! do_soil
354    deallocate(zshift_surf,zlag)
355
356    ! Output the results
357    call print_msg('> Standard outputs')
358    call write_diagevol('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
359    do islope = 1,nslope
360        if (nslope /= 1) then
361            num = '  '
362            write(num,'(i2.2)') islope
363            num = '_slope'//num
364        end if
365        call write_diagevol('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
366        call write_diagevol('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
367        call write_diagevol('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
368        call write_diagevol('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
369        if (co2ice_flow) then
370            call write_diagevol('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
371            deallocate(is_co2ice_flow)
372        end if
373        if (h2oice_flow) then
374            call write_diagevol('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
375            deallocate(is_h2oice_flow)
376        end if
377        call write_diagevol('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
378        if (do_soil) then
379            if (icetable_equilibrium) then
380                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
381                call write_diagevol('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
382            else if (icetable_dynamic) then
383                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
384                call write_diagevol('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
385            end if
386            call write_diagevol('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
387            call write_diagevol('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
388            if (do_sorption) then
389                call write_diagevol('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
390                call write_diagevol('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
391            end if
392        end if
393    end do
394
395    ! Checking mass balance for CO2
396    if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then
397        totmass_co2ice = 0._dp
398        totmass_atmco2 = 0._dp
399        do i = 1,ngrid
400            totmass_atmco2 = totmass_atmco2 + cell_area(i)*ps_avg(i)/g
401            do islope = 1,nslope
402                totmass_co2ice = totmass_co2ice + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
403            end do
404        end do
405        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
406        call print_msg(" > Relative total CO2 mass balance = "//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %')
407        if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then
408            call print_msg('  /!\ Warning: mass balance is not conserved!')
409            totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
410            call print_msg('       Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %')
411            totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0
412            call print_msg('       CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %')
413            totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
414            call print_msg('       Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %')
415        end if
416    end if
417
418    ! Evolve the tendencies
419    call evolve_tend_co2ice(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
420!~     if (do_layering) then
421!~         do i = 1,ngrid
422!~             do islope = 1,nslope
423!~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2oice(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
424!~             end do
425!~         end do
426!    else
427!        do i = 1,ngrid
428!            do islope = 1,nslope
429!                call evolve_tend_h2oice(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
430!            end do
431!        end do
432!~     end if
433
434    ! Increment time
435    n_yr_run = n_yr_run + dt
436    n_yr_sim = n_yr_sim + dt
437    idt = idt + 1
438
439    ! Check the stopping criteria
440    call print_msg("> Checking the stopping criteria")
441    call stopping_crit_pressure(ps_avg_glob_ini,ps_avg_glob,stopcrit)
442    call stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit)
443    call stopping_crit_co2ice(co2ice_sublim_coverage_ini,co2_ice,d_co2ice,stopcrit)
444    call system_clock(c2)
445    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_run) stopcrit%nmax_yr_run_reached = .true.
446    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_runorb) stopcrit%nmax_yr_runorb_reached = .true.
447    if (stopcrit%stop_code() == 0 .and. n_yr_sim >= nmax_yr_sim) stopcrit%nmax_yr_sim_reached = .true.
448    if (stopcrit%stop_code() == 0 .and. timewall .and. real(c2 - c1,dp)/real(cr,dp) >= timelimit - antetime) stopcrit%time_limit_reached = .true.
449    if (stopcrit%is_any_set()) then
450        call print_msg(stopcrit%stop_message())
451        exit
452    else
453        call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.')
454        call print_msg('**** The chained simulation has achieved '//real2str(n_yr_sim)//' Planetary years.')
455        call print_msg('**** This run is continuing!')
456        call print_msg('****')
457    end if
458end do ! End of the evolution loop
459deallocate(is_co2ice_disappeared)
460
461! Finalization
462! ~~~~~~~~~~~~
463call print_msg('')
464call print_msg('********* Finalization *********')
465! Update orbital parameters
466if (evol_orbit) call update_orbit(n_yr_sim,n_yr_run)
467
468! Build ice for the PCM
469allocate(h2o_ice4PCM(ngrid,nslope),co2_ice4PCM(ngrid,nslope),is_h2o_perice(ngrid))
470call build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM)
471
472! Build surface temperature for the PCM
473allocate(tsurf4PCM(ngrid,nslope))
474call build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)
475
476! Build soil for the PCM
477if (do_soil) then
478    allocate(tsoil4PCM(ngrid,nsoil_PCM,nslope),inertiesoil4PCM(ngrid,nsoil_PCM,nslope),flux_geo4PCM(ngrid,nslope))
479    call build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
480end if
481
482! Build atmosphere for the PCM
483allocate(ps4PCM(ngrid),teta4PCM(ngrid,nlayer),air_mass4PCM(ngrid,nlayer))
484call build4PCM_atmosphere(ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini,ps4PCM,pa4PCM,preff4PCM,teta4PCM,air_mass4PCM)
485
486! Build tracers for the PCM
487allocate(q4PCM(ngrid,nlayer,nq))
488call build4PCM_tracers(ps4PCM,q4PCM)
489
490! Build surface radiative properties state for the PCM
491allocate(albedo4PCM(ngrid,nslope),emissivity4PCM(ngrid,nslope))
492call build4PCM_surf_rad_prop(h2o_ice,co2_ice,albedo4PCM,emissivity4PCM)
493
494! Write the "startpem.nc"
495call write_restartpem(h2o_ice,co2_ice,tsoil_avg,TI,icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)
496
497! Write the "startfi.nc"
498call write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM)
499
500! Write the "start.nc"
501call write_restart(ps4PCM,pa4PCM,preff4PCM,q4PCM,teta4PCM,air_mass4PCM)
502
503! Update the duration information of the simulation
504call update_info(n_yr_run,stopcrit%stop_code(),n_yr_sim,nmax_yr_sim)
505
506! Deallocation
507deallocate(emissivity4PCM,albedo4PCM)
508deallocate(q4PCM)
509deallocate(ps4PCM,teta4PCM,air_mass4PCM)
510if (do_soil) deallocate(tsoil4PCM,inertiesoil4PCM,flux_geo4PCM)
511deallocate(tsurf4PCM)
512deallocate(co2_ice4PCM,h2o_ice4PCM,is_h2o_perice)
513deallocate(num)
514if (do_layering) then
515    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
516    do islope = 1,nslope
517        do i = 1,ngrid
518            call del_layering(layerings_map(i,islope))
519        end do
520    end do
521end if
522deallocate(layerings_map)
523deallocate(h2o_ads_reg,co2_ads_reg)
524deallocate(h2o_ice,co2_ice,is_h2oice_ini,is_co2ice_ini)
525deallocate(ps_avg,ps_dev,ps_ts)
526deallocate(tsurf_avg,tsurf_dev,h2o_surfdensity_avg)
527deallocate(tsoil_avg,tsoil_dev,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)
528deallocate(q_h2o_ts,q_co2_ts,q_co2_ts_ini)
529deallocate(d_h2oice,d_co2ice,d_co2ice_ini)
530deallocate(delta_h2o_ads,delta_co2_ads,delta_icetable,icetable_depth_old)
531call end_allocation()
532
533! Footer
534call system_clock(c2)
535call print_end(n_yr_run,real((c2 - c1),dp)/real(cr,dp),r_plnt2earth_yr,pem_ini_date,n_yr_sim)
536
537END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.