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

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

MARS PEM:

  • Add a PEMETAT0 that read "startfi_pem.nc"
  • Add the soil in the model: soil temperature, thermal properties, ice table
  • Add a routine that compute CO2 + H2O adsorption
  • Minor corrections in PEM.F90

LL
File size: 16.2 KB
Line 
1   SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
2                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_yr1_inst,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
3   
4   use iostart_PEM, only:  open_startphy, close_startphy, get_field
5   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
6   use comsoil_h, only:  volcapa,inertiedat
7   use ini_soil_mod, only:  ini_icetable
8   use soil_evolution_mod, only: soil_pem_CN
9   use adsorption_mod, only : regolith_co2adsorption
10  implicit none
11
12 
13  character(len=*), intent(in) :: filename
14  LOGICAL,intent(in) :: startpem_file
15  character*8 :: fichnom
16  integer,intent(in) :: ngrid
17  integer,intent(in) :: nsoil_GCM
18  integer,intent(in) :: nsoil_PEM
19  integer,intent(in) :: nslope
20  integer,intent(in) :: timelen
21  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)                          ! surface temperature at the first year of GCM call
22  real,intent(in) :: tsurf_ave(ngrid,nslope)                  ! surface temperature at the current year
23  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
24  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
25  real,intent(in) :: ps_yr1_inst(ngrid,timelen)
26  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
27  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
28  real,intent(in) :: timestep       ! time step
29  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)
30  real,intent(in) :: tend_co2glaciers(ngrid,nslope)
31  real,intent(in) :: co2ice(ngrid,nslope)
32  real,intent(in) :: waterice(ngrid,nslope)
33  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)
34
35! outputs
36
37
38
39
40
41  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
42  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) temperature (K)
43  real,intent(inout) :: ice_table(ngrid,nslope) ! (m)
44  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature (k)
45  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed (kg/m^2)
46  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption
47! local
48   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)
49   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)
50   LOGICAL :: found
51   integer :: iloop,ig,islope,it,isoil
52   REAL :: TI_breccia = 750.
53   REAL :: TI_bedrock = 2300.
54   real :: kcond
55   real :: delta
56   CHARACTER*2 :: num
57   real :: tsoil_saved(ngrid,nsoil_PEM)
58   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)
59   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
60   real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
61   real :: alph_tmp(nsoil_PEM-1)
62   real :: beta_tmp(nsoil_PEM-1)
63   real :: co2_ads_prev(ngrid)
64   real :: ps_ave_yr1(ngrid)
65   real :: ps_ave_yr2(ngrid)
66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
67!!!
68!!! Purpose: read start_pem. Need a specific iostart_PEM
69!!!
70!!! Order: 1. Thermal Inertia
71!!!        2. Soil Temperature
72!!!        3. Ice table
73!!!        4. Mass of CO2 adsorbed
74!!!
75!!! /!\ This order must be respected !
76!!! Author: LL
77!!!
78!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
79
80
81! 0. Some initializations
82 ps_ave_yr1(:) = sum(ps_yr1_inst(:,:),2)/timelen
83 ps_ave_yr2(:) = sum(ps_inst(:,:),2)/timelen
84
85!1. Run
86
87
88
89
90if (startpem_file) then
91   ! open pem initial state file:
92   call open_startphy(filename)
93 
94!1. Thermal Inertia
95! a. General case
96
97DO islope=1,nslope
98  write(num,fmt='(i2.2)') islope
99  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
100   if(.not.found) then
101      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
102      write(*,*)'will reconstruct the values of TI_PEM'
103
104      do ig = 1,ngrid
105
106          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
107!!! transition
108             delta = 50.
109             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
110            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
111                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
112
113          do iloop=nsoil_GCM+2,n_1km
114            TI_PEM(ig,iloop,islope) = TI_breccia
115         enddo     
116         else ! we keep the high ti values
117           do iloop=nsoil_GCM+1,n_1km
118                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
119           enddo
120
121         endif ! TI PEM and breccia comparison
122
123
124!! transition
125             delta = 1000.
126             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
127            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
128                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
129          do iloop=n_1km+2,nsoil_PEM
130            TI_PEM(ig,iloop,islope) = TI_bedrock
131         enddo   
132      enddo
133
134   else
135     do iloop = nsoil_GCM+1,nsoil_PEM
136       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.
137     enddo
138   endif ! not found
139ENDDO ! islope
140
141print *,'PEMETAT0: THERMAL INERTIA DONE'
142
143
144
145! b. Special case for inertiedat, inertiedat_PEM
146call get_field("inertiedat_PEM",inertiedat_PEM,found)
147if(.not.found) then
148 do iloop = 1,nsoil_GCM
149   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
150 enddo
151
152 
153!!! zone de transition
154delta = 50.
155do ig = 1,ngrid
156if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
157inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
158            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
159                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
160
161
162 do iloop = nsoil_GCM+2,n_1km
163       inertiedat_PEM(ig,iloop) = TI_breccia
164  enddo
165
166else
167   do iloop=nsoil_GCM+1,n_1km
168      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
169   enddo
170endif ! comparison ti breccia
171enddo!ig
172
173!!! zone de transition
174delta = 1000.
175do ig = 1,ngrid
176inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
177            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
178                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
179enddo
180
181
182
183 do iloop = n_1km+2, nsoil_PEM
184   do ig = 1,ngrid
185      inertiedat_PEM(ig,iloop) = TI_bedrock
186   enddo
187 enddo
188endif ! not found
189
190
191!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192
193!2. Soil Temperature
194
195DO islope=1,nslope
196  write(num,fmt='(i2.2)') islope
197   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
198  if(.not.found) then
199      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
200      write(*,*)'will reconstruct the values of Tsoil'
201      do ig = 1,ngrid
202        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
203        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
204
205       do iloop=nsoil_GCM+2,n_1km
206            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
207            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
208      enddo
209        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
210        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) 
211
212       do iloop=n_1km+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,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
215      enddo
216    enddo
217
218   else
219! predictor corrector: restart from year 1 of the GCM and build the evolution of
220! tsoil at depth
221
222    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
223    call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
224    call soil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
225    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
226    call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
227
228
229     do iloop = nsoil_GCM+1,nsoil_PEM
230       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
231     enddo
232   endif
233
234    do it = 1,timelen
235        do isoil = nsoil_GCM+1,nsoil_PEM
236        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
237        enddo
238     enddo
239
240
241
242   
243ENDDO
244print *,'PEMETAT0: SOIL TEMP  DONE'
245
246!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
247
248!3. Ice Table
249
250 
251   call get_field("ice_table",ice_table,found)
252   if(.not.found) then
253      write(*,*)'PEM settings: failed loading <Ice Table>'
254      write(*,*)'will reconstruct the values of ice table'
255
256
257      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
258
259
260   else
261   ! update ice table
262     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
263
264
265   endif
266
267print *,'PEMETAT0: ICE TABLE  DONE'
268
269
270
271!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272
273
274!4. CO2 Adsorption
275DO islope=1,nslope
276  write(num,fmt='(i2.2)') islope
277   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
278   if(.not.found) then
279      write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
280      write(*,*)'will reconstruct the values of co2 adsorbded'
281   m_co2_regolith_phys(:,:,:) = 0.
282   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
283
284   deltam_co2_regolith_phys(:) = 0.
285   exit
286   endif
287ENDDO
288
289if (found) then
290   DO iloop = 1,nsoil_GCM
291      tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
292
293   ENDDO 
294
295   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
296
297
298endif
299print *,'PEMETAT0: CO2 adsorption done  '
300
301
302
303call close_startphy
304
305
306else !No startfi, let's build all by hand
307
308
309
310!a) Thermal inertia
311   do islope = 1,nslope
312      do ig = 1,ngrid
313
314          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
315!!! transition
316             delta = 50.
317             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
318            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
319                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
320
321          do iloop=nsoil_GCM+2,n_1km
322            TI_PEM(ig,iloop,islope) = TI_breccia
323         enddo     
324         else ! we keep the high ti values
325           do iloop=nsoil_GCM+1,n_1km
326                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
327           enddo
328
329         endif
330
331
332!! transition
333             delta = 1000.
334             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
335            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
336                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
337          do iloop=n_1km+2,nsoil_PEM
338            TI_PEM(ig,iloop,islope) = TI_bedrock
339         enddo   
340      enddo
341
342enddo
343
344
345 do iloop = 1,nsoil_GCM
346   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
347 enddo
348!!! zone de transition
349delta = 50.
350do ig = 1,ngrid
351      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
352inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
353            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
354                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
355
356
357 do iloop = nsoil_GCM+2,n_1km
358
359       inertiedat_PEM(ig,iloop) = TI_breccia
360   
361 enddo
362else
363   do iloop = nsoil_GCM+1,n_1km
364       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
365    enddo
366
367endif
368enddo
369
370
371!!! zone de transition
372delta = 1000.
373do ig = 1,ngrid
374inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
375            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
376                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
377enddo
378
379
380
381 do iloop = n_1km+2, nsoil_PEM
382   do ig = 1,ngrid
383      inertiedat_PEM(ig,iloop) = TI_bedrock
384   enddo
385 enddo
386
387print *,'PEMETAT0: TI  DONE'
388
389!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
390
391
392!b) Soil temperature
393    do islope = 1,nslope
394
395     write(*,*) "islope=",islope
396     do ig = 1,ngrid
397        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
398        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
399
400       do iloop=nsoil_GCM+2,n_1km
401            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
402            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
403      enddo
404        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
405        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
406
407       do iloop=n_1km+2,nsoil_PEM
408            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
409            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
410      enddo
411!     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
412     enddo
413   
414    do it = 1,timelen
415        do isoil = nsoil_GCM+1,nsoil_PEM
416        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
417        enddo
418     enddo
419enddo
420print *,'PEMETAT0: TSOIL DONE  '
421!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
422!c) Ice table   
423      do islope = 1,nslope
424      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), & 
425       tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
426      enddo
427!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428
429print *,'PEMETAT0: Ice table DONE  '
430
431
432
433!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
434!d) Regolith adsorbed
435
436   m_co2_regolith_phys(:,:,:) = 0.
437   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
438   deltam_co2_regolith_phys(:) = 0.
439
440print *,'PEMETAT0: CO2 adsorption done  '
441
442
443
444
445
446
447
448
449
450
451endif ! of if (startphy_file)
452
453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
454    DO ig = 1,ngrid
455      DO islope = 1,nslope
456
457     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)
458
459      ENDDO
460    ENDDO
461
462
463
464!! small test
465    DO ig = 1,ngrid
466      DO islope = 1,nslope
467        DO iloop = 1,nsoil_PEM
468           if(isnan(tsoil_PEM(ig,iloop,islope))) then
469         write(*,*) "failed nan construction", ig, iloop, islope
470           stop
471           endif
472        ENDDO
473      ENDDO
474     ENDDO
475
476      write(*,*) "construction ok, no nan"
477     
478
479   END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.