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

Last change on this file since 3784 was 3778, checked in by jbclement, 10 months ago

PEM:

  • New subroutine to detect whether there is subsurface ice or not
  • Rework of the initialization/update/finalization of the situation regarding the layering data structure
  • Introduction of a threshold 'h_patchy_dust' under which the top dust layer is not considered as a stratum
  • 'deposits' is renamed as 'layerings_map'
  • Few cleanings

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