source: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 @ 3154

Last change on this file since 3154 was 3149, checked in by jbclement, 2 years ago

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC
File size: 25.1 KB
Line 
1MODULE pemetat0_mod
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, &
10                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,      &
11                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,              &
12                    m_h2o_regolith_phys,deltam_h2o_regolith_phys)
13
14use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
15use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, index_breccia, index_bedrock
16use comsoil_h,                  only: volcapa, inertiedat
17use adsorption_mod,             only: regolith_adsorption, adsorption_pem
18use ice_table_mod,              only: computeice_table_equilibrium
19use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
20use soil_thermalproperties_mod, only: update_soil_thermalproperties
21use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
22use abort_pem_mod,              only: abort_pem
23use soil_pem_compute_mod,       only: soil_pem_compute
24use soil_pem_ini_mod,           only: soil_pem_ini
25
26#ifndef CPP_STD
27    use comcstfi_h,   only: r, mugaz
28    use surfdat_h,    only: watercaptag, perennial_co2ice
29#else
30    use comcstfi_mod, only: r, mugaz
31#endif
32
33implicit none
34
35include "callkeys.h"
36
37character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
38integer,                        intent(in) :: ngrid               ! # of physical grid points
39integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
40integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
41integer,                        intent(in) :: nslope              ! # of sub-grid slopes
42integer,                        intent(in) :: timelen             ! # time samples
43real,                           intent(in) :: timestep            ! time step [s]
44real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
45real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
46real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
47real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
48real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
49real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
50real, dimension(ngrid,nslope),  intent(in) :: tend_h2o_ice        ! tendencies on h2o ice
51real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! tendencies on co2 ice
52real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
53! outputs
54real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
55real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
56real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
57real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
58real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
59real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
60real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
61real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
62real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
63real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
65real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
66! local
67real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
68real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
69logical                                 :: found                        ! check if variables are found in the start
70logical                                 :: found2                       ! check if variables are found in the start
71integer                                 :: iloop, ig, islope, it, isoil ! index for loops
72real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
73real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
74character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
75real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
76real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr 1 [K]
77real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
78real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
79real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
80logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
81
82#ifdef CPP_STD
83    logical, dimension(ngrid) :: watercaptag
84    watercaptag(:) = .false.
85#endif
86
87!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88!!!
89!!! Purpose: read start_pem. Need a specific iostart_PEM
90!!!
91!!! Order: 0. Previous year of the PEM run
92!!!           Ice initialization
93!!!        1. Thermal Inertia
94!!!        2. Soil Temperature
95!!!        3. Ice table
96!!!        4. Mass of CO2 & H2O adsorbed
97!!!
98!!! /!\ This order must be respected !
99!!! Author: LL
100!!!
101!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
102
103!0. Check if the start_PEM exist.
104
105inquire(file = filename,exist = startpem_file)
106
107write(*,*)'Is start PEM?',startpem_file
108
109!1. Run
110if (startpem_file) then
111    ! open pem initial state file:
112    call open_startphy(filename)
113
114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115#ifndef CPP_STD
116    ! h2o ice
117    h2o_ice = 0.
118    call get_field("h2o_ice",h2o_ice,found)
119    if (.not. found) then
120        write(*,*)'Pemetat0: failed loading <h2o_ice>'
121        write(*,*)'will reconstruct the values from watercaptag'
122        write(*,*)'with default value ''ini_h2o_bigreservoir'''
123        do ig = 1,ngrid
124            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
125        enddo
126    endif
127
128    ! co2 ice
129    co2_ice = perennial_co2ice
130#endif
131
132    if (soil_pem) then
133
134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135!1. Thermal Inertia
136! a. General case
137        do islope = 1,nslope
138            write(num,'(i2.2)') islope
139            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
140            if (.not. found) then
141                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
142                write(*,*)'will reconstruct the values of TI_PEM'
143
144                do ig = 1,ngrid
145                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
146                        !!! transition
147                        delta = depth_breccia
148                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
149                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
150                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
151                        do iloop=index_breccia+2,index_bedrock
152                            TI_PEM(ig,iloop,islope) = TI_breccia
153                        enddo
154                    else ! we keep the high ti values
155                        do iloop = index_breccia + 1,index_bedrock
156                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
157                        enddo
158                    endif ! TI PEM and breccia comparison
159                    !! transition
160                    delta = depth_bedrock
161                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
162                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
163                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
164                    do iloop = index_bedrock + 2,nsoil_PEM
165                        TI_PEM(ig,iloop,islope) = TI_bedrock
166                    enddo
167                enddo
168            else ! found
169                do iloop = nsoil_PCM + 1,nsoil_PEM
170                    TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
171                enddo
172            endif ! not found
173        enddo ! islope
174
175        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
176
177! b. Special case for inertiedat, inertiedat_PEM
178        call get_field("inertiedat_PEM",inertiedat_PEM,found)
179        if (.not.found) then
180            do iloop = 1,nsoil_PCM
181                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
182            enddo
183        !!! zone de transition
184            delta = depth_breccia
185            do ig = 1,ngrid
186                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
187                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
188                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
189                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
190
191                    do iloop = index_breccia+2,index_bedrock
192                        inertiedat_PEM(ig,iloop) = TI_breccia
193                    enddo
194                else
195                    do iloop=index_breccia+1,index_bedrock
196                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
197                    enddo
198                endif ! comparison ti breccia
199            enddo ! ig
200
201            !!! zone de transition
202            delta = depth_bedrock
203            do ig = 1,ngrid
204                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
205                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
206                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
207            enddo
208
209            do iloop = index_bedrock + 2, nsoil_PEM
210                do ig = 1,ngrid
211                    inertiedat_PEM(ig,iloop) = TI_bedrock
212                enddo
213            enddo
214        endif ! not found
215
216!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
217!2. Soil Temperature
218        do islope=1,nslope
219            write(num,fmt='(i2.2)') islope
220            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
221            if (.not.found) then
222                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
223                write(*,*)'will reconstruct the values of Tsoil'
224!                do ig = 1,ngrid
225!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
226!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
227!                    do iloop=index_breccia+2,index_bedrock
228!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
229!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
230!                    enddo
231!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
232!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
233!
234!                    do iloop=index_bedrock+2,nsoil_PEM
235!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
236!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
237!                    enddo
238!                enddo
239                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
240                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
241            else
242! predictor corrector: restart from year 1 of the PCM and build the evolution of
243! tsoil at depth
244                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
245                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
246                call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
247                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
248                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
249
250                do iloop = nsoil_PCM+1,nsoil_PEM
251                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
252                enddo
253            endif !found
254
255            do it = 1,timelen
256                do isoil = nsoil_PCM+1,nsoil_PEM
257                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
258                enddo
259            enddo
260            do isoil = nsoil_PCM+1,nsoil_PEM
261                do ig = 1,ngrid
262                    watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
263                enddo
264            enddo
265        enddo ! islope
266        write(*,*) 'PEMETAT0: SOIL TEMP done'
267
268!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
269!3. Ice Table
270        call get_field("ice_table",ice_table,found)
271        if (.not. found) then
272            write(*,*)'PEM settings: failed loading <ice_table>'
273            write(*,*)'will reconstruct the values of the ice table given the current state'
274            ice_table = -1  ! by default, no ice table
275            ice_table_thickness = -1  ! by default, no ice table
276            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
277            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
278            do islope = 1,nslope
279                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
280            enddo
281        endif
282
283        write(*,*) 'PEMETAT0: ICE TABLE done'
284
285!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
286!4. CO2 & H2O Adsorption
287        if (adsorption_pem) then
288            do islope = 1,nslope
289                write(num,fmt='(i2.2)') islope
290                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
291                if (.not. found) then
292                    m_co2_regolith_phys = 0.
293                    exit
294                endif
295            enddo
296            do islope=1,nslope
297                write(num,fmt='(i2.2)') islope
298                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
299                if (.not. found2) then
300                    m_h2o_regolith_phys = 0.
301                    exit
302                endif
303            enddo
304
305            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
306                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
307
308            if (.not. found) deltam_co2_regolith_phys = 0.
309            if (.not.found2) deltam_h2o_regolith_phys = 0.
310
311            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
312        endif ! adsorption_pem
313    endif ! soil_pem
314
315    call close_startphy
316
317else !No startfi, let's build all by hand
318
319    ! h2o ice
320    h2o_ice = 0.
321    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_h2o_bigreservoir''.'
322        do ig = 1,ngrid
323            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
324        enddo
325
326    ! co2 ice
327    co2_ice = perennial_co2ice
328
329    if (soil_pem) then
330
331!a) Thermal inertia
332        do islope = 1,nslope
333            do ig = 1,ngrid
334                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
335                    !!! transition
336                    delta = depth_breccia
337                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
338                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
339                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
340                    do iloop=index_breccia + 2,index_bedrock
341                        TI_PEM(ig,iloop,islope) = TI_breccia
342                    enddo
343                else ! we keep the high ti values
344                    do iloop=index_breccia + 1,index_bedrock
345                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
346                    enddo
347                endif
348                !! transition
349                delta = depth_bedrock
350                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
351                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
352                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
353                do iloop = index_bedrock + 2,nsoil_PEM
354                    TI_PEM(ig,iloop,islope) = TI_bedrock
355                enddo
356            enddo
357        enddo
358
359        do iloop = 1,nsoil_PCM
360            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
361        enddo
362        !!! zone de transition
363        delta = depth_breccia
364        do ig = 1,ngrid
365            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
366                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
367                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
368                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
369                do iloop = index_breccia + 2,index_bedrock
370                    inertiedat_PEM(ig,iloop) = TI_breccia
371                enddo
372            else
373                do iloop = index_breccia + 1,index_bedrock
374                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
375                enddo
376            endif
377        enddo
378
379        !!! zone de transition
380        delta = depth_bedrock
381        do ig = 1,ngrid
382            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
383                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
384                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
385        enddo
386
387        do iloop = index_bedrock + 2,nsoil_PEM
388            do ig = 1,ngrid
389                inertiedat_PEM(ig,iloop) = TI_bedrock
390            enddo
391        enddo
392
393        write(*,*) 'PEMETAT0: TI done'
394
395!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
396!b) Soil temperature
397        do islope = 1,nslope
398!            do ig = 1,ngrid
399!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
400!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
401!
402!                do iloop=index_breccia+2,index_bedrock
403!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
404!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
405!                enddo
406!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
407!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
408!
409!                do iloop=index_bedrock+2,nsoil_PEM
410!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
411!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
412!                enddo
413!            enddo
414            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
415            call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
416
417! First raw initialization
418            do it = 1,timelen
419                do isoil = nsoil_PCM + 1,nsoil_PEM
420                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
421                enddo
422            enddo
423
424            do it = 1,timelen
425                do isoil = nsoil_PCM + 1,nsoil_PEM
426                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
427                enddo
428            enddo
429
430            do isoil = nsoil_PCM + 1,nsoil_PEM
431                do ig = 1,ngrid
432                    watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
433                enddo
434            enddo
435        enddo !islope
436        write(*,*) 'PEMETAT0: TSOIL done'
437
438!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
439!c) Ice table
440        ice_table = -1.  ! by default, no ice table
441        ice_table_thickness = -1.
442        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
443        call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
444        do islope = 1,nslope
445            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
446        enddo
447        write(*,*) 'PEMETAT0: Ice table done'
448
449!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450!d) Regolith adsorbed
451        if (adsorption_pem) then
452            m_co2_regolith_phys = 0.
453            m_h2o_regolith_phys = 0.
454            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
455                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
456            deltam_co2_regolith_phys = 0.
457            deltam_h2o_regolith_phys = 0.
458        endif
459
460        write(*,*) 'PEMETAT0: CO2 adsorption done'
461    endif !soil_pem
462endif ! of if (startphy_file)
463
464if (soil_pem) then ! Sanity check
465    do ig = 1,ngrid
466        do islope = 1,nslope
467            do iloop = 1,nsoil_PEM
468                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
469            enddo
470        enddo
471    enddo
472endif !soil_pem
473
474END SUBROUTINE pemetat0
475
476END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.