source: trunk/LMDZ.COMMON/libf/evolution/clim_state_init.F90 @ 4174

Last change on this file since 4174 was 4170, checked in by jbclement, 7 days ago

PEM:

  • Deletion of 'flux_ssice' ('zqdsdif_tot') from the "startfi.nc" as it is not needed.
  • Using the yearly average flux for the sublimating subsurface ice instead of the value at last timestep.
  • Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
  • Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).

JBC
File size: 22.3 KB
Line 
1MODULE clim_state_init
2!-----------------------------------------------------------------------
3! NAME
4!     clim_state_init
5!
6! DESCRIPTION
7!     Read the start files to initialize the climate state.
8!
9! AUTHORS & DATE
10!     JB Clement, 12/2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use numerics,  only: dp, di, k4
19use io_netcdf, only: open_nc, close_nc, get_var_nc, get_dim_nc, put_var_nc, start_name, start1D_name, startfi_name, startevo_name
20
21! DECLARATION
22! -----------
23implicit none
24
25contains
26!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27
28!=======================================================================
29SUBROUTINE read_start()
30!-----------------------------------------------------------------------
31! NAME
32!     read_start
33!
34! DESCRIPTION
35!     Read the file "start.nc".
36!
37! AUTHORS & DATE
38!     JB Clement, 12/2025
39!
40! NOTES
41!
42!-----------------------------------------------------------------------
43
44! DEPENDENCIES
45! ------------
46use geometry,   only: nlayer, nlon, nlat, ngrid
47use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM
48use tracers,    only: nq, qnames, set_q_PCM
49use stoppage,   only: stop_clean
50use display,    only: print_msg, LVL_NFO
51
52! DECLARATION
53! -----------
54implicit none
55
56! LOCAL VARIABLES
57! ---------------
58real(dp), dimension(nlayer + 1)           :: tmp1d
59real(dp), dimension(nlon + 1,nlat)        :: tmp2d
60real(dp), dimension(nlon + 1,nlat,nlayer) :: tmp3d
61logical(k4)                               :: here
62integer(di)                               :: i
63
64! CODE
65! ----
66! In case of 1D
67! ~~~~~~~~~~~~~
68if (ngrid == 1) then
69    call read_start1D()
70    return
71end if
72
73! In case of 3D
74! ~~~~~~~~~~~~~
75! Open
76inquire(file = start_name,exist = here)
77if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start_name//'"!',1)
78call print_msg('> Reading "'//start_name//'"',LVL_NFO)
79call open_nc(start_name,'read')
80
81! Get the variables
82call get_var_nc('ap',tmp1d)
83call set_ap(tmp1d)
84
85call get_var_nc('bp',tmp1d)
86call set_bp(tmp1d)
87
88call get_var_nc('ps',tmp2d)
89call set_ps_PCM(tmp2d)
90
91call get_var_nc('teta',tmp3d)
92call set_teta_PCM(tmp3d)
93
94do i = 1,nq
95    call get_var_nc(trim(qnames(i)),tmp3d)
96    call set_q_PCM(tmp3d,i)
97end do
98
99! Close
100call close_nc(start_name)
101
102END SUBROUTINE read_start
103!=======================================================================
104
105!=======================================================================
106SUBROUTINE read_start1D()
107!-----------------------------------------------------------------------
108! NAME
109!     read_start1D
110!
111! DESCRIPTION
112!     Read the file "start1D.txt".
113!
114! AUTHORS & DATE
115!     JB Clement, 12/2025
116!
117! NOTES
118!
119!-----------------------------------------------------------------------
120
121! DEPENDENCIES
122! ------------
123use geometry,   only: nlayer
124use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM, set_u_PCM, set_v_PCM, compute_alt_coord, set_ap, set_bp
125use tracers,    only: nq, set_q_PCM
126use stoppage,   only: stop_clean
127use display,    only: print_msg, LVL_NFO
128
129! DECLARATION
130! -----------
131implicit none
132
133! LOCAL VARIABLES
134! ---------------
135integer(di)                     :: i, k, ierr, funit
136character(30)                   :: header
137real(dp), dimension(1,1)        :: tmp
138real(dp), dimension(1,1,nlayer) :: q_tmp, teta_tmp, wind_tmp
139real(dp), dimension(nlayer + 1) :: ap_tmp, bp_tmp
140real(dp)                        :: pa, preff
141logical(k4)                     :: here
142
143! CODE
144! ----
145! Open "start1D.txt"
146inquire(file = start1D_name,exist = here)
147if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//start1D_name//'"!',1)
148call print_msg('> Reading "'//start1D_name//'"',LVL_NFO)
149open(newunit = funit,file = start1D_name,status = "old",action = "read",iostat = ierr)
150if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//start1D_name//'"!',ierr)
151
152! Get the variables
153read(funit,*) header, tmp, pa, preff
154call set_ps_PCM(tmp)
155
156do i = 1,nq
157    read(funit,*,iostat = ierr) header, (q_tmp(1,1,k),k = 1,nlayer)
158    if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'not enough atmospheric layers defined in the file "'//start1D_name//'" for the tracer "'//trim(header)//'"!',1)
159    call set_q_PCM(q_tmp,i)
160end do
161read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer)
162call set_u_PCM(wind_tmp)
163read(funit,*,iostat = ierr) header, (wind_tmp(1,1,k),k = 1,nlayer)
164call set_v_PCM(wind_tmp)
165read(funit,*,iostat = ierr) header, (teta_tmp(1,1,k),k = 1,nlayer)
166call set_teta_PCM(teta_tmp)
167
168! Close
169close(funit)
170
171! Construct altitude coordinates (not stored in "start1D.txt")
172call compute_alt_coord(pa,preff,ap_tmp,bp_tmp)
173call set_ap(ap_tmp)
174call set_bp(bp_tmp)
175
176END SUBROUTINE read_start1D
177!=======================================================================
178
179!=======================================================================
180SUBROUTINE read_startfi()
181!-----------------------------------------------------------------------
182! NAME
183!     read_startfi
184!
185! DESCRIPTION
186!     Read the file "startfi.nc".
187!
188! AUTHORS & DATE
189!     JB Clement, 12/2025
190!
191! NOTES
192!
193!-----------------------------------------------------------------------
194
195! DEPENDENCIES
196! ------------
197use geometry,  only: ngrid, nslope, nsoil_PCM
198use stoppage,  only: stop_clean
199use config,    only: read_control_data
200use slopes,    only: set_def_slope_mean, set_subslope_dist, set_iflat
201use surface,   only: set_albedodat_PCM, set_albedo_PCM, set_emissivity_PCM
202use surf_temp, only: set_tsurf_PCM
203use soil_temp, only: set_tsoil_PCM, set_flux_geo_PCM
204use frost,     only: set_h2ofrost_PCM, set_co2frost_PCM
205use surf_ice,  only: set_is_h2o_perice_PCM, set_co2_perice_PCM
206use soil,      only: set_TI_PCM, set_inertiedat_PCM
207use ice_table, only: set_coef_ssdif_PCM
208use display,   only: print_msg, LVL_NFO
209
210! DECLARATION
211! -----------
212implicit none
213
214! LOCAL VARIABLES
215! ---------------
216real(dp), dimension(:),   allocatable       :: tmp1d
217real(dp), dimension(:,:), allocatable       :: tmp2d
218real(dp), dimension(ngrid,nsoil_PCM,nslope) :: tmp3d
219logical(k4)                                 :: here
220
221! CODE
222! ----
223inquire(file = startfi_name,exist = here)
224if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//startfi_name//'"!',1)
225
226! Allocate the array to store the variables
227call print_msg('> Reading "'//startfi_name//'"',LVL_NFO)
228allocate(tmp1d(nslope + 1),tmp2d(ngrid,nslope))
229
230! Get control data
231call read_control_data()
232
233! Open
234call open_nc(startfi_name,'read')
235
236! Get the variables
237if (nslope > 1) then
238    call get_var_nc('def_slope_mean',tmp1d)
239    call set_def_slope_mean(tmp1d)
240
241    call get_var_nc('subslope_dist',tmp2d)
242    call set_subslope_dist(tmp2d)
243end if
244
245call get_var_nc('flux_geo',tmp2d)
246call set_flux_geo_PCM(tmp2d)
247
248call get_var_nc('h2o_ice',tmp2d)
249where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
250call set_h2ofrost_PCM(tmp2d)
251
252call get_var_nc('co2',tmp2d)
253where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
254call set_co2frost_PCM(tmp2d)
255
256call get_var_nc('perennial_co2ice',tmp2d)
257where (tmp2d < 0.) tmp2d = 0._dp ! Remove unphysical values of surface tracer
258call set_co2_perice_PCM(tmp2d)
259
260deallocate(tmp1d); allocate(tmp1d(ngrid))
261call get_var_nc('watercaptag',tmp1d)
262call set_is_h2o_perice_PCM(tmp1d)
263
264call get_var_nc('albedodat',tmp1d)
265call set_albedodat_PCM(tmp1d)
266
267call get_var_nc('albedo',tmp2d)
268call set_albedo_PCM(tmp2d)
269
270call get_var_nc('emis',tmp2d)
271call set_emissivity_PCM(tmp2d)
272
273call get_var_nc('tsurf',tmp2d)
274call set_tsurf_PCM(tmp2d)
275
276call get_var_nc('tsoil',tmp3d)
277call set_tsoil_PCM(tmp3d)
278
279call get_var_nc('coef_ssdif',tmp2d)
280call set_coef_ssdif_PCM(tmp2d)
281
282deallocate(tmp2d); allocate(tmp2d(ngrid,nsoil_PCM))
283call get_var_nc('inertiedat',tmp2d)
284call set_inertiedat_PCM(tmp2d)
285
286call get_var_nc('inertiesoil',tmp3d)
287call set_TI_PCM(tmp3d)
288
289! Close/Deallocate
290call close_nc(startfi_name)
291deallocate(tmp1d,tmp2d)
292
293END SUBROUTINE read_startfi
294!=======================================================================
295
296!=======================================================================
297SUBROUTINE read_startevo(tsurf_avg_yr1,tsurf_avg_yr2,ps_avg_glob,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,h2o_ice,co2_ice, &
298                         tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map,       &
299                         h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
300!-----------------------------------------------------------------------
301! NAME
302!     read_startevo
303!
304! DESCRIPTION
305!     Read the file "startevo.nc" which stores the PEM state.
306!
307! AUTHORS & DATE
308!     L. Lange, 2023
309!     JB Clement, 2023-2026
310!
311! NOTES
312!
313!-----------------------------------------------------------------------
314
315! DEPENDENCIES
316! ------------
317use stoppage,           only: stop_clean
318use geometry,           only: ngrid, nslope, nsoil, nsoil_PCM
319use evolution,          only: dt
320use physics,            only: mugaz, r, alpha_clap_h2o, beta_clap_h2o
321use frost,              only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM
322use surf_ice,           only: is_h2o_perice_PCM, h2oice_huge_ini, co2_perice_PCM, threshold_h2oice_cap
323use soil,               only: do_soil, TI, depth_breccia, depth_bedrock, index_breccia, index_bedrock, inertiedat, layer, inertiedat_PCM
324use soil_temp,          only: ini_tsoil_pem, compute_tsoil
325use soil_therm_inertia, only: update_soil_TI, TI_breccia, TI_bedrock
326use ice_table,          only: icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium
327use sorption,           only: do_sorption, evolve_regolith_adsorption
328use tracers,            only: mmol, iPCM_qh2o
329use layered_deposits,   only: layering, do_layering, array2map, ini_layering, add_stratum
330use surf_ice,           only: rho_co2ice, rho_h2oice
331use slopes,             only: subslope_dist, def_slope_mean
332use maths,              only: pi
333use display,            only: print_msg, LVL_WRN, LVL_NFO, LVL_ERR
334use utility,            only: int2str
335
336! DECLARATION
337! -----------
338implicit none
339
340! ARGUMENTS
341! ---------
342real(dp),       dimension(:,:),   intent(in)    :: tsurf_avg_yr1, tsurf_avg_yr2 ! Surface temperature at the last but one PCM run and at the last PCM run
343real(dp),                         intent(in)    :: ps_avg_glob                  ! Mean average pressure on the planet
344real(dp),       dimension(:,:),   intent(in)    :: ps_ts                        ! Surface pressure timeseries
345real(dp),       dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts           ! MMR of CO2 and H2O
346real(dp),       dimension(:,:),   intent(in)    :: h2o_surfdensity_avg          ! Average of surface water density
347real(dp),       dimension(:,:),   intent(inout) :: h2o_ice, co2_ice             ! Surface ice
348real(dp),       dimension(:,:,:), intent(inout) :: tsoil_avg                    ! Soil temperature
349real(dp),       dimension(:,:,:), intent(inout) :: h2o_soildensity_avg          ! Average of soil water soil density
350real(dp),       dimension(:,:),   intent(inout) :: icetable_depth               ! Ice table depth
351real(dp),       dimension(:,:),   intent(inout) :: icetable_thickness           ! Ice table thickness
352real(dp),       dimension(:,:,:), intent(inout) :: ice_porefilling              ! Subsurface ice pore filling
353real(dp),       dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg     ! Mass of H2O and CO2 adsorbed
354type(layering), dimension(:,:),   intent(inout) :: layerings_map                ! Layerings
355real(dp),       dimension(:),     intent(inout) :: delta_h2o_ads, delta_co2_ads ! Mass of H2O and CO2 exchanged due to adsorption/desorption
356
357! LOCAL VARIABLES
358! ---------------
359logical(k4)                               :: here
360integer(di)                               :: i, islope, k, nb_str_max, nsoil_startevo
361real(dp)                                  :: delta           ! Depth of the interface regolith/breccia, breccia/bedrock [m]
362real(dp), dimension(ngrid,nsoil,nslope)   :: TI_startevo     ! Soil thermal inertia saved in the startevo [SI]
363real(dp), dimension(ngrid,nsoil,nslope)   :: tsoil_startevo  ! Soil temperature saved in the startevo [K]
364real(dp), dimension(:,:,:,:), allocatable :: layerings_array ! Array for layerings
365
366! CODE
367! ----
368! Check if the file exists
369inquire(file = startevo_name,exist = here)
370
371! If the file is not here
372! ~~~~~~~~~~~~~~~~~~~~~~~
373if (.not. here) then
374    ! It is possibly because it is the very first PEM run so everything needs to be initalized
375    call print_msg('> The file "'//startevo_name//'" was not found (possibly because this is the first PEM run)',LVL_WRN)
376
377    ! H2O ice
378    call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true and where yearly minimum of frost can be considered as a huge reservoir ('threshold_h2oice_cap').",LVL_WRN)
379    h2o_ice(:,:) = 0._dp
380    do i = 1,ngrid
381        if (is_h2o_perice_PCM(i)) then
382            h2o_ice(i,:) = h2oice_huge_ini
383        else if (sum((h2ofrost_PCM(i,:) - h2o_frost4PCM(i,:))*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
384            h2o_ice(i,:) = h2oice_huge_ini
385        end if
386    end do
387
388    ! CO2 ice
389    call print_msg("'co2_ice' is initialized with 'perennial_co2ice' and yearly minimum of frost found in the PCM.",LVL_WRN)
390    co2_ice(:,:) = co2_perice_PCM(:,:) + co2frost_PCM(:,:) - co2_frost4PCM(:,:)
391
392    if (do_soil) then
393        ! Thermal inertia
394        do islope = 1,nslope
395            do i = 1,ngrid
396                if (TI(i,index_breccia,islope) < TI_breccia) then
397                    !!! transition
398                    delta = depth_breccia
399                    TI(i,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                  &
400                                                          (((delta - layer(index_breccia))/(TI(i,index_breccia,islope)**2)) + &
401                                                          ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
402                    do k = index_breccia + 2,index_bedrock
403                        TI(i,k,islope) = TI_breccia
404                    end do
405                else ! we keep the high ti values
406                    do k = index_breccia + 1,index_bedrock
407                        TI(i,k,islope) = TI(i,index_breccia,islope)
408                    end do
409                end if
410                !! transition
411                delta = depth_bedrock
412                TI(i,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                  &
413                                                      (((delta - layer(index_bedrock))/(TI(i,index_bedrock,islope)**2)) + &
414                                                      ((layer(index_bedrock + 1) - delta)/(TI_breccia**2))))
415                do k = index_bedrock + 2,nsoil
416                    TI(i,k,islope) = TI_bedrock
417                end do
418            end do
419        end do
420
421        do k = 1,nsoil_PCM
422            inertiedat(:,k) = inertiedat_PCM(:,k)
423        end do
424        !!! zone de transition
425        delta = depth_breccia
426        do i = 1,ngrid
427            if (inertiedat(i,index_breccia) < TI_breccia) then
428                inertiedat(i,index_breccia + 1) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                   &
429                                                       (((delta - layer(index_breccia))/(inertiedat(i,index_breccia)**2)) + &
430                                                       ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
431                do k = index_breccia + 2,index_bedrock
432                    inertiedat(i,k) = TI_breccia
433                end do
434            else
435                do k = index_breccia + 1,index_bedrock
436                    inertiedat(i,k) = inertiedat(i,index_breccia)
437                end do
438            end if
439        end do
440        !!! zone de transition
441        delta = depth_bedrock
442        do i = 1,ngrid
443            inertiedat(i,index_bedrock + 1) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                   &
444                                                   (((delta - layer(index_bedrock))/(inertiedat(i,index_bedrock)**2)) + &
445                                                   ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
446        end do
447
448        do k = index_bedrock + 2,nsoil
449            do i = 1,ngrid
450                inertiedat(i,k) = TI_bedrock
451            end do
452        end do
453
454        ! Soil temperature
455        do islope = 1,nslope
456            call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
457            call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
458
459            ! First raw initialization
460            h2o_soildensity_avg(:,nsoil_PCM + 1:nsoil,islope) = exp(beta_clap_h2o/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope) + alpha_clap_h2o)/tsoil_avg(:,nsoil_PCM + 1:nsoil,islope)**mmol(iPCM_qh2o)/(mugaz*r)
461        end do
462
463        ! Ice table
464        if (icetable_equilibrium) then
465            call computeice_table_equilibrium(is_h2o_perice_PCM,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness)
466            call update_soil_TI(h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
467            do islope = 1,nslope
468                call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
469            end do
470        else if (icetable_dynamic) then
471            call update_soil_TI(h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
472            do islope = 1,nslope
473                call ini_tsoil_pem(ngrid,nsoil,TI(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_avg(:,:,islope))
474            end do
475        end if
476
477        ! Absorption in regolith
478        if (do_sorption) then
479            h2o_ads_reg = 0._dp
480            co2_ads_reg = 0._dp
481            call evolve_regolith_adsorption(h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
482        end if ! do_sorption
483    end if ! do_soil
484
485    ! Layering
486    if (do_layering) then
487        call print_msg('layerings_map is initialized with sub-surface strata.',LVL_WRN)
488        call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.",LVL_WRN)
489        do i = 1,ngrid
490            if (is_h2o_perice_PCM(i)) then
491                do islope = 1,nslope
492                    call ini_layering(layerings_map(i,islope))
493                    call add_stratum(layerings_map(i,islope),1.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,h2oice_huge_ini/rho_h2oice,0.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,0._dp)
494                end do
495            else
496                do islope = 1,nslope
497                    call ini_layering(layerings_map(i,islope))
498                    if (co2_perice_PCM(i,islope) > 0.) call add_stratum(layerings_map(i,islope),1.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0._dp)
499                end do
500            end if
501        end do
502    end if ! do_layering
503
504    return
505end if
506
507! If the file is present
508! ~~~~~~~~~~~~~~~~~~~~~~
509call print_msg('> Reading "'//startevo_name//'"',LVL_NFO)
510call open_nc(startevo_name,'read')
511
512! H2O ice
513h2o_ice(:,:) = 0._dp
514call get_var_nc('h2o_ice',h2o_ice(:,:))
515
516! CO2 ice
517co2_ice(:,:) = 0._dp
518call get_var_nc('co2_ice',co2_ice(:,:))
519
520if (do_soil) then
521    ! Check if the number of soil layers is compatible
522    call get_dim_nc('subsurface_layers',nsoil_startevo)
523    if (nsoil_startevo /= nsoil) then
524        call print_msg('nsoil (PEM) = '//int2str(nsoil)//' | nsoil ("'//startevo_name//'") = '//int2str(nsoil_startevo),LVL_ERR)
525        call stop_clean(__FILE__,__LINE__,'nsoil defined in the PEM is different from the one read in "'//startevo_name//'"!',1)
526    end if
527
528    do islope = 1,nslope
529        ! Soil temperature
530        call get_var_nc('tsoil',tsoil_startevo(:,:,islope))
531        ! Predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startevo is adapted firstly
532        !                      for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
533        call compute_tsoil(ngrid,nsoil,.true.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startevo(:,:,islope))
534        call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr1(:,islope),tsoil_startevo(:,:,islope))
535        call compute_tsoil(ngrid,nsoil,.false.,TI(:,:,islope),dt,tsurf_avg_yr2(:,islope),tsoil_startevo(:,:,islope))
536    end do
537    tsoil_avg(:,nsoil_PCM + 1:nsoil,:) = tsoil_startevo(:,nsoil_PCM + 1:nsoil,:)
538    if (any(isnan(tsoil_avg(:,:,:)))) call stop_clean(__FILE__,__LINE__,"NaN detected in 'tsoil_avg'",1)
539
540    ! Thermal inertia
541    call get_var_nc('TI',TI_startevo(:,:,:))
542    TI(:,nsoil_PCM + 1:nsoil,:) = TI_startevo(:,nsoil_PCM + 1:nsoil,:) ! 1st layers can change because of the presence of ice at the surface, so we don't change it here
543    call get_var_nc('inertiedat',inertiedat(:,:))
544
545    ! Ice table
546    call get_var_nc('icetable_depth',icetable_depth(:,:))
547    if (icetable_dynamic) call get_var_nc('ice_porefilling',ice_porefilling(:,:,:))
548    if (icetable_equilibrium) then
549        call get_var_nc('icetable_thickness',icetable_thickness(:,:))
550    else if (icetable_dynamic) then
551        call get_var_nc('ice_porefilling',ice_porefilling(:,:,:))
552    end if
553
554    ! Absorption in regolith
555    if (do_sorption) then
556        call get_var_nc('co2_ads_reg',co2_ads_reg(:,:,:))
557        call get_var_nc('h2o_ads_reg',h2o_ads_reg(:,:,:))
558        call evolve_regolith_adsorption(h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
559    end if ! do_sorption
560end if ! do_soil
561
562! Layering
563if (do_layering) then
564    call get_dim_nc('nb_str_max',nb_str_max)
565    allocate(layerings_array(ngrid,nslope,nb_str_max,6))
566    layerings_array = 0.
567    call get_var_nc('stratif_top_elevation',layerings_array(:,:,:,1))
568    call get_var_nc('stratif_h_co2ice',layerings_array(:,:,:,2))
569    call get_var_nc('stratif_h_h2oice',layerings_array(:,:,:,3))
570    call get_var_nc('stratif_h_dust',layerings_array(:,:,:,4))
571    call get_var_nc('stratif_h_pore',layerings_array(:,:,:,5))
572    call get_var_nc('stratif_poreice_volfrac',layerings_array(:,:,:,6))
573    call array2map(layerings_array,layerings_map)
574    deallocate(layerings_array)
575end if
576
577! Close
578call close_nc(startevo_name)
579
580END SUBROUTINE read_startevo
581!=======================================================================
582
583END MODULE clim_state_init
Note: See TracBrowser for help on using the repository browser.