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

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

PEM:
Few corrections related to r3498 (time step from integer to real) and r3493 (Norbert Schorghofer's subroutines for dynamic ice table) in order to make the code work properly.
JBC
File size: 28.9 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_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,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,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,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,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.