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

Last change on this file since 4074 was 4071, checked in by jbclement, 12 days ago

PEM:

  • Making the computation of ice tendencies more reliable by doing it after 'read_startpem' to know exactly where perennial ice is (no matter if there is a "startpem.nc" or not). Moreover, when there is no "startpem.nc", location of perennial ice depends now on 'watercaptag' and on huge amount of frost. This prevents negative ice tendency while there is no ice which can happen with weird PCM inputs (i.e. 'watercaptag' = F & 'watercap' /= 0 & no "stratpem.nc").
  • Few small safeguards throughout the code.

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