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

Last change on this file since 2962 was 2961, checked in by llange, 19 months ago

PEM

  • Move math functions (findroot, deriv, etc) in a specific module
  • Ice table is no longer infinite, the bottom of the ice table (zend) is set such at d(rho_vap) dz (zend) = 0. Future commit will transfer the loss of ice table to perenial glaciers

LL
File size: 22.2 KB
Line 
1   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, ice_table_thickness, &
2                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
3                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
4                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
5   
6   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
7   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
8   use comsoil_h, only:  volcapa,inertiedat
9   use adsorption_mod, only : regolith_adsorption,adsorption_pem
10   USE temps_mod_evol, ONLY: year_PEM
11   USE ice_table_mod, only: computeice_table_equilibrium
12   USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
13                                   TI_breccia,TI_bedrock
14   use soil_thermalproperties_mod, only: update_soil_thermalproperties
15#ifndef CPP_STD   
16   use surfdat_h, only: watercaptag
17#endif
18
19  implicit none
20 
21  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
22  integer,intent(in) :: ngrid                                 ! # of physical grid points
23  integer,intent(in) :: nsoil_GCM                             ! # of vertical grid points in the GCM
24  integer,intent(in) :: nsoil_PEM                             ! # of vertical grid points in the PEM
25  integer,intent(in) :: nslope                                ! # of sub-grid slopes
26  integer,intent(in) :: timelen                               ! # time samples
27  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)             ! surface temperature at the first year of GCM call [K]
28  real,intent(in) :: tsurf_ave_yr2(ngrid,nslope)              ! surface temperature at the second  year of GCM call [K]
29  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
30  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
31  real,intent(in) :: ps_inst(ngrid,timelen)                   ! surface pressure [Pa]
32  real,intent(in) :: timestep                                 ! time step [s]
33  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)           ! tendencies on h2o glaciers
34  real,intent(in) :: tend_co2glaciers(ngrid,nslope)           ! tendencies on co2 glaciers
35  real,intent(in) :: co2ice(ngrid,nslope)                     ! co2 ice amount [kg/m^2]
36  real,intent(in) :: waterice(ngrid,nslope)                   ! water ice amount [kg/m^2]
37  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)   ! soil temperature during 1st year [K]
38  real, intent(in) :: watersurf_ave(ngrid,nslope)             ! surface water ice density, yearly averaged  (kg/m^3)
39  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3)
40  real, intent(in) :: global_ave_pressure                     ! mean average pressure on the planet [Pa]
41! outputs
42
43  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)                ! soil (mid-layer) thermal inertia in the PEM grid [SI]
44  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)             ! soil (mid-layer) temperature [K]
45  real,intent(inout) :: ice_table(ngrid,nslope)                       ! Ice table depth [m]
46  real,intent(inout) :: ice_table_thickness(ngrid,nslope)             ! Ice table thickness [m]
47  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen)    ! instantaneous soil (mid-layer) temperature [K]
48  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
49  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
50  real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of h2o adsorbed [kg/m^2]
51  real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
52  real,intent(inout) :: water_reservoir(ngrid)                        ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
53! local
54   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
55   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)                        ! soil thermal inertia saved in the start [SI]
56   LOGICAL :: found                                                   ! check if variables are found in the start
57   LOGICAL :: found2                                                  ! check if variables are found in the start
58   integer :: iloop,ig,islope,it,isoil                                ! index for loops
59   real :: kcond                                                      ! Thermal conductivity, intermediate variable [SI]
60   real :: delta                                                      ! Depth of the interface regolith-breccia, breccia -bedrock [m]
61   CHARACTER*2 :: num                                                 ! intermediate string to read PEM start sloped variables
62   real :: tsoil_saved(ngrid,nsoil_PEM)                               ! saved soil temperature [K]
63   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr1[K]
64   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr 2 [K]
65   real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
66   real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
67   real :: year_PEM_read                                              ! Year of the PEM previous run
68   LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
69
70!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
71!!!
72!!! Purpose: read start_pem. Need a specific iostart_PEM
73!!!
74!!! Order: 0. Previous year of the PEM run
75!!!        1. Thermal Inertia
76!!!        2. Soil Temperature
77!!!        3. Ice table
78!!!        4. Mass of CO2 & H2O adsorbed
79!!!        5. Water reservoir
80!!!
81!!! /!\ This order must be respected !
82!!! Author: LL
83!!!
84!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85
86
87!0. Check if the start_PEM exist.
88
89inquire(file=filename,exist =  startpem_file)
90
91write(*,*)'Is start PEM?',startpem_file
92
93!1. Run
94
95if (startpem_file) then
96   ! open pem initial state file:
97   call open_startphy(filename)
98
99   call get_var("Time",year_PEM_read,found)
100   year_PEM=INT(year_PEM_read)
101   if(.not.found) then
102      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
103   else
104      write(*,*)'year_PEM of startpem=', year_PEM
105   endif
106
107    if(soil_pem) then
108 
109!1. Thermal Inertia
110! a. General case
111
112DO islope=1,nslope
113  write(num,fmt='(i2.2)') islope
114  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
115   if(.not.found) then
116      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
117      write(*,*)'will reconstruct the values of TI_PEM'
118
119      do ig = 1,ngrid
120          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
121!!! transition
122             delta = depth_breccia
123             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
124            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
125                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
126            do iloop=index_breccia+2,index_bedrock
127              TI_PEM(ig,iloop,islope) = TI_breccia
128            enddo     
129          else ! we keep the high ti values
130            do iloop=index_breccia+1,index_bedrock
131              TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
132            enddo
133          endif ! TI PEM and breccia comparison
134!! transition
135          delta = depth_bedrock
136          TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
137               (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
138               ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
139          do iloop=index_bedrock+2,nsoil_PEM
140            TI_PEM(ig,iloop,islope) = TI_bedrock
141          enddo   
142      enddo
143   else ! found
144     do iloop = nsoil_GCM+1,nsoil_PEM
145       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.
146     enddo
147   endif ! not found
148ENDDO ! islope
149
150print *,'PEMETAT0: THERMAL INERTIA DONE'
151
152! b. Special case for inertiedat, inertiedat_PEM
153call get_field("inertiedat_PEM",inertiedat_PEM,found)
154if(.not.found) then
155 do iloop = 1,nsoil_GCM
156   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
157 enddo
158!!! zone de transition
159delta = depth_breccia
160do ig = 1,ngrid
161if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
162inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
163            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
164                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
165
166 do iloop = index_breccia+2,index_bedrock
167       inertiedat_PEM(ig,iloop) = TI_breccia
168  enddo
169
170else
171   do iloop=index_breccia+1,index_bedrock
172      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
173   enddo
174endif ! comparison ti breccia
175enddo!ig
176
177!!! zone de transition
178delta = depth_bedrock
179do ig = 1,ngrid
180inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
181            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
182                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
183enddo
184
185 do iloop = index_bedrock+2, nsoil_PEM
186   do ig = 1,ngrid
187      inertiedat_PEM(ig,iloop) = TI_bedrock
188   enddo
189 enddo
190endif ! not found
191
192!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193
194!2. Soil Temperature
195
196DO islope=1,nslope
197  write(num,fmt='(i2.2)') islope
198   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
199  if(.not.found) then
200      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
201      write(*,*)'will reconstruct the values of Tsoil'
202!      do ig = 1,ngrid
203!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
204!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))   
205!       do iloop=index_breccia+2,index_bedrock
206!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
207!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
208!      enddo
209!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
210!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1)) 
211!
212!      do iloop=index_bedrock+2,nsoil_PEM
213!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
214!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
215!      enddo
216!      enddo
217
218
219     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
220     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
221
222   
223   else
224! predictor corrector: restart from year 1 of the GCM and build the evolution of
225! tsoil at depth
226
227    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
228    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
229    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
230    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
231    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
232
233
234     do iloop = nsoil_GCM+1,nsoil_PEM
235       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
236     enddo
237   endif !found
238
239    do it = 1,timelen
240        do isoil = nsoil_GCM+1,nsoil_PEM
241        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
242        enddo
243     enddo
244      do isoil = nsoil_GCM+1,nsoil_PEM
245        do ig = 1,ngrid
246        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)
247        enddo
248      enddo
249
250   
251ENDDO ! islope
252
253print *,'PEMETAT0: SOIL TEMP  DONE'
254
255!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
256
257!3. Ice Table
258
259
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     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
265     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
266     do islope = 1,nslope
267     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
268     enddo
269   endif
270
271print *,'PEMETAT0: ICE TABLE  DONE'
272
273!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
274
275!4. CO2 & H2O Adsorption
276
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
307print *,'PEMETAT0: CO2 & H2O adsorption done  '
308    endif
309endif ! soil_pem
310
311
312!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
313
314!. 5 water reservoir
315
316#ifndef CPP_STD   
317   call get_field("water_reservoir",water_reservoir,found)
318   if(.not.found) then
319      write(*,*)'Pemetat0: failed loading <water_reservoir>'
320      write(*,*)'will reconstruct the values from watercaptag'
321      do ig=1,ngrid
322        if(watercaptag(ig)) then
323           water_reservoir(ig)=water_reservoir_nom
324        else
325           water_reservoir(ig)=0.
326        endif
327      enddo
328   endif
329#endif
330
331
332call close_startphy
333
334else !No startfi, let's build all by hand
335
336    year_PEM=0
337
338    if(soil_pem) then
339
340!a) Thermal inertia
341   do islope = 1,nslope
342      do ig = 1,ngrid
343
344          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
345!!! transition
346             delta = depth_breccia
347             TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
348            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
349                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
350
351          do iloop=index_breccia+2,index_bedrock
352            TI_PEM(ig,iloop,islope) = TI_breccia
353         enddo     
354         else ! we keep the high ti values
355           do iloop=index_breccia+1,index_bedrock
356                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
357           enddo
358         endif
359
360!! transition
361             delta = depth_bedrock
362             TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
363            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
364                      ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
365          do iloop=index_bedrock+2,nsoil_PEM
366            TI_PEM(ig,iloop,islope) = TI_bedrock
367         enddo   
368      enddo
369enddo
370
371 do iloop = 1,nsoil_GCM
372   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
373 enddo
374!!! zone de transition
375delta = depth_breccia
376do ig = 1,ngrid
377      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
378inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
379            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
380                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
381
382
383 do iloop = index_breccia+2,index_bedrock
384
385       inertiedat_PEM(ig,iloop) = TI_breccia
386   
387 enddo
388else
389   do iloop = index_breccia+1,index_bedrock
390       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
391    enddo
392
393endif
394enddo
395
396
397!!! zone de transition
398delta = depth_bedrock
399do ig = 1,ngrid
400inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
401            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
402                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
403enddo
404
405
406
407 do iloop = index_bedrock+2, nsoil_PEM
408   do ig = 1,ngrid
409      inertiedat_PEM(ig,iloop) = TI_bedrock
410   enddo
411 enddo
412
413print *,'PEMETAT0: TI  DONE'
414
415!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
416
417
418!b) Soil temperature
419    do islope = 1,nslope
420!     do ig = 1,ngrid
421!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
422!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
423!
424!       do iloop=index_breccia+2,index_bedrock
425!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
426!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
427!      enddo
428!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
429!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
430
431!       do iloop=index_bedrock+2,nsoil_PEM
432!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
433!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
434!      enddo
435     
436!       enddo
437
438      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
439      call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
440
441   
442    do it = 1,timelen
443        do isoil = nsoil_GCM+1,nsoil_PEM
444        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
445        enddo
446     enddo
447 
448      do isoil = nsoil_GCM+1,nsoil_PEM
449        do ig = 1,ngrid
450        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)
451        enddo
452      enddo
453enddo !islope
454print *,'PEMETAT0: TSOIL DONE  '
455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456!c) Ice table   
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   
463
464
465!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
466
467print *,'PEMETAT0: Ice table DONE  '
468
469
470
471!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
472!d) Regolith adsorbed
473
474 if(adsorption_pem) then
475    m_co2_regolith_phys(:,:,:) = 0.
476    m_h2o_regolith_phys(:,:,:) = 0.
477
478    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
479                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
480   
481    deltam_co2_regolith_phys(:) = 0.
482    deltam_h2o_regolith_phys(:) = 0. 
483 endif
484
485print *,'PEMETAT0: CO2 adsorption done  '
486
487endif !soil_pem
488
489
490
491!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492   
493
494
495!. e) water reservoir
496
497#ifndef CPP_STD   
498
499      do ig=1,ngrid
500        if(watercaptag(ig)) then
501           water_reservoir(ig)=water_reservoir_nom
502        else
503           water_reservoir(ig)=0.
504        endif
505      enddo
506
507#endif
508
509endif ! of if (startphy_file)
510
511if(soil_pem) then
512!! Sanity check
513    DO ig = 1,ngrid
514      DO islope = 1,nslope
515        DO iloop = 1,nsoil_PEM
516           if(isnan(tsoil_PEM(ig,iloop,islope))) then
517         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1)
518           endif
519        ENDDO
520      ENDDO
521     ENDDO
522endif!soil_pem
523
524
525
526   END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.