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

Last change on this file since 2937 was 2937, checked in by llange, 20 months ago

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