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

Last change on this file since 3040 was 3039, checked in by jbclement, 2 years ago

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