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

Last change on this file since 3308 was 3308, checked in by jbclement, 7 months ago

PEM:
Few small corrections to make the PEM work in 3D, in particular concerning the initialization of the planet type and the evolution of ice with slopes.
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,ice_table,ice_table_thickness, &
10                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,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
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
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) :: tend_h2o_ice        ! tendencies on h2o ice
52real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! 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           ! Ice table depth [m]
63real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
64real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
65real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
66real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
67real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
68! local
69real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM                  ! soil temperature saved in the start [K]
70real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                     ! soil thermal inertia saved in the start [SI]
71logical                                 :: found, found2                   ! check if variables are found in the start
72integer                                 :: iloop, ig, islope, it, isoil, k ! index for loops
73real                                    :: kcond                           ! Thermal conductivity, intermediate variable [SI]
74real                                    :: delta                           ! Depth of the interface regolith-breccia, breccia -bedrock [m]
75character(2)                            :: num                             ! intermediate string to read PEM start sloped variables
76real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                   ! intermediate soil temperature during yr 1 [K]
77real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                   ! intermediate soil temperature during yr 2 [K]
78logical                                 :: startpem_file                   ! boolean to check if we read the startfile or not
79real, dimension(:,:,:,:), allocatable   :: stratif_array                   ! Array for stratification (layerings)
80real                                    :: nb_str_max
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 start PEM?',startpem_file
108
109!0.2 Set to default values
110ice_table = -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    ! Stratification (layerings)
120    found = inquire_dimension("nb_str_max")
121    if (.not. found) then
122        write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
123        write(*,*) '''nb_str_max'' is set to 10!'
124        nb_str_max = 10.
125    else
126        nb_str_max = inquire_dimension_length('nb_str_max')
127    endif
128    allocate(stratif_array(ngrid,nslope,int(nb_str_max),6))
129    stratif_array = 0.
130    do islope = 1,nslope
131        write(num,'(i2.2)') islope
132        call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
133        call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
134        found = found .or. found2
135        call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
136        found = found .or. found2
137        call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
138        found = found .or. found2
139        call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
140        found = found .or. found2
141        call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
142        found = found .or. found2
143        if (.not. found) then
144            write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
145        endif ! not found
146    enddo ! islope
147    call array2stratif(stratif_array,ngrid,nslope,stratif)
148    deallocate(stratif_array)
149
150!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151#ifndef CPP_STD
152    ! h2o ice
153    h2o_ice = 0.
154    call get_field("h2o_ice",h2o_ice,found)
155    if (.not. found) then
156        write(*,*)'Pemetat0: failed loading <h2o_ice>'
157        write(*,*)'will reconstruct the values from watercaptag'
158        write(*,*)'with default value ''ini_huge_h2oice'''
159        do ig = 1,ngrid
160            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
161        enddo
162    else
163        ! The variations of infinite reservoirs during the PCM years are taken into account
164        h2o_ice = h2o_ice + watercap
165    endif
166
167    ! co2 ice
168    co2_ice = perennial_co2ice
169#endif
170
171    if (soil_pem) then
172
173!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174!1. Thermal Inertia
175! a. General case
176        do islope = 1,nslope
177            write(num,'(i2.2)') islope
178            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
179            if (.not. found) then
180                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
181                write(*,*)'will reconstruct the values of TI_PEM'
182
183                do ig = 1,ngrid
184                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
185                        !!! transition
186                        delta = depth_breccia
187                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
188                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
189                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
190                        do iloop=index_breccia+2,index_bedrock
191                            TI_PEM(ig,iloop,islope) = TI_breccia
192                        enddo
193                    else ! we keep the high TI values
194                        do iloop = index_breccia + 1,index_bedrock
195                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
196                        enddo
197                    endif ! TI PEM and breccia comparison
198                    !! transition
199                    delta = depth_bedrock
200                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
201                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
202                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
203                    do iloop = index_bedrock + 2,nsoil_PEM
204                        TI_PEM(ig,iloop,islope) = TI_bedrock
205                    enddo
206                enddo
207            else ! found
208                do iloop = nsoil_PCM + 1,nsoil_PEM
209                    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.
210                enddo
211            endif ! not found
212        enddo ! islope
213
214        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
215
216! b. Special case for inertiedat, inertiedat_PEM
217        call get_field("inertiedat_PEM",inertiedat_PEM,found)
218        if (.not.found) then
219            do iloop = 1,nsoil_PCM
220                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
221            enddo
222        !!! zone de transition
223            delta = depth_breccia
224            do ig = 1,ngrid
225                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
226                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
227                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
228                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
229
230                    do iloop = index_breccia+2,index_bedrock
231                        inertiedat_PEM(ig,iloop) = TI_breccia
232                    enddo
233                else
234                    do iloop=index_breccia+1,index_bedrock
235                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
236                    enddo
237                endif ! comparison ti breccia
238            enddo ! ig
239
240            !!! zone de transition
241            delta = depth_bedrock
242            do ig = 1,ngrid
243                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
244                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
245                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
246            enddo
247
248            do iloop = index_bedrock + 2, nsoil_PEM
249                do ig = 1,ngrid
250                    inertiedat_PEM(ig,iloop) = TI_bedrock
251                enddo
252            enddo
253        endif ! not found
254
255!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
256!2. Soil Temperature
257        do islope=1,nslope
258            write(num,fmt='(i2.2)') islope
259            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
260            if (.not.found) then
261                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
262                write(*,*)'will reconstruct the values of Tsoil'
263!                do ig = 1,ngrid
264!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
265!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
266!                    do iloop=index_breccia+2,index_bedrock
267!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
268!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
269!                    enddo
270!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
271!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
272!
273!                    do iloop=index_bedrock+2,nsoil_PEM
274!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
275!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
276!                    enddo
277!                enddo
278                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
279                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
280            else
281! predictor corrector: restart from year 1 of the PCM and build the evolution of
282! tsoil at depth
283                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
284                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
285                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
286                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
287                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
288
289                do iloop = nsoil_PCM+1,nsoil_PEM
290                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
291                enddo
292            endif !found
293
294            do it = 1,timelen
295                do isoil = nsoil_PCM+1,nsoil_PEM
296                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
297                enddo
298            enddo
299            do isoil = nsoil_PCM+1,nsoil_PEM
300                do ig = 1,ngrid
301                    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)
302                enddo
303            enddo
304        enddo ! islope
305        write(*,*) 'PEMETAT0: SOIL TEMP done'
306
307!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308!3. Ice Table
309        if(icetable_equilibrium) then
310            call get_field("ice_table",ice_table,found)
311            if (.not. found) then
312                write(*,*)'PEM settings: failed loading <ice_table>'
313                write(*,*)'will reconstruct the values of the ice table given the current state'
314                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
315                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
316                do islope = 1,nslope
317                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
318                enddo
319            endif
320            write(*,*) 'PEMETAT0: ICE TABLE done'
321        endif
322
323!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324!4. CO2 & H2O Adsorption
325        if (adsorption_pem) then
326            do islope = 1,nslope
327                write(num,fmt='(i2.2)') islope
328                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
329                if (.not. found) then
330                    m_co2_regolith_phys = 0.
331                    exit
332                endif
333            enddo
334            do islope=1,nslope
335                write(num,fmt='(i2.2)') islope
336                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
337                if (.not. found2) then
338                    m_h2o_regolith_phys = 0.
339                    exit
340                endif
341            enddo
342
343            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
344                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
345
346            if (.not. found) deltam_co2_regolith_phys = 0.
347            if (.not.found2) deltam_h2o_regolith_phys = 0.
348
349            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
350        endif ! adsorption_pem
351    endif ! soil_pem
352
353    call close_startphy
354
355else !No startfi, let's build all by hand
356
357    write(*,*)'There is no "startpem.nc"!'
358
359    ! Stratification (layerings)
360    write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'
361    nb_str_max = 3
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''.'
372    co2_ice = perennial_co2ice
373
374    if (soil_pem) then
375
376!a) Thermal inertia
377        do islope = 1,nslope
378            do ig = 1,ngrid
379                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
380                    !!! transition
381                    delta = depth_breccia
382                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
383                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
384                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
385                    do iloop=index_breccia + 2,index_bedrock
386                        TI_PEM(ig,iloop,islope) = TI_breccia
387                    enddo
388                else ! we keep the high ti values
389                    do iloop=index_breccia + 1,index_bedrock
390                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
391                    enddo
392                endif
393                !! transition
394                delta = depth_bedrock
395                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
396                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
397                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
398                do iloop = index_bedrock + 2,nsoil_PEM
399                    TI_PEM(ig,iloop,islope) = TI_bedrock
400                enddo
401            enddo
402        enddo
403
404        do iloop = 1,nsoil_PCM
405            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
406        enddo
407        !!! zone de transition
408        delta = depth_breccia
409        do ig = 1,ngrid
410            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
411                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
412                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
413                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
414                do iloop = index_breccia + 2,index_bedrock
415                    inertiedat_PEM(ig,iloop) = TI_breccia
416                enddo
417            else
418                do iloop = index_breccia + 1,index_bedrock
419                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
420                enddo
421            endif
422        enddo
423
424        !!! zone de transition
425        delta = depth_bedrock
426        do ig = 1,ngrid
427            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
428                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
429                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
430        enddo
431
432        do iloop = index_bedrock + 2,nsoil_PEM
433            do ig = 1,ngrid
434                inertiedat_PEM(ig,iloop) = TI_bedrock
435            enddo
436        enddo
437
438        write(*,*) 'PEMETAT0: TI done'
439
440!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
441!b) Soil temperature
442        do islope = 1,nslope
443!            do ig = 1,ngrid
444!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
445!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
446!
447!                do iloop=index_breccia+2,index_bedrock
448!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
449!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
450!                enddo
451!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
452!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
453!
454!                do iloop=index_bedrock+2,nsoil_PEM
455!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
456!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
457!                enddo
458!            enddo
459            call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
460            call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
461
462! First raw initialization
463            do it = 1,timelen
464                do isoil = nsoil_PCM + 1,nsoil_PEM
465                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
466                enddo
467            enddo
468
469            do it = 1,timelen
470                do isoil = nsoil_PCM + 1,nsoil_PEM
471                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
472                enddo
473            enddo
474
475            do isoil = nsoil_PCM + 1,nsoil_PEM
476                do ig = 1,ngrid
477                    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)
478                enddo
479            enddo
480        enddo !islope
481        write(*,*) 'PEMETAT0: TSOIL done'
482
483!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
484!c) Ice table
485        if(icetable_equilibrium) then
486            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
487            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
488            do islope = 1,nslope
489                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
490            enddo
491            write(*,*) 'PEMETAT0: Ice table done'
492        endif
493       
494!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495!d) Regolith adsorbed
496        if (adsorption_pem) then
497            m_co2_regolith_phys = 0.
498            m_h2o_regolith_phys = 0.
499            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
500                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
501            deltam_co2_regolith_phys = 0.
502            deltam_h2o_regolith_phys = 0.
503        endif
504
505        write(*,*) 'PEMETAT0: CO2 adsorption done'
506    endif !soil_pem
507endif ! of if (startphy_file)
508
509if (soil_pem) then ! Sanity check
510    do ig = 1,ngrid
511        do islope = 1,nslope
512            do iloop = 1,nsoil_PEM
513                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
514            enddo
515        enddo
516    enddo
517endif !soil_pem
518
519END SUBROUTINE pemetat0
520
521END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.