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

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

PEM:
Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.
JBC
File size: 10.8 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! DECLARATION
17! -----------
18implicit none
19
20! MODULE VARIABLES
21! ----------------
22real, dimension(:,:), allocatable :: h2o_ice ! H2O ice surface [kg.m-2]
23real, dimension(:,:), allocatable :: co2_ice ! CO2 ice surface [kg.m-2]
24real :: h2oice_cap_threshold                 ! Threshold to consider H2O ice as infinite reservoir [kg.m-2]
25
26! MODULE PARAMETERS
27! -----------------
28real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
29real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
30
31contains
32!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33
34!=======================================================================
35SUBROUTINE set_perice4PCM(ngrid,nslope,PCMfrost,is_h2o_perice,h2oice_PCM,co2ice_PCM)
36!-----------------------------------------------------------------------
37! NAME
38!     set_perice4PCM
39!
40! DESCRIPTION
41!     Reconstruct perennial ice for PCM from PEM ice computations.
42!
43! AUTHORS & DATE
44!     JB Clement, 12/2025
45!
46! NOTES
47!
48!-----------------------------------------------------------------------
49
50! DEPENDENCIES
51! ------------
52use metamorphism,   only: iPCM_h2ofrost
53use comslope_mod,   only: subslope_dist, def_slope_mean
54use phys_constants, only: pi
55
56! DECLARATION
57! -----------
58implicit none
59
60! ARGUMENTS
61! ---------
62integer,                       intent(in)    :: ngrid, nslope
63real, dimension(:,:,:),        intent(inout) :: PCMfrost               ! PCM frost
64logical, dimension(ngrid),     intent(out)   :: is_h2o_perice          ! H2O perennial ice flag
65real, dimension(ngrid,nslope), intent(out)   :: h2oice_PCM, co2ice_PCM ! Ice for PCM
66
67! LOCAL VARIABLES
68! ---------------
69integer :: i
70
71! CODE
72! ----
73write(*,*) '> Reconstructing perennial ic for the PCM'
74co2ice_PCM(:,:) = co2_ice(:,:)
75h2oice_PCM(:,:) = 0. ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
76do i = 1,ngrid
77    ! Is H2O ice still considered as an infinite reservoir for the PCM?
78    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180.)) > h2oice_cap_threshold) then
79        ! There is enough ice to be considered as an infinite reservoir
80        is_h2o_perice(i) = .true.
81    else
82        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
83        is_h2o_perice(i) = .false.
84        PCMfrost(i,iPCM_h2ofrost,:) = PCMfrost(i,iPCM_h2ofrost,:) + h2o_ice(i,:)
85        h2o_ice(i,:) = 0.
86    endif
87enddo
88
89END SUBROUTINE set_perice4PCM
90!=======================================================================
91
92!=======================================================================
93SUBROUTINE ini_surf_ice(ngrid,nslope)
94!-----------------------------------------------------------------------
95! NAME
96!     ini_surf_ice
97!
98! DESCRIPTION
99!     Initialize surface ice arrays.
100!
101! AUTHORS & DATE
102!     JB Clement 12/2025
103!
104! NOTES
105!
106!-----------------------------------------------------------------------
107
108! DECLARATION
109! -----------
110implicit none
111
112! ARGUMENTS
113! ---------
114integer, intent(in) :: ngrid, nslope
115
116! CODE
117! ----
118if (.not. allocated(h2o_ice)) allocate(h2o_ice(ngrid,nslope))
119if (.not. allocated(co2_ice)) allocate(co2_ice(ngrid,nslope))
120
121END SUBROUTINE ini_surf_ice
122!=======================================================================
123
124!=======================================================================
125SUBROUTINE end_surf_ice()
126!-----------------------------------------------------------------------
127! NAME
128!     end_surf_ice
129!
130! DESCRIPTION
131!     Deallocate surface ice arrays.
132!
133! AUTHORS & DATE
134!     JB Clement, 12/2025
135!
136! NOTES
137!
138!-----------------------------------------------------------------------
139
140! DECLARATION
141! -----------
142implicit none
143
144! CODE
145! ----
146if (allocated(h2o_ice)) deallocate(h2o_ice)
147if (allocated(co2_ice)) deallocate(co2_ice)
148
149END SUBROUTINE end_surf_ice
150!=======================================================================
151
152!=======================================================================
153SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
154!-----------------------------------------------------------------------
155! NAME
156!     evol_co2_ice
157!
158! DESCRIPTION
159!     Compute the evolution of CO2 ice over one year.
160!
161! AUTHORS & DATE
162!     R. Vandemeulebrouck
163!     L. Lange
164!     JB Clement, 2023-2025
165!
166! NOTES
167!
168!-----------------------------------------------------------------------
169
170! DEPENDENCIES
171! ------------
172use evolution, only: dt
173
174! DECLARATION
175! -----------
176implicit none
177
178! ARGUMENTS
179! ---------
180integer,                       intent(in)    :: ngrid, nslope
181real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! CO2 ice
182real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of CO2 ice over one year
183real, dimension(ngrid,nslope), intent(out)   :: zshift_surf ! Elevation shift for the surface [m]
184
185! LOCAL VARIABLES
186! ---------------
187real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
188
189! CODE
190! ----
191! Evolution of CO2 ice for each physical point
192write(*,*) '> Evolution of CO2 ice'
193
194co2_ice_old = co2_ice
195co2_ice = co2_ice + d_co2ice*dt
196where (co2_ice < 0.)
197    co2_ice = 0.
198    d_co2ice = -co2_ice_old/dt
199end where
200
201zshift_surf = d_co2ice*dt/rho_co2ice
202
203END SUBROUTINE evol_co2_ice
204!=======================================================================
205
206!=======================================================================
207SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit)
208!-----------------------------------------------------------------------
209! NAME
210!     evol_h2o_ice
211!
212! DESCRIPTION
213!     Compute the evolution of H2O ice over one year.
214!
215! AUTHORS & DATE
216!     R. Vandemeulebrouck
217!     L. Lange
218!     JB Clement, 2023-2025
219!
220! NOTES
221!
222!-----------------------------------------------------------------------
223
224! DEPENDENCIES
225! ------------
226use evolution,     only: dt
227use stopping_crit, only: stopping_crit_h2o, stopFlags
228
229! DECLARATION
230! -----------
231implicit none
232
233! ARGUMENTS
234! ---------
235integer,                       intent(in)    :: ngrid, nslope
236real, dimension(ngrid),        intent(in)    :: cell_area                ! Area of each mesh grid (m^2)
237real, dimension(ngrid),        intent(in)    :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in soil (kg/m^2)
238real, dimension(ngrid),        intent(in)    :: delta_h2o_icetablesublim ! Mass condensed/sublimated at ice table (kg/m^2)
239real, dimension(ngrid,nslope), intent(inout) :: h2o_ice                  ! H2O ice (kg/m^2)
240real, dimension(ngrid,nslope), intent(inout) :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
241type(stopFlags),               intent(inout) :: stopCrit                 ! Stopping criterion
242real, dimension(ngrid,nslope), intent(out)   :: zshift_surf              ! Elevation shift for the surface [m]
243
244! LOCAL VARIABLES
245! ---------------
246integer                       :: i, islope
247real                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! H2O balance variables
248real, dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
249
250! CODE
251! ----
252write(*,*) '> Evolution of H2O ice'
253
254if (ngrid == 1) then ! In 1D
255    h2o_ice = h2o_ice + d_h2oice*dt
256    zshift_surf = d_h2oice*dt/rho_h2oice
257else ! In 3D
258    call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
259    if (stopCrit%stop_code() > 0) return
260
261    call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
262    h2o_ice = h2o_ice + d_h2oice_new*dt
263    zshift_surf = d_h2oice_new*dt/rho_h2oice
264endif
265
266END SUBROUTINE evol_h2o_ice
267!=======================================================================
268
269!=======================================================================
270SUBROUTINE balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
271!-----------------------------------------------------------------------
272! NAME
273!     balance_h2oice_reservoirs
274!
275! DESCRIPTION
276!     Balance H2O flux from and into atmosphere across reservoirs.
277!
278! AUTHORS & DATE
279!     JB Clement, 2025
280!
281! NOTES
282!
283!-----------------------------------------------------------------------
284
285! DEPENDENCIES
286! ------------
287use evolution, only: dt
288
289! DECLARATION
290! -----------
291implicit none
292
293! ARGUMENTS
294! ---------
295integer,                       intent(in)    :: ngrid, nslope ! # of grid points, # of subslopes
296real,                          intent(in)    :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
297real, dimension(ngrid,nslope), intent(in)    :: h2o_ice       ! H2O ice (kg/m^2)
298real, dimension(ngrid,nslope), intent(inout) :: d_h2oice      ! Tendency of H2O ice (kg/m^2/year)
299real, dimension(ngrid,nslope), intent(out)   :: d_h2oice_new  ! Adjusted tendencies to keep balance between donor and recipient reservoirs
300
301! LOCAL VARIABLES
302! ---------------
303integer :: i, islope
304real    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target ! Balance variables
305
306! CODE
307! ----
308S_target = (S_atm_2_h2o + S_h2o_2_atm)/2.
309S_target_cond_h2oice = S_atm_2_h2oice + S_target - S_atm_2_h2o
310S_target_subl_h2oice = S_h2oice_2_atm + S_target - S_h2o_2_atm
311
312d_h2oice_new = 0.
313S_ghostice = 0.
314do i = 1,ngrid
315    do islope = 1,nslope
316        if (d_h2oice(i,islope) > 0.) then ! Condensing
317            d_h2oice_new(i,islope) = d_h2oice_new(i,islope)*S_target_cond_h2oice/S_atm_2_h2oice
318        else if (d_h2oice(i,islope) < 0.) then ! Sublimating
319            d_target = d_h2oice(i,islope)*S_target_subl_h2oice/S_h2oice_2_atm
320            if (abs(d_target*dt) < h2o_ice(i,islope)) then ! Enough ice to sublimate everything
321                d_h2oice_new(i,islope) = d_target
322            else ! Not enough ice to sublimate everything
323                ! We sublimate what we can
324                d_h2oice_new(i,islope) = h2o_ice(i,islope)/dt
325                ! It means the tendency is zero next time
326                d_h2oice(i,islope) = 0.
327                ! We compute the amount of H2O ice that we could not make sublimate
328                S_ghostice = S_ghostice + abs(d_target*dt) - h2o_ice(i,islope)
329            endif
330        endif
331    enddo
332enddo
333
334! We need to remove this ice unable to sublimate from places where ice condenseds in order to keep balance
335where (d_h2oice_new > 0.) d_h2oice_new = d_h2oice_new*(S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice
336
337END SUBROUTINE balance_h2oice_reservoirs
338!=======================================================================
339
340END MODULE surf_ice
Note: See TracBrowser for help on using the repository browser.