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

Last change on this file since 3876 was 3850, checked in by jbclement, 6 months ago

PEM:

  • Handling the case where the new date is positive (present-day climate) in "recomp_orb_param_mod.F90", especially because of rounding errors.
  • Correction in "visu_evol_layering.py" to get the proper extent of the stratication + replacing 'imshow' by 'pcolormesh' to better process the data.
  • Addition of a graphic to plot the changes in obliquity directly onto the evolution of the stratification over time for better analysis. Obliquity curve is colored according to stratification gain/loss.
  • Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch.
  • Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837.

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