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

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

Mars PEM

  • Following -r 2963, adapt the PEM so that the density of water wapor is correctly computed
  • I've also modified the way Tsoil is initialized, but it might be changed in the future

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