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

Last change on this file since 3190 was 3190, checked in by llange, 10 months ago

PEM
Fixing bug in writediagpem: soil layers written in the diagpem where those of the PCM and not the PEM.
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 compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
24use comslope_mod,               only: def_slope_mean, subslope_dist
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
35include "callkeys.h"
36
37character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
38integer,                        intent(in) :: ngrid               ! # of physical grid points
39integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
40integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
41integer,                        intent(in) :: nslope              ! # of sub-grid slopes
42integer,                        intent(in) :: timelen             ! # time samples
43real,                           intent(in) :: timestep            ! time step [s]
44real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
45real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
46real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
47real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
48real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
49real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
50real, dimension(ngrid,nslope),  intent(in) :: tend_h2o_ice        ! tendencies on h2o ice
51real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! tendencies on co2 ice
52real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
53! outputs
54real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
55real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
56real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
57real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
58real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
59real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
60real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
61real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
62real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
63real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
65real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
66! local
67real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
68real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
69logical                                 :: found                        ! check if variables are found in the start
70logical                                 :: found2                       ! check if variables are found in the start
71integer                                 :: iloop, ig, islope, it, isoil ! index for loops
72real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
73real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
74character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
75real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
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]
78real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
79real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
80logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
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    ! open pem initial state file:
115    call open_startphy(filename)
116
117!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118#ifndef CPP_STD
119    ! h2o ice
120    h2o_ice = 0.
121    call get_field("h2o_ice",h2o_ice,found)
122    if (.not. found) then
123        write(*,*)'Pemetat0: failed loading <h2o_ice>'
124        write(*,*)'will reconstruct the values from watercaptag'
125        write(*,*)'with default value ''ini_huge_h2oice'''
126        do ig = 1,ngrid
127            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
128        enddo
129    else
130        ! The variations of infinite reservoirs during the PCM years are taken into account
131        h2o_ice = h2o_ice + watercap
132    endif
133
134    ! co2 ice
135    co2_ice = perennial_co2ice
136#endif
137
138    if (soil_pem) then
139
140!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141!1. Thermal Inertia
142! a. General case
143        do islope = 1,nslope
144            write(num,'(i2.2)') islope
145            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
146            if (.not. found) then
147                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
148                write(*,*)'will reconstruct the values of TI_PEM'
149
150                do ig = 1,ngrid
151                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
152                        !!! transition
153                        delta = depth_breccia
154                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
155                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
156                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
157                        do iloop=index_breccia+2,index_bedrock
158                            TI_PEM(ig,iloop,islope) = TI_breccia
159                        enddo
160                    else ! we keep the high TI values
161                        do iloop = index_breccia + 1,index_bedrock
162                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
163                        enddo
164                    endif ! TI PEM and breccia comparison
165                    !! transition
166                    delta = depth_bedrock
167                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
168                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
169                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
170                    do iloop = index_bedrock + 2,nsoil_PEM
171                        TI_PEM(ig,iloop,islope) = TI_bedrock
172                    enddo
173                enddo
174            else ! found
175                do iloop = nsoil_PCM + 1,nsoil_PEM
176                    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.
177                enddo
178            endif ! not found
179        enddo ! islope
180
181        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
182
183! b. Special case for inertiedat, inertiedat_PEM
184        call get_field("inertiedat_PEM",inertiedat_PEM,found)
185        if (.not.found) then
186            do iloop = 1,nsoil_PCM
187                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
188            enddo
189        !!! zone de transition
190            delta = depth_breccia
191            do ig = 1,ngrid
192                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
193                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
194                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
195                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
196
197                    do iloop = index_breccia+2,index_bedrock
198                        inertiedat_PEM(ig,iloop) = TI_breccia
199                    enddo
200                else
201                    do iloop=index_breccia+1,index_bedrock
202                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
203                    enddo
204                endif ! comparison ti breccia
205            enddo ! ig
206
207            !!! zone de transition
208            delta = depth_bedrock
209            do ig = 1,ngrid
210                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
211                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
212                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
213            enddo
214
215            do iloop = index_bedrock + 2, nsoil_PEM
216                do ig = 1,ngrid
217                    inertiedat_PEM(ig,iloop) = TI_bedrock
218                enddo
219            enddo
220        endif ! not found
221
222!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
223!2. Soil Temperature
224        do islope=1,nslope
225            write(num,fmt='(i2.2)') islope
226            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
227            if (.not.found) then
228                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
229                write(*,*)'will reconstruct the values of Tsoil'
230!                do ig = 1,ngrid
231!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
232!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
233!                    do iloop=index_breccia+2,index_bedrock
234!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
235!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
236!                    enddo
237!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
238!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
239!
240!                    do iloop=index_bedrock+2,nsoil_PEM
241!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
242!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
243!                    enddo
244!                enddo
245                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
246                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
247            else
248! predictor corrector: restart from year 1 of the PCM and build the evolution of
249! tsoil at depth
250                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
251                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
252                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
253                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
254                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
255               
256                do iloop = nsoil_PCM+1,nsoil_PEM
257                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
258                enddo
259            endif !found
260
261            do it = 1,timelen
262                do isoil = nsoil_PCM+1,nsoil_PEM
263                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
264                enddo
265            enddo
266            do isoil = nsoil_PCM+1,nsoil_PEM
267                do ig = 1,ngrid
268                    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)
269                enddo
270            enddo
271        enddo ! islope
272        write(*,*) 'PEMETAT0: SOIL TEMP done'
273
274!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275!3. Ice Table
276        if(icetable_equilibrium) then
277            call get_field("ice_table",ice_table,found)
278            if (.not. found) then
279                write(*,*)'PEM settings: failed loading <ice_table>'
280                write(*,*)'will reconstruct the values of the ice table given the current state'
281                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
282                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
283                do islope = 1,nslope
284                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
285                enddo
286            endif
287            write(*,*) 'PEMETAT0: ICE TABLE done'
288        endif
289
290!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291!4. CO2 & H2O Adsorption
292        if (adsorption_pem) then
293            do islope = 1,nslope
294                write(num,fmt='(i2.2)') islope
295                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
296                if (.not. found) then
297                    m_co2_regolith_phys = 0.
298                    exit
299                endif
300            enddo
301            do islope=1,nslope
302                write(num,fmt='(i2.2)') islope
303                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
304                if (.not. found2) then
305                    m_h2o_regolith_phys = 0.
306                    exit
307                endif
308            enddo
309
310            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, &
311                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
312
313            if (.not. found) deltam_co2_regolith_phys = 0.
314            if (.not.found2) deltam_h2o_regolith_phys = 0.
315
316            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
317        endif ! adsorption_pem
318    endif ! soil_pem
319
320    call close_startphy
321
322else !No startfi, let's build all by hand
323
324    ! h2o ice
325    h2o_ice = 0.
326    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
327        do ig = 1,ngrid
328            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
329        enddo
330
331    ! co2 ice
332    co2_ice = perennial_co2ice
333
334    if (soil_pem) then
335
336!a) Thermal inertia
337        do islope = 1,nslope
338            do ig = 1,ngrid
339                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
340                    !!! transition
341                    delta = depth_breccia
342                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
343                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
344                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
345                    do iloop=index_breccia + 2,index_bedrock
346                        TI_PEM(ig,iloop,islope) = TI_breccia
347                    enddo
348                else ! we keep the high ti values
349                    do iloop=index_breccia + 1,index_bedrock
350                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
351                    enddo
352                endif
353                !! transition
354                delta = depth_bedrock
355                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
356                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
357                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
358                do iloop = index_bedrock + 2,nsoil_PEM
359                    TI_PEM(ig,iloop,islope) = TI_bedrock
360                enddo
361            enddo
362        enddo
363
364        do iloop = 1,nsoil_PCM
365            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
366        enddo
367        !!! zone de transition
368        delta = depth_breccia
369        do ig = 1,ngrid
370            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
371                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
372                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
373                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
374                do iloop = index_breccia + 2,index_bedrock
375                    inertiedat_PEM(ig,iloop) = TI_breccia
376                enddo
377            else
378                do iloop = index_breccia + 1,index_bedrock
379                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
380                enddo
381            endif
382        enddo
383
384        !!! zone de transition
385        delta = depth_bedrock
386        do ig = 1,ngrid
387            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
388                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
389                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
390        enddo
391
392        do iloop = index_bedrock + 2,nsoil_PEM
393            do ig = 1,ngrid
394                inertiedat_PEM(ig,iloop) = TI_bedrock
395            enddo
396        enddo
397
398        write(*,*) 'PEMETAT0: TI done'
399
400!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
401!b) Soil temperature
402        do islope = 1,nslope
403!            do ig = 1,ngrid
404!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
405!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
406!
407!                do iloop=index_breccia+2,index_bedrock
408!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
409!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
410!                enddo
411!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
412!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
413!
414!                do iloop=index_bedrock+2,nsoil_PEM
415!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
416!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
417!                enddo
418!            enddo
419            call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
420            call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
421
422! First raw initialization
423            do it = 1,timelen
424                do isoil = nsoil_PCM + 1,nsoil_PEM
425                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
426                enddo
427            enddo
428
429            do it = 1,timelen
430                do isoil = nsoil_PCM + 1,nsoil_PEM
431                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
432                enddo
433            enddo
434
435            do isoil = nsoil_PCM + 1,nsoil_PEM
436                do ig = 1,ngrid
437                    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)
438                enddo
439            enddo
440        enddo !islope
441        write(*,*) 'PEMETAT0: TSOIL done'
442
443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
444!c) Ice table
445        if(icetable_equilibrium) then
446            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
447            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
448            do islope = 1,nslope
449                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
450            enddo
451            write(*,*) 'PEMETAT0: Ice table done'
452        endif
453       
454!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
455!d) Regolith adsorbed
456        if (adsorption_pem) then
457            m_co2_regolith_phys = 0.
458            m_h2o_regolith_phys = 0.
459            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, &
460                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
461            deltam_co2_regolith_phys = 0.
462            deltam_h2o_regolith_phys = 0.
463        endif
464
465        write(*,*) 'PEMETAT0: CO2 adsorption done'
466    endif !soil_pem
467endif ! of if (startphy_file)
468
469if (soil_pem) then ! Sanity check
470    do ig = 1,ngrid
471        do islope = 1,nslope
472            do iloop = 1,nsoil_PEM
473                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
474            enddo
475        enddo
476    enddo
477endif !soil_pem
478
479END SUBROUTINE pemetat0
480
481END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.