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

Last change on this file since 4177 was 4174, checked in by jbclement, 5 days ago

PEM:
Refactor of H2O flux balancing:

  • centralizes flux balancing logic;
  • removes intermediate variables, especially related to of H2O mass bookkeeping;
  • adds new stopping criterion when H2O flux balance fails8 (sanity check);
  • introduces configurable weighting strategies for H2O flux balancing.

JBC
File size: 18.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, minieps
19
20! DECLARATION
21! -----------
22implicit none
23
24! PARAMETERS
25! ----------
26integer(di),                              parameter :: WEIGHT_UNIFORM   = 1_di ! Equal weight for all adjustable points
27integer(di),                              parameter :: WEIGHT_ABSFLUX   = 2_di ! Weight by absolute local H2O surface flux
28integer(di),                              parameter :: WEIGHT_CAPACITY  = 3_di ! Weight by available local adjustment range
29integer(di),                              parameter :: WEIGHT_COMBINED  = 4_di ! Capacity weighted by current absolute flux
30integer(di),                              parameter :: WEIGHT_RESERVOIR = 5_di ! Weight by locally available H2O ice mass
31integer(di),                              parameter :: WEIGHT_DIRECTION = 6_di ! Weight by directional capacity (up/down)
32integer(di),                              private   :: weight_method = WEIGHT_DIRECTION ! Default method for balancing
33real(dp),                                 parameter :: rho_co2ice = 1650._dp ! Density of CO2 ice [kg.m-3]
34real(dp),                                 parameter :: rho_h2oice = 920._dp  ! Density of H2O ice [kg.m-3]
35real(dp),                                 protected :: threshold_h2oice_cap  ! Threshold to consider H2O ice as infinite reservoir [kg/m2]
36real(dp),                                 protected :: h2oice_huge_ini       ! Initial value for huge reservoirs of H2O ice [kg/m^2]
37logical(k4), dimension(:),   allocatable, protected :: is_h2o_perice_PCM     ! Flag for the presence of perennial H2O ice in the PCM at the beginning [kg/m2]
38real(dp),    dimension(:,:), allocatable, protected :: co2_perice_PCM        ! Perennial CO2 ice in the PCM at the beginning [kg/m2]
39
40contains
41!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
42
43!=======================================================================
44SUBROUTINE ini_surf_ice()
45!-----------------------------------------------------------------------
46! NAME
47!     ini_surf_ice
48!
49! DESCRIPTION
50!     Initialize surface ice arrays.
51!
52! AUTHORS & DATE
53!     JB Clement 12/2025
54!
55! NOTES
56!
57!-----------------------------------------------------------------------
58
59! DEPENDENCIES
60! ------------
61use geometry, only: ngrid, nslope
62
63! DECLARATION
64! -----------
65implicit none
66
67! CODE
68! ----
69allocate(is_h2o_perice_PCM(ngrid))
70allocate(co2_perice_PCM(ngrid,nslope))
71is_h2o_perice_PCM(:) = .false.
72co2_perice_PCM(:,:) = 0._dp
73
74END SUBROUTINE ini_surf_ice
75!=======================================================================
76
77!=======================================================================
78SUBROUTINE end_surf_ice()
79!-----------------------------------------------------------------------
80! NAME
81!     end_surf_ice
82!
83! DESCRIPTION
84!     Deallocate surface ice arrays.
85!
86! AUTHORS & DATE
87!     JB Clement, 12/2025
88!
89! NOTES
90!
91!-----------------------------------------------------------------------
92
93! DECLARATION
94! -----------
95implicit none
96
97! CODE
98! ----
99if (allocated(is_h2o_perice_PCM)) deallocate(is_h2o_perice_PCM)
100if (allocated(co2_perice_PCM)) deallocate(co2_perice_PCM)
101
102END SUBROUTINE end_surf_ice
103!=======================================================================
104
105!=======================================================================
106SUBROUTINE set_surf_ice_config(threshold_h2oice_cap_in,h2oice_huge_ini_in)
107!-----------------------------------------------------------------------
108! NAME
109!     set_surf_ice_config
110!
111! DESCRIPTION
112!     Setter for "surf_ice" configuration parameters.
113!
114! AUTHORS & DATE
115!     JB Clement, 02/2026
116!
117! NOTES
118!
119!-----------------------------------------------------------------------
120
121! DEPENDENCIES
122! ------------
123use utility,  only: real2str
124use display,  only: print_msg, LVL_NFO
125use stoppage, only: stop_clean
126
127! DECLARATION
128! -----------
129implicit none
130
131! ARGUMENTS
132! ---------
133real(dp), intent(in) :: threshold_h2oice_cap_in, h2oice_huge_ini_in
134
135! CODE
136! ----
137threshold_h2oice_cap = threshold_h2oice_cap_in
138h2oice_huge_ini = h2oice_huge_ini_in
139call print_msg('threshold_h2oice_cap = '//real2str(threshold_h2oice_cap),LVL_NFO)
140call print_msg('h2oice_huge_ini      = '//real2str(h2oice_huge_ini),LVL_NFO)
141if (threshold_h2oice_cap < 0._dp) call stop_clean(__FILE__,__LINE__,'''threshold_h2oice_cap'' must be positive or zero!',1)
142if (h2oice_huge_ini < 0._dp) call stop_clean(__FILE__,__LINE__,'''h2oice_huge_ini'' must be positive or zero!',1)
143
144END SUBROUTINE set_surf_ice_config
145!=======================================================================
146
147!=======================================================================
148SUBROUTINE set_co2_perice_PCM(co2_perice_PCM_in)
149!-----------------------------------------------------------------------
150! NAME
151!     set_co2_perice_PCM
152!
153! DESCRIPTION
154!     Setter for 'co2_perice_PCM'.
155!
156! AUTHORS & DATE
157!     JB Clement, 12/2025
158!
159! NOTES
160!
161!-----------------------------------------------------------------------
162
163! DECLARATION
164! -----------
165implicit none
166
167! ARGUMENTS
168! ---------
169real(dp), dimension(:,:), intent(in) :: co2_perice_PCM_in
170
171! CODE
172! ----
173co2_perice_PCM(:,:) = co2_perice_PCM_in(:,:)
174
175END SUBROUTINE set_co2_perice_PCM
176!=======================================================================
177
178!=======================================================================
179SUBROUTINE set_is_h2o_perice_PCM(is_h2o_perice_PCM_in)
180!-----------------------------------------------------------------------
181! NAME
182!     set_is_h2o_perice_PCM
183!
184! DESCRIPTION
185!     Setter for 'is_h2o_perice_PCM'.
186!
187! AUTHORS & DATE
188!     JB Clement, 12/2025
189!
190! NOTES
191!
192!-----------------------------------------------------------------------
193
194! DECLARATION
195! -----------
196implicit none
197
198! ARGUMENTS
199! ---------
200real(dp), dimension(:), intent(in) :: is_h2o_perice_PCM_in
201
202! CODE
203! ----
204where(is_h2o_perice_PCM_in > 0.5_dp)
205    is_h2o_perice_PCM = .true.
206else where
207    is_h2o_perice_PCM = .false.
208end where
209
210END SUBROUTINE set_is_h2o_perice_PCM
211!=======================================================================
212
213!=======================================================================
214SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM)
215!-----------------------------------------------------------------------
216! NAME
217!     build4PCM_perice
218!
219! DESCRIPTION
220!     Reconstructs perennial ice and frost for the PCM.
221!
222! AUTHORS & DATE
223!     JB Clement, 12/2025
224!
225! NOTES
226!
227!-----------------------------------------------------------------------
228
229! DEPENDENCIES
230! ------------
231use geometry, only: ngrid
232use frost,    only: h2o_frost4PCM
233use slopes,   only: subslope_dist, def_slope_mean
234use maths,    only: pi
235use display,  only: print_msg, LVL_NFO
236
237! DECLARATION
238! -----------
239implicit none
240
241! ARGUMENTS
242! ---------
243real(dp),    dimension(:,:), intent(inout) :: h2o_ice, co2_ice
244logical(k4), dimension(:),   intent(out)   :: is_h2o_perice            ! H2O perennial ice flag
245real(dp),    dimension(:,:), intent(out)   :: h2o_ice4PCM, co2_ice4PCM ! Ice for PCM
246
247! LOCAL VARIABLES
248! ---------------
249integer(di) :: i
250
251! CODE
252! ----
253call print_msg('> Building surface ice/frost for the PCM',LVL_NFO)
254co2_ice4PCM(:,:) = co2_ice(:,:)
255h2o_ice4PCM(:,:) = 0._dp ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
256do i = 1,ngrid
257    ! Is H2O ice still considered as an infinite reservoir for the PCM?
258    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
259        ! There is enough ice to be considered as an infinite reservoir
260        is_h2o_perice(i) = .true.
261    else
262        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
263        is_h2o_perice(i) = .false.
264        h2o_frost4PCM(i,:) = h2o_frost4PCM(i,:) + h2o_ice(i,:)
265        h2o_ice(i,:) = 0._dp
266    end if
267end do
268
269END SUBROUTINE build4PCM_perice
270!=======================================================================
271
272!=======================================================================
273SUBROUTINE evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
274!-----------------------------------------------------------------------
275! NAME
276!     evolve_co2ice
277!
278! DESCRIPTION
279!     Compute the evolution of CO2 ice over one year.
280!
281! AUTHORS & DATE
282!     R. Vandemeulebrouck
283!     L. Lange
284!     JB Clement, 2023-2025
285!
286! NOTES
287!
288!-----------------------------------------------------------------------
289
290! DEPENDENCIES
291! ------------
292use geometry,  only: ngrid, nslope
293use evolution, only: dt
294use display,   only: print_msg, LVL_NFO
295
296! DECLARATION
297! -----------
298implicit none
299
300! ARGUMENTS
301! ---------
302real(dp), dimension(:,:), intent(inout) :: co2_ice
303real(dp), dimension(:,:), intent(inout) :: d_co2ice
304real(dp), dimension(:,:), intent(out)   :: zshift_surf
305
306! LOCAL VARIABLES
307! ---------------
308real(dp), dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
309
310! CODE
311! ----
312! Evolution of CO2 ice for each physical point
313call print_msg('> Evolution of CO2 ice',LVL_NFO)
314
315co2_ice_old(:,:) = co2_ice(:,:)
316co2_ice(:,:) = co2_ice(:,:) + d_co2ice(:,:)*dt
317where (co2_ice < 0._dp)
318    co2_ice(:,:) = 0._dp
319    d_co2ice(:,:) = -co2_ice_old(:,:)/dt
320end where
321
322zshift_surf(:,:) = zshift_surf(:,:) + d_co2ice(:,:)*dt/rho_co2ice
323
324END SUBROUTINE evolve_co2ice
325!=======================================================================
326
327!=======================================================================
328SUBROUTINE evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
329!-----------------------------------------------------------------------
330! NAME
331!     evolve_h2oice
332!
333! DESCRIPTION
334!     Compute the evolution of H2O ice over one year.
335!
336! AUTHORS & DATE
337!     R. Vandemeulebrouck
338!     L. Lange
339!     JB Clement, 2023-2025
340!
341! NOTES
342!
343!-----------------------------------------------------------------------
344
345! DEPENDENCIES
346! ------------
347use geometry,      only: ngrid
348use evolution,     only: dt
349use stopping_crit, only: stopFlags
350use display,       only: print_msg, LVL_NFO
351
352! DECLARATION
353! -----------
354implicit none
355
356! ARGUMENTS
357! ---------
358real(dp), dimension(:),   intent(in)    :: delta_h2o_ads, delta_icetable
359real(dp), dimension(:,:), intent(inout) :: h2o_ice, d_h2oice
360type(stopFlags),          intent(inout) :: stopcrit
361real(dp), dimension(:,:), intent(out)   :: zshift_surf
362
363! CODE
364! ----
365call print_msg('> Evolution of H2O ice',LVL_NFO)
366
367if (ngrid == 1) then ! In 1D
368    where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything
369        ! We sublimate what we can
370        d_h2oice(:,:) = -h2o_ice(:,:)/dt
371    end where
372else ! In 3D
373    call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit)
374    if (stopcrit%stop_code() > 0) return
375end if
376
377h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice(:,:)*dt
378zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice(:,:)*dt/rho_h2oice
379
380END SUBROUTINE evolve_h2oice
381!=======================================================================
382
383!=======================================================================
384SUBROUTINE balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit,weight_type)
385!-----------------------------------------------------------------------
386! NAME
387!     balance_h2o_fluxes
388!
389! DESCRIPTION
390!     Balance H2O fluxes from and into atmosphere across reservoirs.
391!
392! AUTHORS & DATE
393!     JB Clement, 04/2025
394!
395! NOTES
396!
397!-----------------------------------------------------------------------
398
399! DEPENDENCIES
400! ------------
401use evolution,     only: dt
402use geometry,      only: ngrid, nslope, cell_area
403use stopping_crit, only: stopFlags
404use utility,       only: real2str, int2str
405use display,       only: print_msg, LVL_WRN
406
407! DECLARATION
408! -----------
409implicit none
410
411! ARGUMENTS
412! ---------
413real(dp), dimension(:),   intent(in)           :: delta_h2o_ads, delta_icetable
414real(dp), dimension(:,:), intent(in)           :: h2o_ice
415integer(di),              intent(in), optional :: weight_type
416real(dp), dimension(:,:), intent(inout)        :: d_h2oice
417type(stopFlags),          intent(inout)        :: stopcrit
418
419! LOCAL VARIABLES
420! ---------------
421integer(di), parameter :: max_iter = 50_di ! Maximum number of iterations for the balancing procedure
422real(dp),    parameter :: eps = 1.e-12_dp ! Small number to prevent division by zero in weights normalization
423real(dp),    parameter :: tiny_corr = minieps ! Minimum correction to consider that progress is made in the balancing procedure
424integer(di)                              :: i, islope, iter
425integer(di)                              :: method_used
426real(dp)                                 :: flux_scale, tol_flux
427real(qp)                                 :: B_flux, B_flux_surfice, B_flux_nonsurfice, Delta, Delta_used, wsum
428real(dp),    dimension(:,:), allocatable :: w, d_corr, flux_new
429real(dp),    dimension(:,:), allocatable :: flux_min, flux_max
430real(dp),    dimension(:,:), allocatable :: capacity, capacity_up, capacity_down
431logical(k4), dimension(:,:), allocatable :: active, adjustable
432logical(k4), dimension(:,:), allocatable :: clipped_high, clipped_low
433
434! CODE
435! ----
436! Initialization
437allocate(w(ngrid,nslope),d_corr(ngrid,nslope),flux_new(ngrid,nslope))
438allocate(flux_min(ngrid,nslope),flux_max(ngrid,nslope))
439allocate(capacity(ngrid,nslope),capacity_up(ngrid,nslope),capacity_down(ngrid,nslope))
440allocate(active(ngrid,nslope),adjustable(ngrid,nslope))
441allocate(clipped_high(ngrid,nslope),clipped_low(ngrid,nslope))
442
443! Build physical bounds and geometry weights
444do i = 1,ngrid
445    do islope = 1,nslope
446        flux_min(i,islope) = -h2o_ice(i,islope)/dt
447        flux_max(i,islope) = huge(1._dp)
448        if (d_h2oice(i,islope) > 0._dp) then
449            flux_min(i,islope) = 0._dp
450        else if (d_h2oice(i,islope) < 0._dp) then
451            flux_max(i,islope) = 0._dp
452        end if
453        capacity(i,islope) = flux_max(i,islope) - flux_min(i,islope)
454    end do
455end do
456active(:,:) = .true.
457
458! Choose weighting method
459method_used = weight_method
460if (present(weight_type)) then
461    if (weight_type >= WEIGHT_UNIFORM .and. weight_type <= WEIGHT_DIRECTION) then
462        weight_method = weight_type
463        method_used = weight_type
464    else
465        call print_msg('Unknown H2O flux weighting method. Falling back to WEIGHT_DIRECTION.',LVL_WRN)
466        method_used = WEIGHT_DIRECTION
467    end if
468end if
469
470! Compute initial imbalance
471B_flux_nonsurfice = 0._qp
472B_flux_surfice = 0._qp
473do i = 1,ngrid
474    B_flux_nonsurfice = B_flux_nonsurfice + real(delta_h2o_ads(i) + delta_icetable(i),qp)/real(cell_area(i),qp)/real(dt,qp)
475    do islope = 1,nslope
476        B_flux_surfice = B_flux_surfice + real(d_h2oice(i,islope),qp)
477    end do
478end do
479B_flux = B_flux_nonsurfice + B_flux_surfice
480Delta = -B_flux
481flux_scale = max(1._dp,abs(real(B_flux,dp))) ! Scale flux imbalance to get a relevant tolerance for the balancing procedure (e.g., if imbalance is very small, we don't want to require an even smaller imbalance to stop iterating)
482tol_flux = max(1.e-12_dp,1.e-10_dp*flux_scale) ! Tolerance for closing the flux balance, scaled to the initial imbalance
483
484! Iterative correction
485do iter = 1,max_iter
486    ! Stopping criterion
487    if (abs(Delta) <= real(tol_flux,qp)) exit
488
489    ! Prepare correction
490    w(:,:) = 0._dp
491    ! Compute directional capacities
492    capacity_up(:,:) = flux_max(:,:) - d_h2oice(:,:)
493    capacity_down(:,:) = d_h2oice(:,:) - flux_min(:,:)
494
495    ! Compute adjustable mask
496    adjustable(:,:) = active(:,:)
497    if (Delta > 0._qp) then
498        where (capacity_up <= eps) adjustable = .false.
499    else
500        where (capacity_down <= eps) adjustable = .false.
501    end if
502    if (.not. any(adjustable)) exit
503
504    ! Compute weights
505    select case (method_used)
506        case (WEIGHT_UNIFORM) ! Flat redistribution across adjustable cells
507            where (adjustable) w = 1._dp
508        case (WEIGHT_ABSFLUX) ! Prioritize cells already carrying strong fluxes
509            where (adjustable) w = abs(d_h2oice)
510        case (WEIGHT_CAPACITY) ! Prioritize cells with large admissible adjustment room
511            where (adjustable) w = capacity
512        case (WEIGHT_COMBINED) ! Mix current activity (abs flux) and admissible room
513            where (adjustable) w = capacity*max(abs(d_h2oice),eps)
514        case (WEIGHT_RESERVOIR) ! Favor points with larger local ice reservoirs
515            where (adjustable) w = max(h2o_ice,eps)
516        case (WEIGHT_DIRECTION) ! Follow directional room needed by the sign of Delta
517            if (Delta > 0._qp) then
518                where (adjustable) w = max(capacity_up,eps)
519            else
520                where (adjustable) w = max(capacity_down,eps)
521            end if
522        case default
523            where (adjustable) w = 1._dp
524    end select
525
526    wsum = sum(real(w,qp),mask = adjustable)
527    if (wsum <= 0._qp) exit ! Fallback all weights are zero
528
529    ! Compute correction
530    d_corr(:,:) = 0._dp
531    where (adjustable) d_corr = real(Delta*real(w,qp)/wsum,dp)
532    flux_new(:,:) = d_h2oice(:,:) + d_corr(:,:)
533
534    ! Apply bounds (clipping)
535    clipped_high(:,:) = .false.
536    clipped_low(:,:) = .false.
537    where (adjustable .and. flux_new > flux_max)
538        d_corr = flux_max - d_h2oice
539        d_h2oice = flux_max
540        clipped_high = .true.
541    end where
542    where (adjustable .and. flux_new < flux_min)
543        d_corr = flux_min - d_h2oice
544        d_h2oice = flux_min
545        clipped_low = .true.
546    end where
547    where (adjustable .and. .not. (clipped_high .or. clipped_low))
548        d_h2oice = flux_new
549    end where
550    where (clipped_high .or. clipped_low)
551        active = .false.
552    end where
553
554    ! Update imbalance
555    Delta_used = sum(real(d_corr,qp))
556    Delta = Delta - Delta_used
557
558    if (abs(Delta_used) < tiny_corr) exit
559end do
560
561! Diagnostics
562if (abs(Delta) > real(tol_flux,qp)) then
563    call print_msg('Unable to close H2O flux balance on surface-ice reservoirs.',LVL_WRN)
564    call print_msg('Weighting method used = '//int2str(method_used),LVL_WRN)
565    call print_msg('Residual imbalance [kg/m2/y] = '//real2str(Delta),LVL_WRN)
566    call print_msg('Initial imbalance  [kg/m2/y] = '//real2str(B_flux),LVL_WRN)
567    call print_msg('Initial surface ice flux     [kg/m2/y] = '//real2str(B_flux_surfice),LVL_WRN)
568    call print_msg('Initial non-surface ice flux [kg/m2/y] = '//real2str(B_flux_nonsurfice),LVL_WRN)
569    stopcrit%h2o_flux_balance_unclosed = .true.
570end if
571
572! Cleanup
573deallocate(w,d_corr,flux_new,flux_min,flux_max)
574deallocate(capacity,capacity_up,capacity_down)
575deallocate(active,adjustable,clipped_high,clipped_low)
576
577END SUBROUTINE balance_h2o_fluxes
578!=======================================================================
579
580END MODULE surf_ice
Note: See TracBrowser for help on using the repository browser.