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

Last change on this file since 2852 was 2849, checked in by llange, 3 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
RevLine 
[2794]1   SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
[2849]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)
[2794]4   
[2835]5   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
[2794]6   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
7   use comsoil_h, only:  volcapa,inertiedat
[2835]8   use soil_evolution_mod, only: soil_pem,soil_pem_CN
[2794]9   use adsorption_mod, only : regolith_co2adsorption
[2835]10   USE temps_mod_evol, ONLY: year_PEM
11
[2794]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)
[2849]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
[2794]37! outputs
[2779]38
[2794]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)
[2835]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
[2849]63   real :: alpha_clap_h2o = -6143.7
64   real :: beta_clap_h2o = 28.9074
[2835]65
[2794]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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2779]79
[2794]80!1. Run
[2779]81
[2794]82if (startpem_file) then
83   ! open pem initial state file:
84   call open_startphy(filename)
[2835]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
[2794]95 
96!1. Thermal Inertia
97! a. General case
[2779]98
[2794]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'
[2779]105
[2794]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))))
[2835]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
[2794]121!! transition
[2835]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))))
[2794]126          do iloop=n_1km+2,nsoil_PEM
127            TI_PEM(ig,iloop,islope) = TI_bedrock
[2835]128          enddo   
[2794]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
[2779]136
[2794]137print *,'PEMETAT0: THERMAL INERTIA DONE'
[2779]138
[2794]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)
[2835]144 enddo
[2794]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))))
[2779]152
[2794]153 do iloop = nsoil_GCM+2,n_1km
154       inertiedat_PEM(ig,iloop) = TI_breccia
155  enddo
[2779]156
[2794]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
[2779]163
[2794]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
[2779]171
[2794]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
[2779]183DO islope=1,nslope
184  write(num,fmt='(i2.2)') islope
[2794]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)
[2835]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)
[2794]213    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
[2835]214    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
[2794]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
[2849]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
[2794]233   
[2779]234ENDDO
[2835]235
[2794]236print *,'PEMETAT0: SOIL TEMP  DONE'
[2779]237
[2794]238!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
239
240!3. Ice Table
241
242
243
[2849]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)
[2794]246
247
248print *,'PEMETAT0: ICE TABLE  DONE'
249
250!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
251
252!4. CO2 Adsorption
[2779]253DO islope=1,nslope
254  write(num,fmt='(i2.2)') islope
[2794]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
[2849]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
[2794]267   endif
[2779]268ENDDO
269
270
[2794]271print *,'PEMETAT0: CO2 adsorption done  '
272
[2835]273endif ! soil_pem
[2794]274
275call close_startphy
276
277else !No startfi, let's build all by hand
278
[2835]279    year_PEM=0
[2794]280
[2835]281    if(soil_pem) then
[2794]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
[2779]301         endif
302
[2794]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
[2779]313
[2794]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))))
[2779]324
325
[2794]326 do iloop = nsoil_GCM+2,n_1km
[2779]327
[2794]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
[2779]335
[2794]336endif
337enddo
[2779]338
339
[2794]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
[2779]347
348
349
[2794]350 do iloop = n_1km+2, nsoil_PEM
351   do ig = 1,ngrid
352      inertiedat_PEM(ig,iloop) = TI_bedrock
353   enddo
354 enddo
[2779]355
[2794]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
[2849]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
[2794]390enddo
391print *,'PEMETAT0: TSOIL DONE  '
392!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393!c) Ice table   
[2849]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
[2794]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
[2835]413endif !soil_pem
[2794]414
415endif ! of if (startphy_file)
416
417!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2835]418   
419     if(soil_pem) then
420
[2794]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
[2835]441     endif!soil_pem
442
[2794]443      write(*,*) "construction ok, no nan"
444     
445
446   END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.