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

Last change on this file since 3980 was 3961, checked in by jbclement, 2 weeks ago

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