source: trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 @ 4170

Last change on this file since 4170 was 4170, checked in by jbclement, 10 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: 13.6 KB
Line 
1MODULE surf_ice
2!-----------------------------------------------------------------------
3! NAME
4!     surf_ice
5!
6! DESCRIPTION
7!     Surface ice management.
8!
9! AUTHORS & DATE
10!     JB Clement, 12/2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use numerics, only: dp, qp, di, k4
19
20! DECLARATION
21! -----------
22implicit none
23
24! PARAMETERS
25! ----------
26real(dp),                                 parameter :: rho_co2ice = 1650._dp ! Density of CO2 ice [kg.m-3]
27real(dp),                                 parameter :: rho_h2oice = 920._dp  ! Density of H2O ice [kg.m-3]
28real(dp),                                 protected :: threshold_h2oice_cap  ! Threshold to consider H2O ice as infinite reservoir [kg/m2]
29real(dp),                                 protected :: h2oice_huge_ini       ! Initial value for huge reservoirs of H2O ice [kg/m^2]
30logical(k4), dimension(:),   allocatable, protected :: is_h2o_perice_PCM     ! Flag for the presence of perennial H2O ice in the PCM at the beginning [kg/m2]
31real(dp),    dimension(:,:), allocatable, protected :: co2_perice_PCM        ! Perennial CO2 ice in the PCM at the beginning [kg/m2]
32
33contains
34!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35
36!=======================================================================
37SUBROUTINE ini_surf_ice()
38!-----------------------------------------------------------------------
39! NAME
40!     ini_surf_ice
41!
42! DESCRIPTION
43!     Initialize surface ice arrays.
44!
45! AUTHORS & DATE
46!     JB Clement 12/2025
47!
48! NOTES
49!
50!-----------------------------------------------------------------------
51
52! DEPENDENCIES
53! ------------
54use geometry, only: ngrid, nslope
55
56! DECLARATION
57! -----------
58implicit none
59
60! CODE
61! ----
62allocate(is_h2o_perice_PCM(ngrid))
63allocate(co2_perice_PCM(ngrid,nslope))
64is_h2o_perice_PCM(:) = .false.
65co2_perice_PCM(:,:) = 0._dp
66
67END SUBROUTINE ini_surf_ice
68!=======================================================================
69
70!=======================================================================
71SUBROUTINE end_surf_ice()
72!-----------------------------------------------------------------------
73! NAME
74!     end_surf_ice
75!
76! DESCRIPTION
77!     Deallocate surface ice arrays.
78!
79! AUTHORS & DATE
80!     JB Clement, 12/2025
81!
82! NOTES
83!
84!-----------------------------------------------------------------------
85
86! DECLARATION
87! -----------
88implicit none
89
90! CODE
91! ----
92if (allocated(is_h2o_perice_PCM)) deallocate(is_h2o_perice_PCM)
93if (allocated(co2_perice_PCM)) deallocate(co2_perice_PCM)
94
95END SUBROUTINE end_surf_ice
96!=======================================================================
97
98!=======================================================================
99SUBROUTINE set_surf_ice_config(threshold_h2oice_cap_in,h2oice_huge_ini_in)
100!-----------------------------------------------------------------------
101! NAME
102!     set_surf_ice_config
103!
104! DESCRIPTION
105!     Setter for "surf_ice" configuration parameters.
106!
107! AUTHORS & DATE
108!     JB Clement, 02/2026
109!
110! NOTES
111!
112!-----------------------------------------------------------------------
113
114! DEPENDENCIES
115! ------------
116use utility,  only: real2str
117use display,  only: print_msg, LVL_NFO
118use stoppage, only: stop_clean
119
120! DECLARATION
121! -----------
122implicit none
123
124! ARGUMENTS
125! ---------
126real(dp), intent(in) :: threshold_h2oice_cap_in, h2oice_huge_ini_in
127
128! CODE
129! ----
130threshold_h2oice_cap = threshold_h2oice_cap_in
131h2oice_huge_ini = h2oice_huge_ini_in
132call print_msg('threshold_h2oice_cap = '//real2str(threshold_h2oice_cap),LVL_NFO)
133call print_msg('h2oice_huge_ini      = '//real2str(h2oice_huge_ini),LVL_NFO)
134if (threshold_h2oice_cap < 0._dp) call stop_clean(__FILE__,__LINE__,'''threshold_h2oice_cap'' must be positive or zero!',1)
135if (h2oice_huge_ini < 0._dp) call stop_clean(__FILE__,__LINE__,'''h2oice_huge_ini'' must be positive or zero!',1)
136
137END SUBROUTINE set_surf_ice_config
138!=======================================================================
139
140!=======================================================================
141SUBROUTINE set_co2_perice_PCM(co2_perice_PCM_in)
142!-----------------------------------------------------------------------
143! NAME
144!     set_co2_perice_PCM
145!
146! DESCRIPTION
147!     Setter for 'co2_perice_PCM'.
148!
149! AUTHORS & DATE
150!     JB Clement, 12/2025
151!
152! NOTES
153!
154!-----------------------------------------------------------------------
155
156! DECLARATION
157! -----------
158implicit none
159
160! ARGUMENTS
161! ---------
162real(dp), dimension(:,:), intent(in) :: co2_perice_PCM_in
163
164! CODE
165! ----
166co2_perice_PCM(:,:) = co2_perice_PCM_in(:,:)
167
168END SUBROUTINE set_co2_perice_PCM
169!=======================================================================
170
171!=======================================================================
172SUBROUTINE set_is_h2o_perice_PCM(is_h2o_perice_PCM_in)
173!-----------------------------------------------------------------------
174! NAME
175!     set_is_h2o_perice_PCM
176!
177! DESCRIPTION
178!     Setter for 'is_h2o_perice_PCM'.
179!
180! AUTHORS & DATE
181!     JB Clement, 12/2025
182!
183! NOTES
184!
185!-----------------------------------------------------------------------
186
187! DECLARATION
188! -----------
189implicit none
190
191! ARGUMENTS
192! ---------
193real(dp), dimension(:), intent(in) :: is_h2o_perice_PCM_in
194
195! CODE
196! ----
197where(is_h2o_perice_PCM_in > 0.5_dp)
198    is_h2o_perice_PCM = .true.
199else where
200    is_h2o_perice_PCM = .false.
201end where
202
203END SUBROUTINE set_is_h2o_perice_PCM
204!=======================================================================
205
206!=======================================================================
207SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM)
208!-----------------------------------------------------------------------
209! NAME
210!     build4PCM_perice
211!
212! DESCRIPTION
213!     Reconstructs perennial ice and frost for the PCM.
214!
215! AUTHORS & DATE
216!     JB Clement, 12/2025
217!
218! NOTES
219!
220!-----------------------------------------------------------------------
221
222! DEPENDENCIES
223! ------------
224use geometry, only: ngrid
225use frost,    only: h2o_frost4PCM
226use slopes,   only: subslope_dist, def_slope_mean
227use maths,    only: pi
228use display,  only: print_msg, LVL_NFO
229
230! DECLARATION
231! -----------
232implicit none
233
234! ARGUMENTS
235! ---------
236real(dp),    dimension(:,:), intent(inout) :: h2o_ice, co2_ice
237logical(k4), dimension(:),   intent(out)   :: is_h2o_perice            ! H2O perennial ice flag
238real(dp),    dimension(:,:), intent(out)   :: h2o_ice4PCM, co2_ice4PCM ! Ice for PCM
239
240! LOCAL VARIABLES
241! ---------------
242integer(di) :: i
243
244! CODE
245! ----
246call print_msg('> Building surface ice/frost for the PCM',LVL_NFO)
247co2_ice4PCM(:,:) = co2_ice(:,:)
248h2o_ice4PCM(:,:) = 0._dp ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
249do i = 1,ngrid
250    ! Is H2O ice still considered as an infinite reservoir for the PCM?
251    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
252        ! There is enough ice to be considered as an infinite reservoir
253        is_h2o_perice(i) = .true.
254    else
255        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
256        is_h2o_perice(i) = .false.
257        h2o_frost4PCM(i,:) = h2o_frost4PCM(i,:) + h2o_ice(i,:)
258        h2o_ice(i,:) = 0._dp
259    end if
260end do
261
262END SUBROUTINE build4PCM_perice
263!=======================================================================
264
265!=======================================================================
266SUBROUTINE evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
267!-----------------------------------------------------------------------
268! NAME
269!     evolve_co2ice
270!
271! DESCRIPTION
272!     Compute the evolution of CO2 ice over one year.
273!
274! AUTHORS & DATE
275!     R. Vandemeulebrouck
276!     L. Lange
277!     JB Clement, 2023-2025
278!
279! NOTES
280!
281!-----------------------------------------------------------------------
282
283! DEPENDENCIES
284! ------------
285use geometry,  only: ngrid, nslope
286use evolution, only: dt
287use display,   only: print_msg, LVL_NFO
288
289! DECLARATION
290! -----------
291implicit none
292
293! ARGUMENTS
294! ---------
295real(dp), dimension(:,:), intent(inout) :: co2_ice
296real(dp), dimension(:,:), intent(inout) :: d_co2ice
297real(dp), dimension(:,:), intent(out)   :: zshift_surf
298
299! LOCAL VARIABLES
300! ---------------
301real(dp), dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
302
303! CODE
304! ----
305! Evolution of CO2 ice for each physical point
306call print_msg('> Evolution of CO2 ice',LVL_NFO)
307
308co2_ice_old(:,:) = co2_ice(:,:)
309co2_ice(:,:) = co2_ice(:,:) + d_co2ice(:,:)*dt
310where (co2_ice < 0._dp)
311    co2_ice(:,:) = 0._dp
312    d_co2ice(:,:) = -co2_ice_old(:,:)/dt
313end where
314
315zshift_surf(:,:) = zshift_surf(:,:) + d_co2ice(:,:)*dt/rho_co2ice
316
317END SUBROUTINE evolve_co2ice
318!=======================================================================
319
320!=======================================================================
321SUBROUTINE evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
322!-----------------------------------------------------------------------
323! NAME
324!     evolve_h2oice
325!
326! DESCRIPTION
327!     Compute the evolution of H2O ice over one year.
328!
329! AUTHORS & DATE
330!     R. Vandemeulebrouck
331!     L. Lange
332!     JB Clement, 2023-2025
333!
334! NOTES
335!
336!-----------------------------------------------------------------------
337
338! DEPENDENCIES
339! ------------
340use geometry,      only: ngrid, nslope
341use evolution,     only: dt
342use stopping_crit, only: stopping_crit_h2o, stopFlags
343use display,       only: print_msg, LVL_NFO
344
345! DECLARATION
346! -----------
347implicit none
348
349! ARGUMENTS
350! ---------
351real(dp), dimension(:),   intent(in)    :: delta_h2o_ads, delta_icetable
352real(dp), dimension(:,:), intent(inout) :: h2o_ice, d_h2oice
353type(stopFlags),          intent(inout) :: stopcrit
354real(dp), dimension(:,:), intent(out)   :: zshift_surf
355
356! LOCAL VARIABLES
357! ---------------
358real(qp)                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! H2O balance variables
359real(dp), dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
360
361! CODE
362! ----
363call print_msg('> Evolution of H2O ice',LVL_NFO)
364
365if (ngrid == 1) then ! In 1D
366    where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything
367        ! We sublimate what we can
368        d_h2oice_new(:,:) = h2o_ice(:,:)/dt
369        ! It means the tendency is zero next time
370        d_h2oice(:,:) = 0._dp
371    else where
372        d_h2oice_new(:,:) = d_h2oice(:,:)
373    end where
374else ! In 3D
375    call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
376    if (stopcrit%stop_code() > 0) return
377    call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
378end if
379
380h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt
381zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice_new(:,:)*dt/rho_h2oice
382
383END SUBROUTINE evolve_h2oice
384!=======================================================================
385
386!=======================================================================
387SUBROUTINE balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_balanced)
388!-----------------------------------------------------------------------
389! NAME
390!     balance_h2oice_reservoirs
391!
392! DESCRIPTION
393!     Balance H2O flux from and into atmosphere across reservoirs.
394!
395! AUTHORS & DATE
396!     JB Clement, 2025
397!
398! NOTES
399!
400!-----------------------------------------------------------------------
401
402! DEPENDENCIES
403! ------------
404use evolution, only: dt
405use geometry,  only: ngrid, nslope
406
407! DECLARATION
408! -----------
409implicit none
410
411! ARGUMENTS
412! ---------
413real(qp),                 intent(in)    :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm
414real(dp), dimension(:,:), intent(in)    :: h2o_ice
415real(dp), dimension(:,:), intent(inout) :: d_h2oice
416real(dp), dimension(:,:), intent(out)   :: d_balanced
417
418! LOCAL VARIABLES
419! ---------------
420integer(di) :: i, islope
421real(qp)    :: S_corr_subl, S_corr_cond, S_target_subl, S_target_cond, S_ghostice ! Balance variables
422real(dp)    :: d_target
423
424! CODE
425! ----
426! Compute global targets
427S_corr_cond = (S_h2o_2_atm - S_atm_2_h2o)/2._qp ! Correction = deviation to the mean
428S_corr_subl = (S_atm_2_h2o - S_h2o_2_atm)/2._qp ! Correction = deviation to the mean
429S_target_cond = abs(S_atm_2_h2oice + S_corr_cond)
430S_target_subl = abs(S_h2oice_2_atm + S_corr_subl)
431
432! Compute balanced tendencies with positivity limiter
433d_balanced(:,:) = 0._dp
434S_ghostice = 0._qp
435do i = 1,ngrid
436    do islope = 1,nslope
437        if (d_h2oice(i,islope) > 0._dp) then ! Condensation
438            d_balanced(i,islope) = d_h2oice(i,islope)*real(S_target_cond/S_atm_2_h2oice,dp)
439        else if (d_h2oice(i,islope) < 0._dp) then ! Sublimation
440            d_target = d_h2oice(i,islope)*real(S_target_subl/S_h2oice_2_atm,dp)
441            if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate
442                d_balanced(i,islope) = d_target
443            else ! Not enough ice to sublimate
444                ! Sublimate all the available ice
445                d_balanced(i,islope) = -h2o_ice(i,islope)/dt
446                ! If fully depleted, zero tendency for next time
447                d_h2oice(i,islope) = 0._dp
448                ! Compute the amount of H2O ice unable to sublimate
449                S_ghostice = S_ghostice + real(abs(d_target*dt),qp) - real(h2o_ice(i,islope),qp)
450            end if
451        end if
452    end do
453end do
454
455! Enforce conservation removing the ghost ice from places where ice condensed
456where (d_balanced(:,:) > 0._dp) d_balanced(:,:) = d_balanced(:,:)*real((S_target_cond - S_ghostice)/S_target_cond,dp)
457
458END SUBROUTINE balance_h2oice_reservoirs
459!=======================================================================
460
461END MODULE surf_ice
Note: See TracBrowser for help on using the repository browser.