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

Last change on this file since 3525 was 3525, checked in by jbclement, 6 days ago

PEM:
Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
JBC
File size: 29.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_depth,ice_table_thickness, &
10                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,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,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, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, 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, icetable_depth, icetable_thickness, 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 comslope_mod,               only: def_slope_mean, subslope_dist
25use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo
26
27#ifndef CPP_STD
28    use comcstfi_h,   only: r, mugaz, pi
29    use surfdat_h,    only: watercaptag, watercap, perennial_co2ice
30#else
31    use comcstfi_mod, only: r, mugaz, pi
32#endif
33
34implicit none
35
36include "callkeys.h"
37
38character(*),                   intent(in) :: filename            ! name of the startfi_PEM.nc
39integer,                        intent(in) :: ngrid               ! # of physical grid points
40integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
41integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
42integer,                        intent(in) :: nslope              ! # of sub-grid slopes
43integer,                        intent(in) :: timelen             ! # time samples
44real,                           intent(in) :: timestep            ! time step [s]
45real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
46real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
47real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
48real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
49real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
50real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
51real, dimension(ngrid,nslope),  intent(in) :: d_h2oice        ! tendencies on h2o ice
52real, dimension(ngrid,nslope),  intent(in) :: d_co2ice        ! tendencies on co2 ice
53real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
54! outputs
55real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
56real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
57real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
58real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
59type(layering), dimension(ngrid,nslope),         intent(inout) :: stratif             ! stratification (layerings)
60real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
61real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
62real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_depth     ! Ice table depth [m]
63real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
65real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
66real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
67real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
68real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
69! local
70real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM                  ! soil temperature saved in the start [K]
71real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                     ! soil thermal inertia saved in the start [SI]
72logical                                 :: found, found2                   ! check if variables are found in the start
73integer                                 :: iloop, ig, islope, it, isoil, k ! index for loops
74real                                    :: kcond                           ! Thermal conductivity, intermediate variable [SI]
75real                                    :: delta                           ! Depth of the interface regolith-breccia, breccia -bedrock [m]
76character(2)                            :: num                             ! intermediate string to read PEM start sloped variables
77real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                   ! intermediate soil temperature during yr 1 [K]
78real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                   ! intermediate soil temperature during yr 2 [K]
79logical                                 :: startpem_file                   ! boolean to check if we read the startfile or not
80real, dimension(:,:,:,:), allocatable   :: stratif_array                   ! Array for stratification (layerings)
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.1 Check if the start_PEM exist.
104
105inquire(file = filename,exist = startpem_file)
106
107write(*,*)'Is there a "startpem" file?',startpem_file
108
109!0.2 Set to default values
110ice_table_depth = -1. ! by default, no ice table
111ice_table_thickness = -1.
112!1. Run
113if (startpem_file) then
114    write(*,*)'There is a "startpem.nc"!'
115    ! open pem initial state file:
116    call open_startphy(filename)
117
118!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119#ifndef CPP_STD
120    ! h2o ice
121    h2o_ice = 0.
122    call get_field("h2o_ice",h2o_ice,found)
123    if (.not. found) then
124        write(*,*)'Pemetat0: failed loading <h2o_ice>'
125        write(*,*)'will reconstruct the values from watercaptag'
126        write(*,*)'with default value ''ini_huge_h2oice'''
127        do ig = 1,ngrid
128            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
129        enddo
130    else
131        ! The variations of infinite reservoirs during the PCM years are taken into account
132        h2o_ice = h2o_ice + watercap
133    endif
134
135    ! co2 ice
136    co2_ice = perennial_co2ice
137#endif
138
139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140    ! Stratification (layerings)
141    nb_str_max = 1
142    if (layering_algo) then
143        found = inquire_dimension("nb_str_max")
144        if (.not. found) then
145            write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
146            write(*,*) '''nb_str_max'' is set to 3!'
147        else
148            nb_str_max = int(inquire_dimension_length('nb_str_max'))
149        endif
150        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
151        stratif_array = 0.
152        do islope = 1,nslope
153            write(num,'(i2.2)') islope
154            call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
155            call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
156            found = found .or. found2
157            call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
158            found = found .or. found2
159            call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
160            found = found .or. found2
161            call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
162            found = found .or. found2
163            call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
164            found = found .or. found2
165            if (.not. found) then
166                write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
167            endif ! not found
168        enddo ! islope
169        call array2stratif(stratif_array,ngrid,nslope,stratif)
170        deallocate(stratif_array)
171    endif
172
173    if (soil_pem) then
174
175!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176!1. Thermal Inertia
177! a. General case
178        do islope = 1,nslope
179            write(num,'(i2.2)') islope
180            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
181            if (.not. found) then
182                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
183                write(*,*)'will reconstruct the values of TI_PEM'
184
185                do ig = 1,ngrid
186                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
187                        !!! transition
188                        delta = depth_breccia
189                        TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
190                                                                  (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
191                                                                  ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
192                        do iloop=index_breccia + 2,index_bedrock
193                            TI_PEM(ig,iloop,islope) = TI_breccia
194                        enddo
195                    else ! we keep the high TI values
196                        do iloop = index_breccia + 1,index_bedrock
197                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
198                        enddo
199                    endif ! TI PEM and breccia comparison
200                    !! transition
201                    delta = depth_bedrock
202                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
203                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
204                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
205                    do iloop = index_bedrock + 2,nsoil_PEM
206                        TI_PEM(ig,iloop,islope) = TI_bedrock
207                    enddo
208                enddo
209            else ! found
210                do iloop = nsoil_PCM + 1,nsoil_PEM
211                    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.
212                enddo
213            endif ! not found
214        enddo ! islope
215
216        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
217
218! b. Special case for inertiedat, inertiedat_PEM
219        call get_field("inertiedat_PEM",inertiedat_PEM,found)
220        if (.not. found) then
221            do iloop = 1,nsoil_PCM
222                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
223            enddo
224        !!! zone de transition
225            delta = depth_breccia
226            do ig = 1,ngrid
227                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
228                    inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
229                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
230                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
231
232                    do iloop = index_breccia + 2,index_bedrock
233                        inertiedat_PEM(ig,iloop) = TI_breccia
234                    enddo
235                else
236                    do iloop=index_breccia + 1,index_bedrock
237                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
238                    enddo
239                endif ! comparison ti breccia
240            enddo ! ig
241
242            !!! zone de transition
243            delta = depth_bedrock
244            do ig = 1,ngrid
245                inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
246                                                           (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
247                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
248            enddo
249
250            do iloop = index_bedrock + 2, nsoil_PEM
251                do ig = 1,ngrid
252                    inertiedat_PEM(ig,iloop) = TI_bedrock
253                enddo
254            enddo
255        endif ! not found
256
257!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258!2. Soil Temperature
259        do islope=1,nslope
260            write(num,fmt='(i2.2)') islope
261            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
262            if (.not. found) then
263                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
264                write(*,*)'will reconstruct the values of Tsoil'
265!                do ig = 1,ngrid
266!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
267!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
268!                    do iloop=index_breccia+2,index_bedrock
269!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
270!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
271!                    enddo
272!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
273!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
274!
275!                    do iloop=index_bedrock+2,nsoil_PEM
276!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
277!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
278!                    enddo
279!                enddo
280                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
281                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
282            else
283! predictor corrector: restart from year 1 of the PCM and build the evolution of
284! tsoil at depth
285                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
286                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
287                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
288                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
289                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
290
291                do iloop = nsoil_PCM+1,nsoil_PEM
292                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
293                enddo
294            endif !found
295
296            do it = 1,timelen
297                do isoil = nsoil_PCM+1,nsoil_PEM
298                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
299                enddo
300            enddo
301            do isoil = nsoil_PCM+1,nsoil_PEM
302                do ig = 1,ngrid
303                    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)
304                enddo
305            enddo
306        enddo ! islope
307        write(*,*) 'PEMETAT0: SOIL TEMP done'
308
309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310!3. Ice Table
311        if (icetable_equilibrium) then
312            call get_field("ice_table_depth",ice_table_depth,found)
313            if (.not. found) then
314                write(*,*)'PEM settings: failed loading <ice_table_depth>'
315                write(*,*)'will reconstruct the values of the ice table given the current state'
316                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
317                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
318                do islope = 1,nslope
319                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
320                enddo
321            endif
322            write(*,*) 'PEMETAT0: ICE TABLE done'
323        else if (icetable_dynamic) then
324            call get_field("ice_table_depth",ice_table_depth,found)
325            if (.not. found) then
326                write(*,*)'PEM settings: failed loading <ice_table_depth>'
327                write(*,*)'will reconstruct the values of the ice table given the current state'
328                ice_table_depth = -9999.
329                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
330                do islope = 1,nslope
331                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
332                enddo
333            endif
334            call get_field("ice_porefilling",ice_porefilling,found)
335            if (.not. found) then
336                write(*,*)'PEM settings: failed loading <ice_porefilling>'
337                ice_porefilling = 1.
338            endif
339            write(*,*) 'PEMETAT0: ICE TABLE done'
340        endif
341
342!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
343!4. CO2 & H2O Adsorption
344        if (adsorption_pem) then
345            do islope = 1,nslope
346                write(num,fmt='(i2.2)') islope
347                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
348                if (.not. found) then
349                    m_co2_regolith_phys = 0.
350                    exit
351                endif
352            enddo
353            do islope=1,nslope
354                write(num,fmt='(i2.2)') islope
355                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
356                if (.not. found2) then
357                    m_h2o_regolith_phys = 0.
358                    exit
359                endif
360            enddo
361
362            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
363                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
364
365            if (.not. found) deltam_co2_regolith_phys = 0.
366            if (.not. found2) deltam_h2o_regolith_phys = 0.
367
368            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
369        endif ! adsorption_pem
370    endif ! soil_pem
371
372    call close_startphy
373
374else !No startpem, let's build all by hand
375
376    write(*,*)'There is no "startpem.nc"!'
377
378    ! h2o ice
379    h2o_ice = 0.
380    write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
381    do ig = 1,ngrid
382        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
383    enddo
384
385    ! co2 ice
386    write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
387    co2_ice = perennial_co2ice
388
389    ! Stratification (layerings)
390    nb_str_max = 1
391    if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
392
393    if (soil_pem) then
394
395!a) Thermal inertia
396        do islope = 1,nslope
397            do ig = 1,ngrid
398                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
399                    !!! transition
400                    delta = depth_breccia
401                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
402                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
403                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
404                    do iloop=index_breccia + 2,index_bedrock
405                        TI_PEM(ig,iloop,islope) = TI_breccia
406                    enddo
407                else ! we keep the high ti values
408                    do iloop=index_breccia + 1,index_bedrock
409                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
410                    enddo
411                endif
412                !! transition
413                delta = depth_bedrock
414                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
415                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
416                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
417                do iloop = index_bedrock + 2,nsoil_PEM
418                    TI_PEM(ig,iloop,islope) = TI_bedrock
419                enddo
420            enddo
421        enddo
422
423        do iloop = 1,nsoil_PCM
424            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
425        enddo
426        !!! zone de transition
427        delta = depth_breccia
428        do ig = 1,ngrid
429            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
430                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
431                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
432                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
433                do iloop = index_breccia + 2,index_bedrock
434                    inertiedat_PEM(ig,iloop) = TI_breccia
435                enddo
436            else
437                do iloop = index_breccia + 1,index_bedrock
438                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
439                enddo
440            endif
441        enddo
442
443        !!! zone de transition
444        delta = depth_bedrock
445        do ig = 1,ngrid
446            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
447                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
448                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
449        enddo
450
451        do iloop = index_bedrock + 2,nsoil_PEM
452            do ig = 1,ngrid
453                inertiedat_PEM(ig,iloop) = TI_bedrock
454            enddo
455        enddo
456
457        write(*,*) 'PEMETAT0: TI done'
458
459!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
460!b) Soil temperature
461        do islope = 1,nslope
462!            do ig = 1,ngrid
463!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
464!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
465!
466!                do iloop=index_breccia+2,index_bedrock
467!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
468!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
469!                enddo
470!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
471!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
472!
473!                do iloop=index_bedrock+2,nsoil_PEM
474!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
475!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
476!                enddo
477!            enddo
478            call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
479            call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
480
481! First raw initialization
482            do it = 1,timelen
483                do isoil = nsoil_PCM + 1,nsoil_PEM
484                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
485                enddo
486            enddo
487
488            do it = 1,timelen
489                do isoil = nsoil_PCM + 1,nsoil_PEM
490                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
491                enddo
492            enddo
493
494            do isoil = nsoil_PCM + 1,nsoil_PEM
495                do ig = 1,ngrid
496                    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)
497                enddo
498            enddo
499        enddo !islope
500        write(*,*) 'PEMETAT0: TSOIL done'
501
502!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503!c) Ice table
504        if (icetable_equilibrium) then
505            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
506            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
507            do islope = 1,nslope
508                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
509            enddo
510            write(*,*) 'PEMETAT0: Ice table done'
511        else if (icetable_dynamic) then
512            ice_porefilling = 1.
513            ice_table_depth = -9999.
514            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
515            do islope = 1,nslope
516                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
517            enddo
518            write(*,*) 'PEMETAT0: Ice table done'
519        endif
520       
521!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522!d) Regolith adsorbed
523        if (adsorption_pem) then
524            m_co2_regolith_phys = 0.
525            m_h2o_regolith_phys = 0.
526            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
527                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
528            deltam_co2_regolith_phys = 0.
529            deltam_h2o_regolith_phys = 0.
530        endif
531
532        write(*,*) 'PEMETAT0: CO2 adsorption done'
533    endif !soil_pem
534endif ! of if (startphy_file)
535
536if (soil_pem) then ! Sanity check
537    do ig = 1,ngrid
538        do islope = 1,nslope
539            do iloop = 1,nsoil_PEM
540                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
541            enddo
542        enddo
543    enddo
544endif !soil_pem
545
546END SUBROUTINE pemetat0
547
548END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.