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

Last change on this file since 3983 was 3983, checked in by jbclement, 10 days ago

PEM:

  • Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
  • Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
  • Ice reservoirs representation in the PEM is modernized.

JBC
File size: 27.7 KB
Line 
1MODULE pemetat0_mod
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
10                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice,                         &
11                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,layerings_map)
12
13use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
14use comsoil_h_PEM,              only: soil_pem, layer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
15use comsoil_h,                  only: inertiedat
16use adsorption_mod,             only: regolith_adsorption, adsorption_pem
17use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
18use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
19use soil_thermalproperties_mod, only: update_soil_thermalproperties
20use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
21use abort_pem_mod,              only: abort_pem
22use compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
23use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo, add_stratum, ini_layering
24use callkeys_mod,               only: startphy_file
25use glaciers_mod,               only: rho_co2ice, rho_h2oice
26use comcstfi_h,                 only: r, mugaz, pi
27use surfdat_h,                  only: watercaptag, perennial_co2ice, qsurf
28use metamorphism,               only: frost4PCM, iPCM_h2ofrost, iPCM_co2frost
29
30implicit none
31
32character(*),                   intent(in) :: filename      ! name of the startfi_PEM.nc
33integer,                        intent(in) :: ngrid         ! # of physical grid points
34integer,                        intent(in) :: nsoil_PCM     ! # of vertical grid points in the PCM
35integer,                        intent(in) :: nsoil_PEM     ! # of vertical grid points in the PEM
36integer,                        intent(in) :: nslope        ! # of sub-grid slopes
37integer,                        intent(in) :: timelen       ! # time samples
38real,                           intent(in) :: timestep      ! time step [s]
39real,                           intent(in) :: ps_avg_global ! mean average pressure on the planet [Pa]
40real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
41real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
42real, dimension(ngrid,timelen), intent(in) :: q_co2         ! MMR tracer co2 [kg/kg]
43real, dimension(ngrid,timelen), intent(in) :: q_h2o         ! MMR tracer h2o [kg/kg]
44real, dimension(ngrid,timelen), intent(in) :: ps_timeseries ! surface pressure [Pa]
45real, dimension(ngrid,nslope),  intent(in) :: d_h2oice      ! tendencies on h2o ice
46real, dimension(ngrid,nslope),  intent(in) :: d_co2ice      ! tendencies on co2 ice
47real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
48! outputs
49real, dimension(ngrid),                  intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
50real, dimension(ngrid),                  intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
51real, dimension(ngrid,nslope),           intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
52real, dimension(ngrid,nslope),           intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
53type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map       ! Layerings
54real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
55real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
56real, dimension(ngrid,nslope),           intent(inout) :: icetable_depth      ! Ice table depth [m]
57real, dimension(ngrid,nslope),           intent(inout) :: icetable_thickness  ! Ice table thickness [m]
58real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
59real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
60real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
61real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
62! local
63real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem    ! soil temperature saved in the start [K]
64real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM       ! soil thermal inertia saved in the start [SI]
65logical                                 :: found, found2     ! check if variables are found in the start
66integer                                 :: iloop, ig, islope ! index for loops
67real                                    :: delta             ! Depth of the interface regolith-breccia, breccia -bedrock [m]
68character(2)                            :: num               ! intermediate string to read PEM start sloped variables
69logical                                 :: startpem_file     ! boolean to check if we read the startfile or not
70real, dimension(:,:,:,:), allocatable   :: stratif_array     ! Array for layerings
71
72!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73!!!
74!!! Purpose: read start_pem. Need a specific iostart_PEM
75!!!
76!!! Order: 0. Previous year of the PEM run
77!!!           Ice initialization
78!!!        1. Thermal Inertia
79!!!        2. Soil Temperature
80!!!        3. Ice table
81!!!        4. Mass of CO2 & H2O adsorbed
82!!!
83!!! /!\ This order must be respected!
84!!! Author: LL
85!!!
86!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87
88!0.1 Check if the start_PEM exist.
89
90inquire(file = filename,exist = startpem_file)
91
92write(*,*)'Is there a "startpem" file?',startpem_file
93
94!0.2 Set to default values
95icetable_depth = -1. ! by default, no ice table
96icetable_thickness = -1.
97ice_porefilling = 0.
98!1. Run
99if (startpem_file) then
100    write(*,*)'There is a "startpem.nc"!'
101    ! open pem initial state file:
102    call open_startphy(filename)
103
104!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105    ! H2O ice
106    h2o_ice(:,:) = 0.
107    call get_field("h2o_ice",h2o_ice,found)
108    if (.not. found) then
109        write(*,*)'Pemetat0: failed loading <h2o_ice>'
110        write(*,*)'will reconstruct the values from watercaptag'
111        write(*,*)'with default value ''ini_huge_h2oice'''
112        do ig = 1,ngrid
113            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice + qsurf(ig,iPCM_h2ofrost,:) - frost4PCM(ig,:)%h2o
114        enddo
115    endif
116
117    ! CO2 ice
118    co2_ice(:,:) = 0.
119    call get_field("co2_ice",co2_ice,found)
120    if (.not. found) then
121        write(*,*)'Pemetat0: failed loading <co2_ice>'
122        write(*,*)'will reconstruct the values from perennial_co2ice'
123        co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - frost4PCM(:,:)%co2
124    endif
125
126!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
127    ! Layerings
128    nb_str_max = 68
129    if (layering_algo) then
130        found = inquire_dimension("nb_str_max")
131        if (.not. found) then
132            write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
133            write(*,*) '''nb_str_max'' is set to 68 (sub-surface layers)!'
134        else
135            nb_str_max = int(inquire_dimension_length('nb_str_max'))
136        endif
137        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
138        stratif_array = 0.
139        do islope = 1,nslope
140            write(num,'(i2.2)') islope
141            call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,1),found)
142            if (.not. found) then
143                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_top_elevation>!'
144                exit
145            endif
146            call get_field('stratif_slope'//num//'_h_co2ice',stratif_array(:,islope,:,2),found)
147            if (.not. found) then
148                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_co2ice>!'
149                exit
150            endif
151            call get_field('stratif_slope'//num//'_h_h2oice',stratif_array(:,islope,:,3),found)
152            if (.not. found) then
153                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_h2oice>!'
154                exit
155            endif
156            call get_field('stratif_slope'//num//'_h_dust',stratif_array(:,islope,:,4),found)
157            if (.not. found) then
158                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_dust>!'
159                exit
160            endif
161            call get_field('stratif_slope'//num//'_h_pore',stratif_array(:,islope,:,5),found)
162            if (.not. found) then
163                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_pore>!'
164                exit
165            endif
166            call get_field('stratif_slope'//num//'_poreice_volfrac',stratif_array(:,islope,:,6),found)
167            if (.not. found) then
168                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_poreice_volfrac>!'
169                exit
170            endif
171        enddo ! islope
172        if (.not. found) then
173            write(*,*) 'So the layerings are initialized with sub-surface strata.'
174            write(*,*) 'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
175            do ig = 1,ngrid
176                if (watercaptag(ig)) then
177                    do islope = 1,nslope
178                        call ini_layering(layerings_map(ig,islope))
179                        call add_stratum(layerings_map(ig,islope),1.05*ini_huge_h2oice/rho_h2oice,0.,ini_huge_h2oice/rho_h2oice,0.05*ini_huge_h2oice/rho_h2oice,0.,0.)
180                    enddo
181                else
182                    do islope = 1,nslope
183                        call ini_layering(layerings_map(ig,islope))
184                        if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05*perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05*perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
185                    enddo
186                endif
187            enddo
188        else
189            call array2stratif(stratif_array,ngrid,nslope,layerings_map)
190        endif
191        deallocate(stratif_array)
192    endif
193
194    if (soil_pem) then
195
196!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197!1. Thermal Inertia
198! a. General case
199        do islope = 1,nslope
200            write(num,'(i2.2)') islope
201            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
202            if (.not. found) then
203                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
204                write(*,*)'will reconstruct the values of TI_PEM'
205
206                do ig = 1,ngrid
207                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
208                        !!! transition
209                        delta = depth_breccia
210                        TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
211                                                                   (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
212                                                                   ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
213                        do iloop = index_breccia + 2,index_bedrock
214                            TI_PEM(ig,iloop,islope) = TI_breccia
215                        enddo
216                    else ! we keep the high TI values
217                        do iloop = index_breccia + 1,index_bedrock
218                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
219                        enddo
220                    endif ! TI PEM and breccia comparison
221                    !! transition
222                    delta = depth_bedrock
223                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
224                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
225                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
226                    do iloop = index_bedrock + 2,nsoil_PEM
227                        TI_PEM(ig,iloop,islope) = TI_bedrock
228                    enddo
229                enddo
230            else ! found
231                do iloop = nsoil_PCM + 1,nsoil_PEM
232                    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.
233                enddo
234            endif ! not found
235        enddo ! islope
236
237        write(*,*) 'PEMETAT0: TI done'
238
239! b. Special case for inertiedat, inertiedat_PEM
240        call get_field("inertiedat_PEM",inertiedat_PEM,found)
241        if (.not. found) then
242            do iloop = 1,nsoil_PCM
243                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
244            enddo
245        !!! zone de transition
246            delta = depth_breccia
247            do ig = 1,ngrid
248                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
249                    inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
250                                                                (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
251                                                                ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
252
253                    do iloop = index_breccia + 2,index_bedrock
254                        inertiedat_PEM(ig,iloop) = TI_breccia
255                    enddo
256                else
257                    do iloop = index_breccia + 1,index_bedrock
258                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
259                    enddo
260                endif ! comparison ti breccia
261            enddo ! ig
262
263            !!! zone de transition
264            delta = depth_bedrock
265            do ig = 1,ngrid
266                inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
267                                                            (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
268                                                            ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
269            enddo
270
271            do iloop = index_bedrock + 2,nsoil_PEM
272                do ig = 1,ngrid
273                    inertiedat_PEM(ig,iloop) = TI_bedrock
274                enddo
275            enddo
276        endif ! not found
277
278!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
279!2. Soil Temperature
280        do islope = 1,nslope
281            write(num,fmt='(i2.2)') islope
282            call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
283            if (.not. found) then
284                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
285                write(*,*)'will reconstruct the values of Tsoil'
286                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
287                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
288            else
289! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
290!                      for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
291                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
292                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
293                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
294
295                tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) = tsoil_startpem(:,nsoil_PCM + 1:nsoil_PEM,islope)
296            endif !found
297
298            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
299        enddo ! islope
300        write(*,*) 'PEMETAT0: TSOIL done'
301
302!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303!3. Ice Table
304        if (icetable_equilibrium) then
305            call get_field("icetable_depth",icetable_depth,found)
306            if (.not. found) then
307                write(*,*)'PEM settings: failed loading <icetable_depth>'
308                write(*,*)'will reconstruct the values of the ice table given the current state'
309                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
310                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
311                do islope = 1,nslope
312                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
313                enddo
314            endif
315            write(*,*) 'PEMETAT0: ICE TABLE done'
316        else if (icetable_dynamic) then
317            call get_field("ice_porefilling",ice_porefilling,found)
318            if (.not. found) write(*,*)'PEM settings: failed loading <ice_porefilling>'
319            call get_field("icetable_depth",icetable_depth,found)
320            if (.not. found) then
321                write(*,*)'PEM settings: failed loading <icetable_depth>'
322                write(*,*)'will reconstruct the values of the ice table given the current state'
323                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
324                do islope = 1,nslope
325                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
326                enddo
327            endif
328            write(*,*) 'PEMETAT0: ICE TABLE done'
329        endif
330
331!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
332!4. CO2 & H2O Adsorption
333        if (adsorption_pem) then
334            do islope = 1,nslope
335                write(num,fmt='(i2.2)') islope
336                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
337                if (.not. found) then
338                    m_co2_regolith_phys = 0.
339                    exit
340                endif
341            enddo
342            do islope = 1,nslope
343                write(num,fmt='(i2.2)') islope
344                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
345                if (.not. found2) then
346                    m_h2o_regolith_phys = 0.
347                    exit
348                endif
349            enddo
350
351            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
352                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
353
354            if (.not. found) deltam_co2_regolith_phys = 0.
355            if (.not. found2) deltam_h2o_regolith_phys = 0.
356
357            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
358        endif ! adsorption_pem
359    endif ! soil_pem
360
361    call close_startphy
362
363else ! No startpem, let's build all by hand
364
365    write(*,*)'There is no "startpem.nc"!'
366
367    ! h2o ice
368    h2o_ice(:,:) = 0.
369    write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
370    do ig = 1,ngrid
371        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice + qsurf(ig,iPCM_h2ofrost,:) - frost4PCM(ig,:)%h2o
372    enddo
373
374    ! co2 ice
375    write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice'' found in the PCM.'
376    co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - frost4PCM(:,:)%co2
377
378    ! Layerings
379    nb_str_max = 68
380    if (layering_algo) then
381        write(*,*)'So the layerings are initialized with sub-surface strata.'
382        write(*,*)'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
383        do ig = 1,ngrid
384            if (watercaptag(ig)) then
385                do islope = 1,nslope
386                    call ini_layering(layerings_map(ig,islope))
387                    call add_stratum(layerings_map(ig,islope),1.05*ini_huge_h2oice/rho_h2oice,0.,ini_huge_h2oice/rho_h2oice,0.05*ini_huge_h2oice/rho_h2oice,0.,0.)
388                enddo
389            else
390                do islope = 1,nslope
391                    call ini_layering(layerings_map(ig,islope))
392                    if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05*perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05*perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
393                enddo
394            endif
395        enddo
396    endif
397
398    if (soil_pem) then
399
400!a) Thermal inertia
401        do islope = 1,nslope
402            do ig = 1,ngrid
403                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
404                    !!! transition
405                    delta = depth_breccia
406                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                   &
407                                                               (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
408                                                               ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
409                    do iloop = index_breccia + 2,index_bedrock
410                        TI_PEM(ig,iloop,islope) = TI_breccia
411                    enddo
412                else ! we keep the high ti values
413                    do iloop = index_breccia + 1,index_bedrock
414                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
415                    enddo
416                endif
417                !! transition
418                delta = depth_bedrock
419                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                   &
420                                                           (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
421                                                           ((layer_PEM(index_bedrock + 1) - delta)/(TI_breccia**2))))
422                do iloop = index_bedrock + 2,nsoil_PEM
423                    TI_PEM(ig,iloop,islope) = TI_bedrock
424                enddo
425            enddo
426        enddo
427
428        do iloop = 1,nsoil_PCM
429            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
430        enddo
431        !!! zone de transition
432        delta = depth_breccia
433        do ig = 1,ngrid
434            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
435                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
436                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
437                                                            ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
438                do iloop = index_breccia + 2,index_bedrock
439                    inertiedat_PEM(ig,iloop) = TI_breccia
440                enddo
441            else
442                do iloop = index_breccia + 1,index_bedrock
443                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
444                enddo
445            endif
446        enddo
447
448        !!! zone de transition
449        delta = depth_bedrock
450        do ig = 1,ngrid
451            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
452                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
453                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
454        enddo
455
456        do iloop = index_bedrock + 2,nsoil_PEM
457            do ig = 1,ngrid
458                inertiedat_PEM(ig,iloop) = TI_bedrock
459            enddo
460        enddo
461
462        write(*,*) 'PEMETAT0: TI done'
463
464!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
465!b) Soil temperature
466        do islope = 1,nslope
467            call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
468            call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
469
470! First raw initialization
471            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
472        enddo !islope
473        write(*,*) 'PEMETAT0: TSOIL done'
474
475!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
476!c) Ice table
477        if (icetable_equilibrium) then
478            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
479            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
480            do islope = 1,nslope
481                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
482            enddo
483            write(*,*) 'PEMETAT0: ICE TABLE done'
484        else if (icetable_dynamic) then
485            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
486            do islope = 1,nslope
487                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
488            enddo
489            write(*,*) 'PEMETAT0: ICE TABLE done'
490        endif
491
492!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493!d) Regolith adsorbed
494        if (adsorption_pem) then
495            m_co2_regolith_phys = 0.
496            m_h2o_regolith_phys = 0.
497            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
498                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
499            deltam_co2_regolith_phys = 0.
500            deltam_h2o_regolith_phys = 0.
501        endif
502
503        write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
504    endif !soil_pem
505endif ! of if (startphy_file)
506
507if (soil_pem) then ! Sanity check
508    do ig = 1,ngrid
509        do islope = 1,nslope
510            do iloop = 1,nsoil_PEM
511                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
512            enddo
513        enddo
514    enddo
515endif !soil_pem
516
517END SUBROUTINE pemetat0
518
519END MODULE pemetat0_mod
Note: See TracBrowser for help on using the repository browser.