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

Last change on this file since 3984 was 3984, checked in by jbclement, 10 days ago

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