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

Last change on this file since 2889 was 2888, checked in by llange, 2 years ago

PEM

  • Following r-2886-5, debugging some issues ( conservation of H2O, water_reservoir not initialized, info_pem which was not working)
  • Cleaning of some redundant routines (e.g., stopping criterion of the PEM) and variables renamed (e.g., alpha_criterion now transofrmed in ice_criterion)
  • Ice table subroutine is now in a module and recalled "compute_ice_table_equilibrium" (v.s. dynamic ice table, efforts ongoing)
  • Abort_pem introduced to avoid just stopping the pem with raw "stop" in the codes
  • conf_pem has been improved so that values are not hard coded (e.g., geothermal flux, mass of water_reservoir)
  • Regolith thermal properties updated in a cleaner way

(Not atomic commit, sorry)
LL & RV
File size: 21.5 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
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
92
93startpem_file = .false.
94!1. Run
95
96if (startpem_file) then
97   ! open pem initial state file:
98   call open_startphy(filename)
99
100    if(soil_pem) then
101
102   call get_var("Time",year_PEM_read,found)
103   year_PEM=INT(year_PEM_read)
104   if(.not.found) then
105      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
106   else
107      write(*,*)'year_PEM of startpem=', year_PEM
108   endif
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,nsoil_GCM,islope).lt.TI_breccia) then
122!!! transition
123             delta = 50.
124             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
125            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
126                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
127            do iloop=nsoil_GCM+2,n_1km
128              TI_PEM(ig,iloop,islope) = TI_breccia
129            enddo     
130          else ! we keep the high ti values
131            do iloop=nsoil_GCM+1,n_1km
132              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
133            enddo
134          endif ! TI PEM and breccia comparison
135!! transition
136          delta = 1000.
137          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
138               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
139               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
140          do iloop=n_1km+2,nsoil_PEM
141            TI_PEM(ig,iloop,islope) = TI_bedrock
142          enddo   
143      enddo
144   else
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 = 50.
161do ig = 1,ngrid
162if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
163inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
164            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
165                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
166
167 do iloop = nsoil_GCM+2,n_1km
168       inertiedat_PEM(ig,iloop) = TI_breccia
169  enddo
170
171else
172   do iloop=nsoil_GCM+1,n_1km
173      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
174   enddo
175endif ! comparison ti breccia
176enddo!ig
177
178!!! zone de transition
179delta = 1000.
180do ig = 1,ngrid
181inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
182            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
183                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
184enddo
185
186 do iloop = n_1km+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,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
205        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
206
207       do iloop=nsoil_GCM+2,n_1km
208            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
209            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
210      enddo
211        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
212        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) 
213
214       do iloop=n_1km+2,nsoil_PEM
215            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
216            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
217      enddo
218          do iloop = nsoil_GCM+1,nsoil_PEM
219             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
220          enddo
221      enddo
222!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
223!     do iloop = 1,2000000
224!      write(*,*) 'iloop,islope',islope,iloop
225!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
226!     enddo
227   
228   else
229! predictor corrector: restart from year 1 of the GCM and build the evolution of
230! tsoil at depth
231
232    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
233    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
234    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
235    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
236    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
237
238
239     do iloop = nsoil_GCM+1,nsoil_PEM
240       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
241     enddo
242   endif
243
244    do it = 1,timelen
245        do isoil = nsoil_GCM+1,nsoil_PEM
246        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
247        enddo
248     enddo
249      do isoil = nsoil_GCM+1,nsoil_PEM
250        do ig = 1,ngrid
251        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
252        enddo
253      enddo
254
255   
256ENDDO ! islope
257
258print *,'PEMETAT0: SOIL TEMP  DONE'
259
260!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261
262!3. Ice Table
263
264
265
266     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
267     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
268
269
270print *,'PEMETAT0: ICE TABLE  DONE'
271
272!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273
274!4. CO2 & H2O Adsorption
275 if(adsorption_pem) then
276  DO islope=1,nslope
277   write(num,fmt='(i2.2)') islope
278   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
279    if((.not.found)) then
280       m_co2_regolith_phys(:,:,:) = 0.
281    endif
282    exit
283  ENDDO
284
285  DO islope=1,nslope
286   write(num,fmt='(i2.2)') islope
287   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
288    if((.not.found2)) then
289       m_h2o_regolith_phys(:,:,:) = 0.
290    endif
291    exit
292  ENDDO
293
294    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
295                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
296
297    if((.not.found)) then
298      deltam_co2_regolith_phys(:) = 0.
299    endif
300    if((.not.found2)) then
301       deltam_h2o_regolith_phys(:) = 0. 
302    endif
303print *,'PEMETAT0: CO2 & H2O adsorption done  '
304    endif
305endif ! soil_pem
306
307
308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
309
310!. 5 water reservoir
311
312#ifndef CPP_STD   
313   call get_field("water_reservoir",water_reservoir,found)
314   if(.not.found) then
315      write(*,*)'Pemetat0: failed loading <water_reservoir>'
316      write(*,*)'will reconstruct the values from watercaptag'
317      do ig=1,ngrid
318        if(watercaptag(ig)) then
319           water_reservoir=water_reservoir_nom
320        else
321           water_reservoir=0.
322        endif
323      enddo
324   endif
325#endif
326
327
328call close_startphy
329
330else !No startfi, let's build all by hand
331
332    year_PEM=0
333
334    if(soil_pem) then
335
336!a) Thermal inertia
337   do islope = 1,nslope
338      do ig = 1,ngrid
339
340          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
341!!! transition
342             delta = 50.
343             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
344            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
345                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
346
347          do iloop=nsoil_GCM+2,n_1km
348            TI_PEM(ig,iloop,islope) = TI_breccia
349         enddo     
350         else ! we keep the high ti values
351           do iloop=nsoil_GCM+1,n_1km
352                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
353           enddo
354         endif
355
356!! transition
357             delta = 1000.
358             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
359            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
360                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
361          do iloop=n_1km+2,nsoil_PEM
362            TI_PEM(ig,iloop,islope) = TI_bedrock
363         enddo   
364      enddo
365enddo
366
367 do iloop = 1,nsoil_GCM
368   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
369 enddo
370!!! zone de transition
371delta = 50.
372do ig = 1,ngrid
373      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
374inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
375            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
376                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
377
378
379 do iloop = nsoil_GCM+2,n_1km
380
381       inertiedat_PEM(ig,iloop) = TI_breccia
382   
383 enddo
384else
385   do iloop = nsoil_GCM+1,n_1km
386       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
387    enddo
388
389endif
390enddo
391
392
393!!! zone de transition
394delta = 1000.
395do ig = 1,ngrid
396inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
397            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
398                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
399enddo
400
401
402
403 do iloop = n_1km+2, nsoil_PEM
404   do ig = 1,ngrid
405      inertiedat_PEM(ig,iloop) = TI_bedrock
406   enddo
407 enddo
408
409print *,'PEMETAT0: TI  DONE'
410
411!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
412
413
414!b) Soil temperature
415    do islope = 1,nslope
416     do ig = 1,ngrid
417        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
418        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
419
420       do iloop=nsoil_GCM+2,n_1km
421            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
422            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
423      enddo
424        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
425        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
426
427       do iloop=n_1km+2,nsoil_PEM
428            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
429            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
430      enddo
431     
432           do iloop = nsoil_GCM+1,nsoil_PEM
433             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
434          enddo
435       enddo
436!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
437!     do iloop = 1,2000000
438!      write(*,*) 'islope,iloop',islope,iloop
439!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
440!     enddo
441 
442   
443    do it = 1,timelen
444        do isoil = nsoil_GCM+1,nsoil_PEM
445        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
446        enddo
447     enddo
448      do isoil = nsoil_GCM+1,nsoil_PEM
449        do ig = 1,ngrid
450        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
451        enddo
452      enddo
453enddo
454print *,'PEMETAT0: TSOIL DONE  '
455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456!c) Ice table   
457     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
458     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
459
460
461!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
462
463print *,'PEMETAT0: Ice table DONE  '
464
465
466
467!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
468!d) Regolith adsorbed
469
470 if(adsorption_pem) then
471    m_co2_regolith_phys(:,:,:) = 0.
472    m_h2o_regolith_phys(:,:,:) = 0.
473
474    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
475                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
476   
477    deltam_co2_regolith_phys(:) = 0.
478    deltam_h2o_regolith_phys(:) = 0. 
479 endif
480
481print *,'PEMETAT0: CO2 adsorption done  '
482
483endif !soil_pem
484
485endif ! of if (startphy_file)
486
487!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
488   
489
490
491!. e) water reservoir
492
493#ifndef CPP_STD   
494
495      do ig=1,ngrid
496        if(watercaptag(ig)) then
497           water_reservoir=water_reservoir_nom
498        else
499           water_reservoir=0.
500        endif
501      enddo
502
503#endif
504
505if(soil_pem) then
506!! Sanity check
507    DO ig = 1,ngrid
508      DO islope = 1,nslope
509        DO iloop = 1,nsoil_PEM
510           if(isnan(tsoil_PEM(ig,iloop,islope))) then
511         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1)
512           endif
513        ENDDO
514      ENDDO
515     ENDDO
516endif!soil_pem
517
518
519
520   END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.