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

Last change on this file since 3137 was 3123, checked in by llange, 2 years ago

MARS PEM
1) Adapting the Mars PEM soil grid to the one of the PCM. The first layers in the PEM follow those from the PCM (PYM grid), and then, for layers at depth, we use the classical power low grid.
2) Correction when reading the soil temperature, water density at the surface, and initialization of the ice table.
LL
File size: 22.5 KB
Line 
1MODULE pemetat0_mod
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness,   &
10                    tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
11                    global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,                &
12                    m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
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, water_reservoir_nom, 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
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
25
26#ifndef CPP_STD
27    use comcstfi_h,   only: r, mugaz
28    use surfdat_h,    only: watercaptag
29#else
30    use comcstfi_mod, only: r, mugaz
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_GCM           ! # of vertical grid points in the GCM
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_ave_pressure ! mean average pressure on the planet [Pa]
45real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
46real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM 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_h2oglaciers    ! tendencies on h2o glaciers
51real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
52real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
53real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
54real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
55! outputs
56real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
57real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [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),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
66real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! 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 yr1[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!!!        1. Thermal Inertia
94!!!        2. Soil Temperature
95!!!        3. Ice table
96!!!        4. Mass of CO2 & H2O adsorbed
97!!!        5. Water reservoir
98!!!
99!!! /!\ This order must be respected !
100!!! Author: LL
101!!!
102!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
103
104!0. Check if the start_PEM exist.
105
106inquire(file=filename,exist =  startpem_file)
107
108write(*,*)'Is start PEM?',startpem_file
109
110!1. Run
111if (startpem_file) then
112    ! open pem initial state file:
113    call open_startphy(filename)
114
115    if (soil_pem) then
116
117!1. Thermal Inertia
118! a. General case
119    do islope = 1,nslope
120        write(num,'(i2.2)') islope
121        call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
122        if (.not. found) then
123            write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
124            write(*,*)'will reconstruct the values of TI_PEM'
125
126            do ig = 1,ngrid
127                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
128                    !!! transition
129                    delta = depth_breccia
130                    TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
131                                                        (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
132                                                        ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
133                    do iloop=index_breccia+2,index_bedrock
134                        TI_PEM(ig,iloop,islope) = TI_breccia
135                    enddo
136                else ! we keep the high ti values
137                    do iloop = index_breccia + 1,index_bedrock
138                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
139                    enddo
140                endif ! TI PEM and breccia comparison
141                !! transition
142                delta = depth_bedrock
143                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
144                                                        (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
145                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
146                do iloop = index_bedrock + 2,nsoil_PEM
147                    TI_PEM(ig,iloop,islope) = TI_bedrock
148                enddo
149            enddo
150        else ! found
151            do iloop = nsoil_GCM + 1,nsoil_PEM
152                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.
153            enddo
154        endif ! not found
155    enddo ! islope
156
157    write(*,*) 'PEMETAT0: THERMAL INERTIA done'
158
159! b. Special case for inertiedat, inertiedat_PEM
160call get_field("inertiedat_PEM",inertiedat_PEM,found)
161if(.not.found) then
162 do iloop = 1,nsoil_GCM
163   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
164 enddo
165!!! zone de transition
166delta = depth_breccia
167do ig = 1,ngrid
168if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
169inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
170            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
171                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
172
173 do iloop = index_breccia+2,index_bedrock
174       inertiedat_PEM(ig,iloop) = TI_breccia
175  enddo
176
177else
178   do iloop=index_breccia+1,index_bedrock
179      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
180   enddo
181endif ! comparison ti breccia
182enddo!ig
183
184!!! zone de transition
185delta = depth_bedrock
186do ig = 1,ngrid
187inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
188            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
189                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
190enddo
191
192 do iloop = index_bedrock+2, nsoil_PEM
193   do ig = 1,ngrid
194      inertiedat_PEM(ig,iloop) = TI_bedrock
195   enddo
196 enddo
197endif ! not found
198
199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200!2. Soil Temperature
201do islope=1,nslope
202  write(num,fmt='(i2.2)') islope
203   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
204  if(.not.found) then
205      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
206      write(*,*)'will reconstruct the values of Tsoil'
207!      do ig = 1,ngrid
208!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
209!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
210!       do iloop=index_breccia+2,index_bedrock
211!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
212!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
213!      enddo
214!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
215!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
216!
217!      do iloop=index_bedrock+2,nsoil_PEM
218!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
219!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
220!      enddo
221!      enddo
222
223
224     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
225     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
226
227
228   else
229! predictor corrector: restart from year 1 of the GCM and build the evolution of
230! tsoil at depth
231
232    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
233    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
234    call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
235    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
236    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
237
238
239     do iloop = nsoil_GCM+1,nsoil_PEM
240       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
241     enddo
242   endif !found
243
244    do it = 1,timelen
245        do isoil = nsoil_GCM+1,nsoil_PEM
246        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
247        enddo
248     enddo
249      do isoil = nsoil_GCM+1,nsoil_PEM
250        do ig = 1,ngrid
251        watersoil_ave(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)
252        enddo
253      enddo
254enddo ! islope
255
256write(*,*) 'PEMETAT0: SOIL TEMP done'
257
258!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259!3. Ice Table
260  call get_field("ice_table",ice_table,found)
261   if(.not.found) then
262      write(*,*)'PEM settings: failed loading <ice_table>'
263      write(*,*)'will reconstruct the values of the ice table given the current state'
264      ice_table(:,:) = -1  ! by default, no ice table 
265      ice_table_thickness(:,:) = -1  ! by default, no ice table
266     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
267     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
268     do islope = 1,nslope
269     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
270     enddo
271   endif
272
273write(*,*) 'PEMETAT0: ICE TABLE done'
274
275!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
276!4. CO2 & H2O Adsorption
277 if(adsorption_pem) then
278  do islope=1,nslope
279   write(num,fmt='(i2.2)') islope
280   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
281    if((.not.found)) then
282       m_co2_regolith_phys(:,:,:) = 0.
283       exit
284    endif
285
286  enddo
287
288  do islope=1,nslope
289   write(num,fmt='(i2.2)') islope
290   call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
291    if((.not.found2)) then
292       m_h2o_regolith_phys(:,:,:) = 0.
293      exit
294    endif
295
296  enddo
297
298    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
299                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
300
301    if((.not.found)) then
302      deltam_co2_regolith_phys(:) = 0.
303    endif
304    if((.not.found2)) then
305       deltam_h2o_regolith_phys(:) = 0.
306    endif
307write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
308    endif
309endif ! soil_pem
310
311!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312!. 5 water reservoir
313#ifndef CPP_STD
314   call get_field("water_reservoir",water_reservoir,found)
315   if(.not.found) then
316      write(*,*)'Pemetat0: failed loading <water_reservoir>'
317      write(*,*)'will reconstruct the values from watercaptag'
318      do ig=1,ngrid
319        if(watercaptag(ig)) then
320           water_reservoir(ig)=water_reservoir_nom
321        else
322           water_reservoir(ig)=0.
323        endif
324      enddo
325   endif
326#endif
327
328call close_startphy
329
330else !No startfi, let's build all by hand
331
332    if (soil_pem) then
333
334!a) Thermal inertia
335   do islope = 1,nslope
336      do ig = 1,ngrid
337
338          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
339!!! transition
340             delta = depth_breccia
341             TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
342            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
343                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
344
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
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
363enddo
364
365 do iloop = 1,nsoil_GCM
366   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
367 enddo
368!!! zone de transition
369delta = depth_breccia
370do ig = 1,ngrid
371      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
372inertiedat_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
376
377 do iloop = index_breccia+2,index_bedrock
378
379       inertiedat_PEM(ig,iloop) = TI_breccia
380
381 enddo
382else
383   do iloop = index_breccia+1,index_bedrock
384       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
385    enddo
386
387endif
388enddo
389
390!!! zone de transition
391delta = depth_bedrock
392do ig = 1,ngrid
393inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
394            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
395                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
396enddo
397
398
399
400 do iloop = index_bedrock+2, nsoil_PEM
401   do ig = 1,ngrid
402      inertiedat_PEM(ig,iloop) = TI_bedrock
403   enddo
404 enddo
405
406write(*,*) 'PEMETAT0: TI done'
407
408!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
409!b) Soil temperature
410    do islope = 1,nslope
411!     do ig = 1,ngrid
412!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
413!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
414!
415!       do iloop=index_breccia+2,index_bedrock
416!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
417!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
418!      enddo
419!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
420!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
421
422!       do iloop=index_bedrock+2,nsoil_PEM
423!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
424!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
425!      enddo
426
427!       enddo
428
429      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
430      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
431 
432! First raw initialization
433      do it = 1,timelen
434        do isoil = nsoil_GCM+1,nsoil_PEM
435          tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
436        enddo
437      enddo
438
439      do it = 1,timelen
440        do isoil = nsoil_GCM+1,nsoil_PEM
441          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
442        enddo
443       enddo
444
445       do isoil = nsoil_GCM+1,nsoil_PEM
446         do ig = 1,ngrid
447           watersoil_ave(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)
448         enddo
449       enddo
450enddo !islope
451write(*,*) 'PEMETAT0: TSOIL done'
452
453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
454!c) Ice table
455       ice_table(:,:) = -1.  ! by default, no ice table
456       ice_table_thickness(:,:) = -1.
457       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
458       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
459       do islope = 1,nslope
460         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
461       enddo
462       write(*,*) 'PEMETAT0: Ice table done'
463
464!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
465!d) Regolith adsorbed
466if (adsorption_pem) then
467    m_co2_regolith_phys(:,:,:) = 0.
468    m_h2o_regolith_phys(:,:,:) = 0.
469
470    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
471                             m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
472
473    deltam_co2_regolith_phys(:) = 0.
474    deltam_h2o_regolith_phys(:) = 0.
475endif
476
477write(*,*) 'PEMETAT0: CO2 adsorption done'
478endif !soil_pem
479
480!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481!. e) water reservoir
482#ifndef CPP_STD
483    do ig = 1,ngrid
484        if (watercaptag(ig)) then
485            water_reservoir(ig) = water_reservoir_nom
486        else
487            water_reservoir(ig) = 0.
488        endif
489    enddo
490#endif
491
492endif ! of if (startphy_file)
493
494if (soil_pem) then ! Sanity check
495    do ig = 1,ngrid
496        do islope = 1,nslope
497            do iloop = 1,nsoil_PEM
498                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
499            enddo
500        enddo
501    enddo
502endif !soil_pem
503
504END SUBROUTINE pemetat0
505
506END MODULE pemetat0_mod
507
Note: See TracBrowser for help on using the repository browser.