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

Last change on this file since 4074 was 4074, checked in by jbclement, 9 days ago

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout 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, ntot_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 layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering
35use maths,              only: pi
36use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
37use orbit,              only: evol_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
38use output,             only: write_diagevol, dim_ngrid, dim_nsoil
39use physics,            only: g
40use slopes,             only: subslope_dist, def_slope_mean
41use soil,               only: do_soil, set_soil, TI, build4PCM_soil
42use soil_temp,          only: tsoil_PCM, shift_tsoil2surf, evolve_soil_temp
43use soil_therm_inertia, only: update_soil_TI
44use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
45use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice
46use surface,            only: emissivity_PCM, build4PCM_surf_rad_prop
47use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice
48use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf
49use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
50use tracers,            only: adapt_tracers2pressure, build4PCM_tracers, nq
51use utility,            only: real2str
52use workflow_status,    only: read_workflow_status, update_workflow_status
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 workflow
166call read_workflow_status()
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_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)
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 < ntot_yr_sim)
286    call print_msg('**** Iteration of the PEM run [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 >= ntot_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 workflow 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 workflow
504call update_workflow_status(n_yr_run,stopcrit%stop_code(),n_yr_sim,ntot_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.