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

Last change on this file since 2857 was 2855, checked in by llange, 3 years ago

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