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

Last change on this file since 3174 was 3170, checked in by llange, 2 years ago

Mars PEM
Fixing a small bug: the subroutine compute_icetable was always called, even if tthe option 'icetable_equilibrium' was set to false in the run_PEM.def. It is now fixed by adding a flag before the call.
LL
File size: 25.4 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)
13
14use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
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 soil_pem_compute_mod,       only: soil_pem_compute
24use soil_pem_ini_mod,           only: soil_pem_ini
25use comslope_mod,               only: def_slope_mean, subslope_dist
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(len=*),               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]
59real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
60real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
61real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
62real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
63real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
65real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
66real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
67! local
68real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
69real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
70logical                                 :: found                        ! check if variables are found in the start
71logical                                 :: found2                       ! check if variables are found in the start
72integer                                 :: iloop, ig, islope, it, isoil ! 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)        :: tsoil_saved                  ! saved soil temperature [K]
77real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr 1 [K]
78real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
79real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
80real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
81logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
82
83#ifdef CPP_STD
84    logical, dimension(ngrid) :: watercaptag
85    watercaptag = .false.
86#endif
87
88!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
89!!!
90!!! Purpose: read start_pem. Need a specific iostart_PEM
91!!!
92!!! Order: 0. Previous year of the PEM run
93!!!           Ice initialization
94!!!        1. Thermal Inertia
95!!!        2. Soil Temperature
96!!!        3. Ice table
97!!!        4. Mass of CO2 & H2O adsorbed
98!!!
99!!! /!\ This order must be respected !
100!!! Author: LL
101!!!
102!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
103
104!0.1 Check if the start_PEM exist.
105
106inquire(file = filename,exist = startpem_file)
107
108write(*,*)'Is start PEM?',startpem_file
109
110!0.2 Set to default values
111ice_table = -1.  ! by default, no ice table
112ice_table_thickness = -1.
113!1. Run
114if (startpem_file) then
115    ! open pem initial state file:
116    call open_startphy(filename)
117
118!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119#ifndef CPP_STD
120    ! h2o ice
121    h2o_ice = 0.
122    call get_field("h2o_ice",h2o_ice,found)
123    if (.not. found) then
124        write(*,*)'Pemetat0: failed loading <h2o_ice>'
125        write(*,*)'will reconstruct the values from watercaptag'
126        write(*,*)'with default value ''ini_huge_h2oice'''
127        do ig = 1,ngrid
128            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
129        enddo
130    else
131        ! The variations of infinite reservoirs during the PCM years are taken into account
132        h2o_ice = h2o_ice + watercap
133    endif
134
135    ! co2 ice
136    co2_ice = perennial_co2ice
137#endif
138
139    if (soil_pem) then
140
141!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142!1. Thermal Inertia
143! a. General case
144        do islope = 1,nslope
145            write(num,'(i2.2)') islope
146            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
147            if (.not. found) then
148                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
149                write(*,*)'will reconstruct the values of TI_PEM'
150
151                do ig = 1,ngrid
152                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
153                        !!! transition
154                        delta = depth_breccia
155                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
156                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
157                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
158                        do iloop=index_breccia+2,index_bedrock
159                            TI_PEM(ig,iloop,islope) = TI_breccia
160                        enddo
161                    else ! we keep the high TI values
162                        do iloop = index_breccia + 1,index_bedrock
163                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
164                        enddo
165                    endif ! TI PEM and breccia comparison
166                    !! transition
167                    delta = depth_bedrock
168                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
169                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
170                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
171                    do iloop = index_bedrock + 2,nsoil_PEM
172                        TI_PEM(ig,iloop,islope) = TI_bedrock
173                    enddo
174                enddo
175            else ! found
176                do iloop = nsoil_PCM + 1,nsoil_PEM
177                    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.
178                enddo
179            endif ! not found
180        enddo ! islope
181
182        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
183
184! b. Special case for inertiedat, inertiedat_PEM
185        call get_field("inertiedat_PEM",inertiedat_PEM,found)
186        if (.not.found) then
187            do iloop = 1,nsoil_PCM
188                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
189            enddo
190        !!! zone de transition
191            delta = depth_breccia
192            do ig = 1,ngrid
193                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
194                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
195                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
196                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
197
198                    do iloop = index_breccia+2,index_bedrock
199                        inertiedat_PEM(ig,iloop) = TI_breccia
200                    enddo
201                else
202                    do iloop=index_breccia+1,index_bedrock
203                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
204                    enddo
205                endif ! comparison ti breccia
206            enddo ! ig
207
208            !!! zone de transition
209            delta = depth_bedrock
210            do ig = 1,ngrid
211                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
212                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
213                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
214            enddo
215
216            do iloop = index_bedrock + 2, nsoil_PEM
217                do ig = 1,ngrid
218                    inertiedat_PEM(ig,iloop) = TI_bedrock
219                enddo
220            enddo
221        endif ! not found
222
223!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224!2. Soil Temperature
225        do islope=1,nslope
226            write(num,fmt='(i2.2)') islope
227            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
228            if (.not.found) then
229                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
230                write(*,*)'will reconstruct the values of Tsoil'
231!                do ig = 1,ngrid
232!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
233!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
234!                    do iloop=index_breccia+2,index_bedrock
235!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
236!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
237!                    enddo
238!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
239!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
240!
241!                    do iloop=index_bedrock+2,nsoil_PEM
242!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
243!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
244!                    enddo
245!                enddo
246                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
247                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
248            else
249! predictor corrector: restart from year 1 of the PCM and build the evolution of
250! tsoil at depth
251                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
252                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
253                call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
254                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
255                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
256
257                do iloop = nsoil_PCM+1,nsoil_PEM
258                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
259                enddo
260            endif !found
261
262            do it = 1,timelen
263                do isoil = nsoil_PCM+1,nsoil_PEM
264                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
265                enddo
266            enddo
267            do isoil = nsoil_PCM+1,nsoil_PEM
268                do ig = 1,ngrid
269                    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)
270                enddo
271            enddo
272        enddo ! islope
273        write(*,*) 'PEMETAT0: SOIL TEMP done'
274
275!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
276!3. Ice Table
277        if(icetable_equilibrium) then
278            call get_field("ice_table",ice_table,found)
279            if (.not. found) then
280                write(*,*)'PEM settings: failed loading <ice_table>'
281                write(*,*)'will reconstruct the values of the ice table given the current state'
282                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
283                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
284                do islope = 1,nslope
285                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
286                enddo
287            endif
288            write(*,*) 'PEMETAT0: ICE TABLE done'
289        endif
290
291!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
292!4. CO2 & H2O Adsorption
293        if (adsorption_pem) then
294            do islope = 1,nslope
295                write(num,fmt='(i2.2)') islope
296                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
297                if (.not. found) then
298                    m_co2_regolith_phys = 0.
299                    exit
300                endif
301            enddo
302            do islope=1,nslope
303                write(num,fmt='(i2.2)') islope
304                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
305                if (.not. found2) then
306                    m_h2o_regolith_phys = 0.
307                    exit
308                endif
309            enddo
310
311            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, &
312                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
313
314            if (.not. found) deltam_co2_regolith_phys = 0.
315            if (.not.found2) deltam_h2o_regolith_phys = 0.
316
317            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
318        endif ! adsorption_pem
319    endif ! soil_pem
320
321    call close_startphy
322
323else !No startfi, let's build all by hand
324
325    ! h2o ice
326    h2o_ice = 0.
327    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
328        do ig = 1,ngrid
329            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
330        enddo
331
332    ! co2 ice
333    co2_ice = perennial_co2ice
334
335    if (soil_pem) then
336
337!a) Thermal inertia
338        do islope = 1,nslope
339            do ig = 1,ngrid
340                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
341                    !!! transition
342                    delta = depth_breccia
343                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
344                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
345                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
346                    do iloop=index_breccia + 2,index_bedrock
347                        TI_PEM(ig,iloop,islope) = TI_breccia
348                    enddo
349                else ! we keep the high ti values
350                    do iloop=index_breccia + 1,index_bedrock
351                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
352                    enddo
353                endif
354                !! transition
355                delta = depth_bedrock
356                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
357                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
358                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
359                do iloop = index_bedrock + 2,nsoil_PEM
360                    TI_PEM(ig,iloop,islope) = TI_bedrock
361                enddo
362            enddo
363        enddo
364
365        do iloop = 1,nsoil_PCM
366            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
367        enddo
368        !!! zone de transition
369        delta = depth_breccia
370        do ig = 1,ngrid
371            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
372                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
373                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
374                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
375                do iloop = index_breccia + 2,index_bedrock
376                    inertiedat_PEM(ig,iloop) = TI_breccia
377                enddo
378            else
379                do iloop = index_breccia + 1,index_bedrock
380                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
381                enddo
382            endif
383        enddo
384
385        !!! zone de transition
386        delta = depth_bedrock
387        do ig = 1,ngrid
388            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
389                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
390                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
391        enddo
392
393        do iloop = index_bedrock + 2,nsoil_PEM
394            do ig = 1,ngrid
395                inertiedat_PEM(ig,iloop) = TI_bedrock
396            enddo
397        enddo
398
399        write(*,*) 'PEMETAT0: TI done'
400
401!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
402!b) Soil temperature
403        do islope = 1,nslope
404!            do ig = 1,ngrid
405!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
406!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
407!
408!                do iloop=index_breccia+2,index_bedrock
409!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
410!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
411!                enddo
412!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
413!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
414!
415!                do iloop=index_bedrock+2,nsoil_PEM
416!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
417!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
418!                enddo
419!            enddo
420            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
421            call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
422
423! First raw initialization
424            do it = 1,timelen
425                do isoil = nsoil_PCM + 1,nsoil_PEM
426                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
427                enddo
428            enddo
429
430            do it = 1,timelen
431                do isoil = nsoil_PCM + 1,nsoil_PEM
432                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
433                enddo
434            enddo
435
436            do isoil = nsoil_PCM + 1,nsoil_PEM
437                do ig = 1,ngrid
438                    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)
439                enddo
440            enddo
441        enddo !islope
442        write(*,*) 'PEMETAT0: TSOIL done'
443
444!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
445!c) Ice table
446        if(icetable_equilibrium) then
447            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
448            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
449            do islope = 1,nslope
450                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
451            enddo
452            write(*,*) 'PEMETAT0: Ice table done'
453        endif
454       
455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456!d) Regolith adsorbed
457        if (adsorption_pem) then
458            m_co2_regolith_phys = 0.
459            m_h2o_regolith_phys = 0.
460            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, &
461                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
462            deltam_co2_regolith_phys = 0.
463            deltam_h2o_regolith_phys = 0.
464        endif
465
466        write(*,*) 'PEMETAT0: CO2 adsorption done'
467    endif !soil_pem
468endif ! of if (startphy_file)
469
470if (soil_pem) then ! Sanity check
471    do ig = 1,ngrid
472        do islope = 1,nslope
473            do iloop = 1,nsoil_PEM
474                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
475            enddo
476        enddo
477    enddo
478endif !soil_pem
479
480END SUBROUTINE pemetat0
481
482END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.