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

Last change on this file since 3015 was 2985, checked in by romain.vande, 2 years ago

Generic PEM :

Adapt PEM to run with the generic model.
(CPP_STD keyword to exclude some part of the code at the compilation)

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