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

Last change on this file since 4174 was 4174, checked in by jbclement, 3 days ago

PEM:
Refactor of H2O flux balancing:

  • centralizes flux balancing logic;
  • removes intermediate variables, especially related to of H2O mass bookkeeping;
  • adds new stopping criterion when H2O flux balance fails8 (sanity check);
  • introduces configurable weighting strategies for H2O flux balancing.

JBC
File size: 21.3 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 with r2778 & r2779
11!     L. Lange, 22/07/2022
12!     JB Clement, 2023-2025
13!
14! NOTES
15!     Ownership criterion for declarations:
16!       - Declare in module "planet" variables that are persistent climate
17!         state or persistent references used across multiple time steps.
18!       - Declare in program "pem.F90" transient workflow/control varaibles,
19!         one-shot initialization buffers and per-iteration working buffers.
20!-----------------------------------------------------------------------
21
22! DEPENDENCIES
23! ------------
24! Common modules
25use job_mod,        only: timelimit, antetime, timewall
26use parse_args_mod, only: parse_args
27! PEM modules
28use allocation,         only: ini_allocation, end_allocation
29use atmosphere,         only: ps_PCM, evolve_pressure, CO2cond_ps_PCM
30use backup,             only: save_clim_state, backup_rate
31use clim_state_init,    only: read_start, read_startfi, read_startevo
32use config,             only: read_rundef, read_display_config
33use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
34use evolution,          only: n_yr_run, n_yr_sim, ntot_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
35use geometry,           only: ngrid, nslope, nsoil_PCM, nsoil, cell_area, total_surface
36use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
37use ice_table,          only: icetable_equilibrium, icetable_dynamic, evolve_ice_table
38use layered_deposits,   only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
39use maths,              only: pi
40use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
41use orbit,              only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
42use output,             only: write_diagevo, dim_ngrid, dim_nsoil
43use planet
44use physics,            only: g
45use slopes,             only: subslope_dist, def_slope_mean
46use soil,               only: do_soil, set_soil, TI
47use soil_temp,          only: tsoil_PCM, shift_tsoil2surf, evolve_soil_temp
48use soil_therm_inertia, only: update_soil_TI
49use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
50use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2oice, stopping_crit_co2ice
51use surface,            only: emissivity_PCM
52use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2o_fluxes
53use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice
54use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice
55use tracers,            only: adapt_tracers2pressure
56use utility,            only: real2str
57use workflow_status,    only: i_pem_run, read_workflow_status, update_workflow_status
58use xios_data,          only: load_xios_data
59
60! DECLARATION
61! -----------
62implicit none
63
64! LOCAL VARIABLES
65! ---------------
66! Utility-related:
67integer(li)  :: cr     ! Number of clock ticks per second (count rate)
68integer(li)  :: c1, c2 ! Counts of processor clock
69character(8) :: num    ! Slope suffix to ouput variables
70integer(di)  :: i, islope
71! Ice-related:
72logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow   ! Flag for location of CO2 glacier flow
73logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow   ! Flag for location of H2O glacier flow
74real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg/m2]
75real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg/m2]
76real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost  ! Minimum of H2O frost over the last PCM year [kg/m2]
77real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost  ! Minimum of CO2 frost over the last PCM year [kg/m2]
78! Surface-related:
79real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K]
80real(dp), dimension(:,:), allocatable :: zshift_surf   ! Elevation shift for the surface [m]
81real(dp), dimension(:,:), allocatable :: zlag          ! Newly built lag thickness [m]
82! Layering-related:
83type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
84real(dp),       dimension(:,:), allocatable :: h2oice_depth     ! Depth of subsurface ice layer
85real(dp),       dimension(:,:), allocatable :: h2oice_depth_old ! Old depth of subsurface ice layer
86logical(k4),    dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm
87! Evolution-related:
88real(dp)        :: nmax_yr_runorb ! Maximum number of years for the run due to orbital parameters
89type(stopFlags) :: stopcrit       ! Stopping criteria
90! Balance-related
91real(qp) :: totmass_co2ice, totmass_atmco2, totmass_adsco2             ! Current total CO2 masses (surface ice|atmospheric|adsorbed)
92real(qp) :: totmass_co2ice_ini, totmass_atmco2_ini, totmass_adsco2_ini ! Initial total CO2 masses (surface ice|atmospheric|adsorbed)
93real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
94real(qp) :: totmass_adsh2o                                             ! Current total adsorbed H2O mass [kg]
95
96! CODE
97! ----
98! Pre-processing step
99!~~~~~~~~~~~~~~~~~~~~
100! Elapsed time with system clock
101call system_clock(count_rate = cr)
102call system_clock(c1)
103
104! Parse command-line options
105call parse_args()
106
107! Initialization
108! ~~~~~~~~~~~~~~
109! Header
110call read_display_config()
111call print_ini()
112
113! Allocate module arrays
114call ini_allocation()
115
116! Read the duration information of the workflow
117call read_workflow_status()
118
119! Read the PEM parameters
120call read_rundef()
121
122! Read the orbital parameters
123call read_orbitpm(n_yr_sim)
124
125! Read the "start.nc"
126call read_start()
127
128! Read the "startfi.nc"
129call read_startfi()
130
131! Read the PCM data given by XIOS
132call allocate_xios_state()
133allocate(tsurf_avg_yr1(ngrid,nslope))
134allocate(minPCM_h2operice(ngrid,nslope,2))
135allocate(minPCM_co2perice(ngrid,nslope,2))
136allocate(minPCM_h2ofrost(ngrid,nslope,2))
137allocate(minPCM_co2frost(ngrid,nslope,2))
138
139call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
140                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost,flux_ssice_avg)
141
142! Initiate soil settings and TI
143if (do_soil) call set_soil(TI)
144
145! Compute the deviation from the average
146call allocate_deviation_state()
147ps_dev(:) = ps_PCM(:) - ps_avg(:)
148tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:)
149tsoil_dev(:,:,:) = tsoil_PCM(:,:,:) - tsoil_avg(:,1:nsoil_PCM,:)
150
151! Compute global surface pressure
152ps_avg_glob = sum(cell_area*ps_avg)/total_surface
153
154! Compute the accepted maximum number of years due to orbital parameters (if needed)
155call compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)
156
157! Read the "startevo.nc"
158call allocate_startevo_state()
159call read_startevo(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,h2o_ice,co2_ice, &
160                   tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map,       &
161                   h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
162deallocate(tsurf_avg_yr1)
163
164! Compute ice tendencies from yearly minima
165call allocate_tendencies()
166call print_msg('> Computing surface ice tendencies',LVL_NFO)
167call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice)
168call print_msg('H2O ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
169call compute_tendice(minPCM_co2perice,minPCM_co2frost,co2_ice,d_co2ice)
170call print_msg('CO2 ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
171deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
172
173! Save initial set-up useful for the next computations
174call print_msg('> Saving some initial climate state variables',LVL_NFO)
175call allocate_initial_snapshots()
176ps_avg_glob_ini = ps_avg_glob
177d_co2ice_ini(:,:) = d_co2ice(:,:)
178q_co2_ts_ini(:,:) = q_co2_ts(:,:)
179h2oice_sublim_coverage_ini = 0._dp
180co2ice_sublim_coverage_ini = 0._dp
181totmass_co2ice_ini = 0._qp
182totmass_atmco2_ini = 0._qp
183totmass_adsco2_ini = 0._qp
184totmass_adsh2o = 0._qp
185do i = 1,ngrid
186    totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
187    do islope = 1,nslope
188        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)
189        if (co2_ice(i,islope) > 0._dp) then
190            is_co2ice_ini(i,islope) = .true.
191            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)
192        end if
193        if (h2o_ice(i,islope) > 0._dp) then
194            is_h2oice_ini(i,islope) = .true.
195            if (d_h2oice(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)
196        end if
197    end do
198end do
199call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini),LVL_NFO)
200call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini),LVL_NFO)
201call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini),LVL_NFO)
202if (do_sorption) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2_ini,totmass_adsh2o)
203
204! Main evolution loop
205! ~~~~~~~~~~~~~~~~~~~
206call print_msg('',LVL_NFO)
207call print_msg('********* Evolution *********',LVL_NFO)
208if (do_layering) then
209    allocate(h2oice_depth(ngrid,nslope),h2oice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
210    new_str = .true.
211    new_lag = .true.
212    do islope = 1,nslope
213        do i = 1,ngrid
214            current(i,islope)%p => layerings_map(i,islope)%top
215        end do
216    end do
217    ! Conversion to surface ice
218    call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
219end if
220call allocate_loop_state()
221n_yr_run = 0
222idt = 0
223do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim)
224    call print_msg('**** Iteration of the PEM run [plnt y]: '//real2str(n_yr_run + dt),LVL_NFO)
225    ! Evolve global surface pressure
226    call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)
227
228    ! Adapt tracers according to global surface pressure
229    call adapt_tracers2pressure(ps_avg_glob_old,ps_avg_glob,ps_ts,q_h2o_ts,q_co2_ts)
230
231    ! Evolve surface ice
232    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
233    if (do_layering) then
234        h2oice_depth_old(:,:) = h2oice_depth(:,:)
235        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
236        where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:)
237        do islope = 1,nslope
238            do i = 1,ngrid
239                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
240                call print_layering(layerings_map(i,islope))
241            end do
242        end do
243        ! Conversion to surface ice
244        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
245        ! Balance H2O ice reservoirs
246        call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit)
247        if (stopcrit%stop_code() > 0) then
248            deallocate(zshift_surf,zlag)
249            call print_msg(stopcrit%stop_message(),LVL_NFO)
250            exit
251        end if
252    else
253        zlag(:,:) = 0._dp
254        zshift_surf(:,:) = 0._dp
255        call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
256        call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
257    end if
258
259    ! Flow glaciers according to surface ice
260    if (co2ice_flow) then
261        allocate(is_co2ice_flow(ngrid,nslope))
262        call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
263    end if
264    if (h2oice_flow) then
265        allocate(is_h2oice_flow(ngrid,nslope))
266        call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
267    end if
268    if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map)
269
270    ! Adapt surface temperature if surface ice disappeared
271    call adapt_tsurf2disappearedice(co2_ice,is_co2ice_ini,is_co2ice_disappeared,tsurf_avg)
272    !call adapt_tsurf2disappearedice(h2o_ice,is_h2oice_ini,tsurf_avg)
273
274    if (do_soil) then
275        ! Shift soil temperature to surface
276        call shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf_avg,tsoil_avg)
277
278        ! Evolve soil temperature
279        call evolve_soil_temp(tsoil_avg,tsurf_avg,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)
280
281        ! Evolve ice table
282        call evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable,q_h2o_ts, &
283                              ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)
284
285        ! Update soil thermal properties according to ice table and soil temperature
286        call update_soil_TI(h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
287
288        ! Evolve adsorbed species
289        if (do_sorption) then
290            call evolve_regolith_adsorption(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)
291            call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
292        else
293            totmass_adsh2o = 0._dp
294            totmass_adsco2 = 0._dp
295        end if
296    end if ! do_soil
297    deallocate(zshift_surf,zlag)
298
299    ! Output the results
300    call print_msg('> Standard outputs',LVL_NFO)
301    call write_diagevo('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
302    do islope = 1,nslope
303        if (nslope == 1) then
304            num = ''
305        else
306            write(num,'("_slope",i2.2)') islope
307        end if
308        call write_diagevo('h2oice'//trim(num),'H2O ice','kg/m2',h2o_ice(:,islope),(/dim_ngrid/))
309        call write_diagevo('co2ice'//trim(num),'CO2 ice','kg/m2',co2_ice(:,islope),(/dim_ngrid/))
310        call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg/m2/y',d_h2oice(:,islope),(/dim_ngrid/))
311        call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg/m2/y',d_co2ice(:,islope),(/dim_ngrid/))
312        if (co2ice_flow) then
313            call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
314        end if
315        if (h2oice_flow) then
316            call write_diagevo('Flow_h2oice'//trim(num),'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
317        end if
318        call write_diagevo('tsurf'//trim(num),'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
319        if (do_soil) then
320            if (icetable_equilibrium) then
321                call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
322                call write_diagevo('icetable_thick'//trim(num),'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
323            else if (icetable_dynamic) then
324                call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
325                call write_diagevo('ice_porefilling'//trim(num),'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
326            end if
327            call write_diagevo('tsoil_avg'//trim(num),'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
328            call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
329            if (do_sorption) then
330                call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg/m2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
331                call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg/m2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
332            end if
333        end if
334    end do
335    if (co2ice_flow .and. allocated(is_co2ice_flow)) deallocate(is_co2ice_flow)
336    if (h2oice_flow .and. allocated(is_h2oice_flow)) deallocate(is_h2oice_flow)
337
338    ! Checking mass balance for CO2
339    if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then
340        totmass_co2ice = 0._dp
341        totmass_atmco2 = 0._dp
342        do i = 1,ngrid
343            totmass_atmco2 = totmass_atmco2 + cell_area(i)*ps_avg(i)/g
344            do islope = 1,nslope
345                totmass_co2ice = totmass_co2ice + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
346            end do
347        end do
348        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
349        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)//' %',LVL_NFO)
350        if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then
351            call print_msg('Mass balance is not conserved!',LVL_WRN)
352            totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
353            call print_msg('    Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %',LVL_WRN)
354            totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0
355            call print_msg('    CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %',LVL_WRN)
356            totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
357            call print_msg('    Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %',LVL_WRN)
358        end if
359    end if
360
361    ! Evolve the tendencies
362    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)
363    if (do_layering) then
364        call evolve_flux_ssice(h2oice_depth_old,h2oice_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
365    else if (icetable_equilibrium .or. icetable_dynamic) then
366        call evolve_flux_ssice(icetable_depth_old,icetable_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
367    end if
368
369    ! Increment time
370    n_yr_run = n_yr_run + dt
371    n_yr_sim = n_yr_sim + dt
372    idt = idt + 1
373
374    ! Save periodic backups of restart files
375    if (backup_rate > 0) then
376        if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
377                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
378    end if
379
380    ! Check the stopping criteria
381    call print_msg("> Checking the stopping criteria",LVL_NFO)
382    call stopping_crit_pressure(ps_avg_glob_ini,ps_avg_glob,stopcrit)
383    call stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit)
384    call stopping_crit_co2ice(co2ice_sublim_coverage_ini,co2_ice,d_co2ice,stopcrit)
385    call system_clock(c2)
386    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_run) stopcrit%nmax_yr_run_reached = .true.
387    if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_runorb) stopcrit%nmax_yr_runorb_reached = .true.
388    if (stopcrit%stop_code() == 0 .and. n_yr_sim >= ntot_yr_sim) stopcrit%nmax_yr_sim_reached = .true.
389    if (stopcrit%stop_code() == 0 .and. timewall .and. real(c2 - c1,dp)/real(cr,dp) >= timelimit - antetime) stopcrit%time_limit_reached = .true.
390    if (stopcrit%is_any_set()) then
391        call print_msg(stopcrit%stop_message(),LVL_NFO)
392        exit
393    else
394        call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.',LVL_NFO)
395        call print_msg('**** The workflow has achieved '//real2str(n_yr_sim)//' Planetary years.',LVL_NFO)
396        call print_msg('**** This run is continuing!',LVL_NFO)
397        call print_msg('****',LVL_NFO)
398    end if
399end do ! End of the evolution loop
400
401! Finalization
402! ~~~~~~~~~~~~
403call print_msg('',LVL_NFO)
404call print_msg('********* Finalization *********',LVL_NFO)
405! Update orbital parameters
406if (evo_orbit) call update_orbit(n_yr_sim,n_yr_run)
407
408call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
409                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map)
410
411! Update the duration information of the workflow
412call update_workflow_status(n_yr_run,stopcrit%stop_code(),n_yr_sim,ntot_yr_sim)
413
414! Deallocation
415if (do_layering) then
416    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
417    do islope = 1,nslope
418        do i = 1,ngrid
419            call del_layering(layerings_map(i,islope))
420        end do
421    end do
422end if
423call end_allocation()
424
425! Footer
426call system_clock(c2)
427call print_end(i_pem_run,n_yr_run,n_yr_sim,real((c2 - c1),dp)/real(cr,dp),pem_ini_date,r_plnt2earth_yr)
428
429END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.