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

Last change on this file since 3583 was 3571, checked in by jbclement, 11 months ago

PEM:

  • New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
  • Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
  • Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
  • Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
  • Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
  • Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
  • Cosmetic cleanings for readability.

JBC
File size: 25.0 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,icetable_depth,icetable_thickness, &
10                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,                           &
11                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                                      &
12                    d_h2oice,d_co2ice,co2_ice,h2o_ice,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
13
14use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
15use comsoil_h_PEM,              only: soil_pem, layer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
16use comsoil_h,                  only: inertiedat
17use adsorption_mod,             only: regolith_adsorption, adsorption_pem
18use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
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 compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
24use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo
25
26#ifndef CPP_STD
27    use comcstfi_h,   only: r, mugaz, pi
28    use surfdat_h,    only: watercaptag, watercap, perennial_co2ice
29#else
30    use comcstfi_mod, only: r, mugaz, pi
31#endif
32
33implicit none
34
35include "callkeys.h"
36
37character(*),                   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) :: ps_avg_global ! 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_timeseries ! surface pressure [Pa]
50real, dimension(ngrid,nslope),  intent(in) :: d_h2oice      ! tendencies on h2o ice
51real, dimension(ngrid,nslope),  intent(in) :: d_co2ice      ! 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]
58type(layering), dimension(ngrid,nslope), intent(inout) :: stratif             ! stratification (layerings)
59real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
60real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
61real, dimension(ngrid,nslope),           intent(inout) :: icetable_depth      ! Ice table depth [m]
62real, dimension(ngrid,nslope),           intent(inout) :: icetable_thickness  ! Ice table thickness [m]
63real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
64real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
65real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
66real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
67! local
68real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem           ! soil temperature saved in the start [K]
69real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM              ! soil thermal inertia saved in the start [SI]
70logical                                 :: found, found2            ! check if variables are found in the start
71integer                                 :: iloop, ig, islope, isoil ! index for loops
72real                                    :: delta                    ! Depth of the interface regolith-breccia, breccia -bedrock [m]
73character(2)                            :: num                      ! intermediate string to read PEM start sloped variables
74logical                                 :: startpem_file            ! boolean to check if we read the startfile or not
75real, dimension(:,:,:,:), allocatable   :: stratif_array            ! Array for stratification (layerings)
76
77#ifdef CPP_STD
78    logical, dimension(ngrid) :: watercaptag
79    watercaptag = .false.
80#endif
81
82!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
83!!!
84!!! Purpose: read start_pem. Need a specific iostart_PEM
85!!!
86!!! Order: 0. Previous year of the PEM run
87!!!           Ice initialization
88!!!        1. Thermal Inertia
89!!!        2. Soil Temperature
90!!!        3. Ice table
91!!!        4. Mass of CO2 & H2O adsorbed
92!!!
93!!! /!\ This order must be respected !
94!!! Author: LL
95!!!
96!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97
98!0.1 Check if the start_PEM exist.
99
100inquire(file = filename,exist = startpem_file)
101
102write(*,*)'Is there a "startpem" file?',startpem_file
103
104!0.2 Set to default values
105icetable_depth = -1. ! by default, no ice table
106icetable_thickness = -1.
107ice_porefilling = 0.
108!1. Run
109if (startpem_file) then
110    write(*,*)'There is a "startpem.nc"!'
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_huge_h2oice'''
123        do ig = 1,ngrid
124            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
125        enddo
126    else
127        ! The variations of infinite reservoirs during the PCM years are taken into account
128        h2o_ice = h2o_ice + watercap
129    endif
130
131    ! co2 ice
132    co2_ice = perennial_co2ice
133#endif
134
135!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
136    ! Stratification (layerings)
137    nb_str_max = 1
138    if (layering_algo) then
139        found = inquire_dimension("nb_str_max")
140        if (.not. found) then
141            write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
142            write(*,*) '''nb_str_max'' is set to 3!'
143        else
144            nb_str_max = int(inquire_dimension_length('nb_str_max'))
145        endif
146        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
147        stratif_array = 0.
148        do islope = 1,nslope
149            write(num,'(i2.2)') islope
150            call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
151            call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
152            found = found .or. found2
153            call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
154            found = found .or. found2
155            call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
156            found = found .or. found2
157            call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
158            found = found .or. found2
159            call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
160            found = found .or. found2
161            if (.not. found) then
162                write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
163            endif ! not found
164        enddo ! islope
165        call array2stratif(stratif_array,ngrid,nslope,stratif)
166        deallocate(stratif_array)
167    endif
168
169    if (soil_pem) then
170
171!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
172!1. Thermal Inertia
173! a. General case
174        do islope = 1,nslope
175            write(num,'(i2.2)') islope
176            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
177            if (.not. found) then
178                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
179                write(*,*)'will reconstruct the values of TI_PEM'
180
181                do ig = 1,ngrid
182                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
183                        !!! transition
184                        delta = depth_breccia
185                        TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
186                                                                  (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
187                                                                  ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
188                        do iloop=index_breccia + 2,index_bedrock
189                            TI_PEM(ig,iloop,islope) = TI_breccia
190                        enddo
191                    else ! we keep the high TI values
192                        do iloop = index_breccia + 1,index_bedrock
193                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
194                        enddo
195                    endif ! TI PEM and breccia comparison
196                    !! transition
197                    delta = depth_bedrock
198                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
199                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
200                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
201                    do iloop = index_bedrock + 2,nsoil_PEM
202                        TI_PEM(ig,iloop,islope) = TI_bedrock
203                    enddo
204                enddo
205            else ! found
206                do iloop = nsoil_PCM + 1,nsoil_PEM
207                    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.
208                enddo
209            endif ! not found
210        enddo ! islope
211
212        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
213
214! b. Special case for inertiedat, inertiedat_PEM
215        call get_field("inertiedat_PEM",inertiedat_PEM,found)
216        if (.not. found) then
217            do iloop = 1,nsoil_PCM
218                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
219            enddo
220        !!! zone de transition
221            delta = depth_breccia
222            do ig = 1,ngrid
223                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
224                    inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
225                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
226                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
227
228                    do iloop = index_breccia + 2,index_bedrock
229                        inertiedat_PEM(ig,iloop) = TI_breccia
230                    enddo
231                else
232                    do iloop=index_breccia + 1,index_bedrock
233                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
234                    enddo
235                endif ! comparison ti breccia
236            enddo ! ig
237
238            !!! zone de transition
239            delta = depth_bedrock
240            do ig = 1,ngrid
241                inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
242                                                           (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
243                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
244            enddo
245
246            do iloop = index_bedrock + 2, nsoil_PEM
247                do ig = 1,ngrid
248                    inertiedat_PEM(ig,iloop) = TI_bedrock
249                enddo
250            enddo
251        endif ! not found
252
253!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254!2. Soil Temperature
255        do islope = 1,nslope
256            write(num,fmt='(i2.2)') islope
257            call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
258            if (.not. found) then
259                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
260                write(*,*)'will reconstruct the values of Tsoil'
261                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
262                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
263            else
264! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
265!                      for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
266                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
267                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
268                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
269
270                do iloop = nsoil_PCM + 1,nsoil_PEM
271                    tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
272                enddo
273            endif !found
274
275            do isoil = nsoil_PCM + 1,nsoil_PEM
276                watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
277            enddo
278        enddo ! islope
279        write(*,*) 'PEMETAT0: SOIL TEMP done'
280
281!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282!3. Ice Table
283        if (icetable_equilibrium) then
284            call get_field("icetable_depth",icetable_depth,found)
285            if (.not. found) then
286                write(*,*)'PEM settings: failed loading <icetable_depth>'
287                write(*,*)'will reconstruct the values of the ice table given the current state'
288                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
289                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
290                do islope = 1,nslope
291                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
292                enddo
293            endif
294            write(*,*) 'PEMETAT0: ICE TABLE done'
295        else if (icetable_dynamic) then
296            call get_field("ice_porefilling",ice_porefilling,found)
297            if (.not. found) write(*,*)'PEM settings: failed loading <ice_porefilling>'
298            call get_field("icetable_depth",icetable_depth,found)
299            if (.not. found) then
300                write(*,*)'PEM settings: failed loading <icetable_depth>'
301                write(*,*)'will reconstruct the values of the ice table given the current state'
302                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
303                do islope = 1,nslope
304                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
305                enddo
306            endif
307            write(*,*) 'PEMETAT0: ICE TABLE done'
308        endif
309
310!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
311!4. CO2 & H2O Adsorption
312        if (adsorption_pem) then
313            do islope = 1,nslope
314                write(num,fmt='(i2.2)') islope
315                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
316                if (.not. found) then
317                    m_co2_regolith_phys = 0.
318                    exit
319                endif
320            enddo
321            do islope = 1,nslope
322                write(num,fmt='(i2.2)') islope
323                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
324                if (.not. found2) then
325                    m_h2o_regolith_phys = 0.
326                    exit
327                endif
328            enddo
329
330            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
331                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
332
333            if (.not. found) deltam_co2_regolith_phys = 0.
334            if (.not. found2) deltam_h2o_regolith_phys = 0.
335
336            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
337        endif ! adsorption_pem
338    endif ! soil_pem
339
340    call close_startphy
341
342else !No startpem, let's build all by hand
343
344    write(*,*)'There is no "startpem.nc"!'
345
346    ! h2o ice
347    h2o_ice = 0.
348    write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
349    do ig = 1,ngrid
350        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
351    enddo
352
353    ! co2 ice
354    write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
355    co2_ice = perennial_co2ice
356
357    ! Stratification (layerings)
358    nb_str_max = 1
359    if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
360
361    if (soil_pem) then
362
363!a) Thermal inertia
364        do islope = 1,nslope
365            do ig = 1,ngrid
366                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
367                    !!! transition
368                    delta = depth_breccia
369                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
370                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
371                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
372                    do iloop=index_breccia + 2,index_bedrock
373                        TI_PEM(ig,iloop,islope) = TI_breccia
374                    enddo
375                else ! we keep the high ti values
376                    do iloop=index_breccia + 1,index_bedrock
377                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
378                    enddo
379                endif
380                !! transition
381                delta = depth_bedrock
382                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
383                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
384                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
385                do iloop = index_bedrock + 2,nsoil_PEM
386                    TI_PEM(ig,iloop,islope) = TI_bedrock
387                enddo
388            enddo
389        enddo
390
391        do iloop = 1,nsoil_PCM
392            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
393        enddo
394        !!! zone de transition
395        delta = depth_breccia
396        do ig = 1,ngrid
397            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
398                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
399                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
400                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
401                do iloop = index_breccia + 2,index_bedrock
402                    inertiedat_PEM(ig,iloop) = TI_breccia
403                enddo
404            else
405                do iloop = index_breccia + 1,index_bedrock
406                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
407                enddo
408            endif
409        enddo
410
411        !!! zone de transition
412        delta = depth_bedrock
413        do ig = 1,ngrid
414            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
415                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
416                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
417        enddo
418
419        do iloop = index_bedrock + 2,nsoil_PEM
420            do ig = 1,ngrid
421                inertiedat_PEM(ig,iloop) = TI_bedrock
422            enddo
423        enddo
424
425        write(*,*) 'PEMETAT0: TI done'
426
427!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428!b) Soil temperature
429        do islope = 1,nslope
430            call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
431            call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
432
433! First raw initialization
434            do isoil = nsoil_PCM + 1,nsoil_PEM
435                do ig = 1,ngrid
436                    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)
437                enddo
438            enddo
439        enddo !islope
440        write(*,*) 'PEMETAT0: TSOIL done'
441
442!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
443!c) Ice table
444        if (icetable_equilibrium) then
445            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
446            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
447            do islope = 1,nslope
448                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
449            enddo
450            write(*,*) 'PEMETAT0: Ice table done'
451        else if (icetable_dynamic) then
452            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
453            do islope = 1,nslope
454                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
455            enddo
456            write(*,*) 'PEMETAT0: Ice table done'
457        endif
458
459!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
460!d) Regolith adsorbed
461        if (adsorption_pem) then
462            m_co2_regolith_phys = 0.
463            m_h2o_regolith_phys = 0.
464            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
465                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
466            deltam_co2_regolith_phys = 0.
467            deltam_h2o_regolith_phys = 0.
468        endif
469
470        write(*,*) 'PEMETAT0: CO2 adsorption done'
471    endif !soil_pem
472endif ! of if (startphy_file)
473
474if (soil_pem) then ! Sanity check
475    do ig = 1,ngrid
476        do islope = 1,nslope
477            do iloop = 1,nsoil_PEM
478                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
479            enddo
480        enddo
481    enddo
482endif !soil_pem
483
484END SUBROUTINE pemetat0
485
486END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.