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

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

PEM
Update the codes for subsurface evolution to run with XIOS
1) Water density at the surface and in the soil is now read in the XIOS file
2) Reshape routine created as XIOS output have one element less in the longitude grid
3) The ice table is now computed using these water densities
4) Minor fixs in the main, adsorption module, and tendencies evolutions

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