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

Last change on this file since 2864 was 2863, checked in by llange, 3 years ago

PEM
H2O adsorption is now computed. Total mass of H2o adsorbded is written in the restart_PEM.
+Minor fixes and edit
LL & RV

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