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

Last change on this file since 2841 was 2835, checked in by romain.vande, 3 years ago

Mars PEM:
Introduction of the possibility to follow an orbital forcing.
Introduction of new control parameters.
Cleaning of the PEM (removing unused files, add comments and new files)

A file named run_PEM.def can be added to the run.def. It contains the following variables:

_ evol_orbit_pem: Boolean. Do you want to follow an orbital forcing predefined (read in ob_ex_lsp.asc for example)? (default=false)
_ year_bp_ini: Integer. Number of year before present to start the pem run if evol_orbit_pem=.true. , default=0
_ Max_iter_pem: Integer. Maximal number of iteration if none of the stopping criterion is reached and if evol_orbit_pem=.false., default=99999999
_ dt_pem: Integer. Time step of the PEM in year, default=1
_ alpha_criterion: Real. Acceptance rate of sublimating ice surface change, default=0.2
_ soil_pem: Boolean. Do you want to run with subsurface physical processes in the PEM? default=.true.

RV
File size: 16.4 KB
Line 
1   SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
2                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
3   
4   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
5   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
6   use comsoil_h, only:  volcapa,inertiedat
7   use ini_soil_mod, only:  ini_icetable
8   use soil_evolution_mod, only: soil_pem,soil_pem_CN
9   use adsorption_mod, only : regolith_co2adsorption
10   USE temps_mod_evol, ONLY: year_PEM
11
12  implicit none
13 
14  character(len=*), intent(in) :: filename
15  LOGICAL,intent(in) :: startpem_file
16  character*8 :: fichnom
17  integer,intent(in) :: ngrid
18  integer,intent(in) :: nsoil_GCM
19  integer,intent(in) :: nsoil_PEM
20  integer,intent(in) :: nslope
21  integer,intent(in) :: timelen
22  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)                          ! surface temperature at the first year of GCM call
23  real,intent(in) :: tsurf_ave(ngrid,nslope)                  ! surface temperature at the current year
24  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
25  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
26  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
27  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
28  real,intent(in) :: timestep       ! time step
29  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)
30  real,intent(in) :: tend_co2glaciers(ngrid,nslope)
31  real,intent(in) :: co2ice(ngrid,nslope)
32  real,intent(in) :: waterice(ngrid,nslope)
33  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)
34
35! outputs
36
37  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
38  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) temperature (K)
39  real,intent(inout) :: ice_table(ngrid,nslope) ! (m)
40  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature (k)
41  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed (kg/m^2)
42  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption
43! local
44   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)
45   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)
46   LOGICAL :: found
47   integer :: iloop,ig,islope,it,isoil
48   REAL :: TI_breccia = 750.
49   REAL :: TI_bedrock = 2300.
50   real :: kcond
51   real :: delta
52   CHARACTER*2 :: num
53   real :: tsoil_saved(ngrid,nsoil_PEM)
54   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)
55   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
56   real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
57   real :: alph_tmp(ngrid,nsoil_PEM-1)
58   real :: beta_tmp(ngrid,nsoil_PEM-1)
59   real :: co2_ads_prev(ngrid)
60   real :: year_PEM_read
61
62!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
63!!!
64!!! Purpose: read start_pem. Need a specific iostart_PEM
65!!!
66!!! Order: 1. Thermal Inertia
67!!!        2. Soil Temperature
68!!!        3. Ice table
69!!!        4. Mass of CO2 adsorbed
70!!!
71!!! /!\ This order must be respected !
72!!! Author: LL
73!!!
74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75
76!1. Run
77
78if (startpem_file) then
79   ! open pem initial state file:
80   call open_startphy(filename)
81
82    if(soil_pem) then
83
84   call get_var("Time",year_PEM_read,found)
85   year_PEM=INT(year_PEM_read)
86   if(.not.found) then
87      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
88   else
89      write(*,*)'year_PEM of startpem=', year_PEM
90   endif
91 
92!1. Thermal Inertia
93! a. General case
94
95DO islope=1,nslope
96  write(num,fmt='(i2.2)') islope
97  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
98   if(.not.found) then
99      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
100      write(*,*)'will reconstruct the values of TI_PEM'
101
102      do ig = 1,ngrid
103          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
104!!! transition
105             delta = 50.
106             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
107            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
108                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
109            do iloop=nsoil_GCM+2,n_1km
110              TI_PEM(ig,iloop,islope) = TI_breccia
111            enddo     
112          else ! we keep the high ti values
113            do iloop=nsoil_GCM+1,n_1km
114              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
115            enddo
116          endif ! TI PEM and breccia comparison
117!! transition
118          delta = 1000.
119          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
120               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
121               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
122          do iloop=n_1km+2,nsoil_PEM
123            TI_PEM(ig,iloop,islope) = TI_bedrock
124          enddo   
125      enddo
126   else
127     do iloop = nsoil_GCM+1,nsoil_PEM
128       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.
129     enddo
130   endif ! not found
131ENDDO ! islope
132
133print *,'PEMETAT0: THERMAL INERTIA DONE'
134
135! b. Special case for inertiedat, inertiedat_PEM
136call get_field("inertiedat_PEM",inertiedat_PEM,found)
137if(.not.found) then
138 do iloop = 1,nsoil_GCM
139   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
140 enddo
141!!! zone de transition
142delta = 50.
143do ig = 1,ngrid
144if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
145inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
146            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
147                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
148
149 do iloop = nsoil_GCM+2,n_1km
150       inertiedat_PEM(ig,iloop) = TI_breccia
151  enddo
152
153else
154   do iloop=nsoil_GCM+1,n_1km
155      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
156   enddo
157endif ! comparison ti breccia
158enddo!ig
159
160!!! zone de transition
161delta = 1000.
162do ig = 1,ngrid
163inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
164            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
165                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
166enddo
167
168 do iloop = n_1km+2, nsoil_PEM
169   do ig = 1,ngrid
170      inertiedat_PEM(ig,iloop) = TI_bedrock
171   enddo
172 enddo
173endif ! not found
174
175!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176
177!2. Soil Temperature
178
179DO islope=1,nslope
180  write(num,fmt='(i2.2)') islope
181   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
182  if(.not.found) then
183      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
184      write(*,*)'will reconstruct the values of Tsoil'
185      do ig = 1,ngrid
186        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
187        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
188
189       do iloop=nsoil_GCM+2,n_1km
190            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
191            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
192      enddo
193        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
194        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1)) 
195
196       do iloop=n_1km+2,nsoil_PEM
197            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
198            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
199      enddo
200    enddo
201
202   else
203! predictor corrector: restart from year 1 of the GCM and build the evolution of
204! tsoil at depth
205
206    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
207    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
208    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
209    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
210    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
211
212
213     do iloop = nsoil_GCM+1,nsoil_PEM
214       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
215     enddo
216   endif
217
218    do it = 1,timelen
219        do isoil = nsoil_GCM+1,nsoil_PEM
220        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
221        enddo
222     enddo
223   
224ENDDO
225
226print *,'PEMETAT0: SOIL TEMP  DONE'
227
228!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229
230!3. Ice Table
231
232   call get_field("ice_table",ice_table,found)
233   if(.not.found) then
234      write(*,*)'PEM settings: failed loading <Ice Table>'
235      write(*,*)'will reconstruct the values of ice table'
236
237      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
238
239   else
240   ! update ice table
241     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
242
243   endif
244
245print *,'PEMETAT0: ICE TABLE  DONE'
246
247!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
248
249!4. CO2 Adsorption
250DO islope=1,nslope
251  write(num,fmt='(i2.2)') islope
252   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
253   if(.not.found) then
254      write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
255      write(*,*)'will reconstruct the values of co2 adsorbded'
256   m_co2_regolith_phys(:,:,:) = 0.
257   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
258
259   deltam_co2_regolith_phys(:) = 0.
260   exit
261   endif
262ENDDO
263
264if (found) then
265   DO iloop = 1,nsoil_GCM
266      tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
267
268   ENDDO 
269
270   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
271
272
273endif
274print *,'PEMETAT0: CO2 adsorption done  '
275
276endif ! soil_pem
277
278call close_startphy
279
280else !No startfi, let's build all by hand
281
282    year_PEM=0
283
284    if(soil_pem) then
285
286!a) Thermal inertia
287   do islope = 1,nslope
288      do ig = 1,ngrid
289
290          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then
291!!! transition
292             delta = 50.
293             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
294            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
295                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
296
297          do iloop=nsoil_GCM+2,n_1km
298            TI_PEM(ig,iloop,islope) = TI_breccia
299         enddo     
300         else ! we keep the high ti values
301           do iloop=nsoil_GCM+1,n_1km
302                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
303           enddo
304         endif
305
306!! transition
307             delta = 1000.
308             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
309            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
310                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
311          do iloop=n_1km+2,nsoil_PEM
312            TI_PEM(ig,iloop,islope) = TI_bedrock
313         enddo   
314      enddo
315enddo
316
317 do iloop = 1,nsoil_GCM
318   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
319 enddo
320!!! zone de transition
321delta = 50.
322do ig = 1,ngrid
323      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then
324inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
325            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
326                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
327
328
329 do iloop = nsoil_GCM+2,n_1km
330
331       inertiedat_PEM(ig,iloop) = TI_breccia
332   
333 enddo
334else
335   do iloop = nsoil_GCM+1,n_1km
336       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
337    enddo
338
339endif
340enddo
341
342
343!!! zone de transition
344delta = 1000.
345do ig = 1,ngrid
346inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
347            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
348                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
349enddo
350
351
352
353 do iloop = n_1km+2, nsoil_PEM
354   do ig = 1,ngrid
355      inertiedat_PEM(ig,iloop) = TI_bedrock
356   enddo
357 enddo
358
359print *,'PEMETAT0: TI  DONE'
360
361!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
362
363
364!b) Soil temperature
365    do islope = 1,nslope
366
367     write(*,*) "islope=",islope
368     do ig = 1,ngrid
369        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
370        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
371
372       do iloop=nsoil_GCM+2,n_1km
373            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
374            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
375      enddo
376        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
377        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
378
379       do iloop=n_1km+2,nsoil_PEM
380            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
381            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
382      enddo
383!     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
384     enddo
385   
386    do it = 1,timelen
387        do isoil = nsoil_GCM+1,nsoil_PEM
388        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
389        enddo
390     enddo
391enddo
392print *,'PEMETAT0: TSOIL DONE  '
393!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
394!c) Ice table   
395      do islope = 1,nslope
396      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), & 
397       tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
398      enddo
399!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400
401print *,'PEMETAT0: Ice table DONE  '
402
403
404
405!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
406!d) Regolith adsorbed
407
408   m_co2_regolith_phys(:,:,:) = 0.
409   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
410   deltam_co2_regolith_phys(:) = 0.
411
412print *,'PEMETAT0: CO2 adsorption done  '
413
414endif !soil_pem
415
416endif ! of if (startphy_file)
417
418!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
419   
420     if(soil_pem) then
421
422    DO ig = 1,ngrid
423      DO islope = 1,nslope
424
425     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)
426
427      ENDDO
428    ENDDO
429
430!! small test
431    DO ig = 1,ngrid
432      DO islope = 1,nslope
433        DO iloop = 1,nsoil_PEM
434           if(isnan(tsoil_PEM(ig,iloop,islope))) then
435         write(*,*) "failed nan construction", ig, iloop, islope
436           stop
437           endif
438        ENDDO
439      ENDDO
440     ENDDO
441
442     endif!soil_pem
443
444      write(*,*) "construction ok, no nan"
445     
446
447   END SUBROUTINE
Note: See TracBrowser for help on using the repository browser.