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

Last change on this file since 3143 was 3143, checked in by jbclement, 12 months ago

PEM:
'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code.

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