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

Last change on this file since 4068 was 4068, checked in by jbclement, 2 weeks ago

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