source: trunk/LMDZ.COMMON/libf/evolution/glaciers.F90 @ 4002

Last change on this file since 4002 was 3991, checked in by jbclement, 5 weeks 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: 12.0 KB
Line 
1MODULE glaciers
2!-----------------------------------------------------------------------
3! NAME
4!     glaciers
5!
6! DESCRIPTION
7!     Compute flow and transfer of CO2 and H2O ice glaciers on slopes
8!     based on maximum ice thickness and basal drag conditions.
9!
10! AUTHORS & DATE
11!     L. Lange
12!     JB Clement, 2023-2025
13!
14! NOTES
15!
16!-----------------------------------------------------------------------
17
18! DECLARATION
19! -----------
20implicit none
21
22! PARAMETERS
23! ----------
24logical :: h2oice_flow ! True by default, flag to compute H2O ice flow
25logical :: co2ice_flow ! True by default, flag to compute CO2 ice flow
26
27contains
28!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29
30!=======================================================================
31SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow)
32!-----------------------------------------------------------------------
33! NAME
34!     flow_co2glaciers
35!
36! DESCRIPTION
37!     Compute maximum thickness and transfer CO2 ice between subslopes.
38!
39! AUTHORS & DATE
40!     L. Lange
41!     JB Clement, 2023-2025
42!
43! NOTES
44!
45!-----------------------------------------------------------------------
46
47! DECLARATION
48! -----------
49implicit none
50
51! ARGUMENTS
52! ---------
53integer,                                  intent(in)    :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
54real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
55real, dimension(ngrid),                   intent(in)    :: def_slope_mean                ! Grid points: values of the sub grid slope angles
56real, dimension(ngrid,timelen),           intent(in)    :: vmr_co2_PEM                   ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
57real, dimension(ngrid,timelen),           intent(in)    :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
58real,                                     intent(in)    :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
59real,                                     intent(in)    :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
60real, dimension(ngrid,nslope),            intent(inout) :: co2ice                        ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
61integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_co2flow                  ! Flag to see if there is flow on the subgrid slopes
62
63! LOCAL VARIABLES
64! ---------------
65real, dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K]
66real, dimension(ngrid,nslope) :: hmax  ! Maximum thickness before flow
67
68! CODE
69! ----
70write(*,*) "> Flow of CO2 glaciers"
71call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
72call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
73call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow)
74
75END SUBROUTINE flow_co2glaciers
76!=======================================================================
77
78!=======================================================================
79SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow)
80!-----------------------------------------------------------------------
81! NAME
82!     flow_h2oglaciers
83!
84! DESCRIPTION
85!     Compute maximum thickness and transfer H2O ice between subslopes.
86!
87! AUTHORS & DATE
88!     L. Lange
89!     JB Clement, 2023-2025
90!
91! NOTES
92!
93!-----------------------------------------------------------------------
94
95! DECLARATION
96! -----------
97implicit none
98
99! ARGUMENTS
100! ---------
101integer,                                  intent(in)    ::  ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
102real, dimension(ngrid,nslope),            intent(in)    ::  subslope_dist        ! Distribution of the subgrid slopes
103real, dimension(ngrid),                   intent(in)    ::  def_slope_mean       ! Values of the sub grid slope angles
104real, dimension(ngrid,nslope),            intent(in)    ::  Tice                 ! Ice Temperature [K]
105real, dimension(ngrid,nslope),            intent(inout) ::  h2oice               ! H2O ice on the subgrid slopes [kg/m^2]
106integer(kind=1), dimension(ngrid,nslope), intent(out)   ::  flag_h2oflow         ! Flow flag on subgrid slopes
107
108! LOCAL VARIABLES
109! ---------------
110real, dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow
111
112! CODE
113! ----
114write(*,*) "> Flow of H2O glaciers"
115call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
116call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow)
117
118END SUBROUTINE flow_h2oglaciers
119!=======================================================================
120
121!=======================================================================
122SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
123!-----------------------------------------------------------------------
124! NAME
125!     compute_hmaxglaciers
126!
127! DESCRIPTION
128!     Compute the maximum thickness of CO2 and H2O glaciers given a
129!     slope angle before initiating flow.
130!
131! AUTHORS & DATE
132!     L. Lange
133!     JB Clement, 2023-2025
134!
135! NOTES
136!     Based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
137!
138!-----------------------------------------------------------------------
139
140! DEPENDENCIES
141! ------------
142use ice_table,      only: rho_ice
143use stop_pem,       only: stop_clean
144use phys_constants, only: pi, g
145
146! DECLARATION
147! -----------
148implicit none
149
150! ARGUMENTS
151! ---------
152integer,                       intent(in)  :: ngrid, nslope  ! # of grid points and subslopes
153integer,                       intent(in)  :: iflat          ! index of the flat subslope
154real, dimension(nslope),       intent(in)  :: def_slope_mean ! Values of the subgrid slope angles [deg]
155real, dimension(ngrid,nslope), intent(in)  :: Tice           ! Ice temperature [K]
156character(3),                  intent(in)  :: name_ice       ! Nature of ice
157real, dimension(ngrid,nslope), intent(out) :: hmax           ! Maximum thickness before flow [m]
158
159! LOCAL VARIABLES
160! ---------------
161real    :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
162integer :: ig, islope
163real    :: slo_angle
164
165! CODE
166! ----
167select case (trim(adjustl(name_ice)))
168    case('h2o')
169        tau_d = 1.e5
170    case('co2')
171        tau_d = 5.e3
172    case default
173        call stop_clean("compute_hmaxglaciers","Type of ice not known!",1)
174end select
175
176do ig = 1,ngrid
177    do islope = 1,nslope
178        if (islope == iflat) then
179            hmax(ig,islope) = 1.e8
180        else
181            slo_angle = abs(def_slope_mean(islope)*pi/180.)
182            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
183        endif
184    enddo
185enddo
186
187END SUBROUTINE compute_hmaxglaciers
188!=======================================================================
189
190!=======================================================================
191SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow)
192!-----------------------------------------------------------------------
193! NAME
194!     transfer_ice_duringflow
195!
196! DESCRIPTION
197!     Transfer the excess of ice from one subslope to another.
198!     No transfer between mesh at the time.
199!
200! AUTHORS & DATE
201!     L. Lange
202!     JB Clement, 2023-2025
203!
204! NOTES
205!
206!-----------------------------------------------------------------------
207
208! DEPENDENCIES
209! ------------
210use ice_table,      only: rho_ice
211use stop_pem,       only: stop_clean
212use phys_constants, only: pi
213
214! DECLARATION
215! -----------
216implicit none
217
218! ARGUMENTS
219! ---------
220integer,                                  intent(in)    :: ngrid, nslope  ! # of physical points and subslope
221integer,                                  intent(in)    :: iflat          ! Index of the flat subslope
222real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist  ! Distribution of the subgrid slopes
223real, dimension(nslope),                  intent(in)    :: def_slope_mean ! Values of the subgrid slopes
224real, dimension(ngrid,nslope),            intent(in)    :: hmax           ! Maximum height before initiating flow [m]
225real, dimension(ngrid,nslope),            intent(in)    :: Tice           ! Ice temperature [K]
226real, dimension(ngrid,nslope),            intent(inout) :: qice           ! Ice in the subslope [kg/m^2]
227integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_flow      ! Flow flag on subgrid slopes
228
229! LOCAL VARIABLES
230! ---------------
231integer :: ig, islope
232integer :: iaval ! index where ice is transferred
233
234! CODE
235! ----
236flag_flow = 0
237
238do ig = 1,ngrid
239    do islope = 1,nslope
240        if (islope /= iflat) then ! ice can be infinite on flat ground
241! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
242            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
243! Second: determine the flatest slopes possible:
244                if (islope > iflat) then
245                    iaval = islope - 1
246                else
247                    iaval = islope + 1
248                endif
249                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
250                    if (iaval > iflat) then
251                        iaval = iaval - 1
252                    else
253                        iaval = iaval + 1
254                    endif
255                enddo
256                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) &
257                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
258
259                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
260                flag_flow(ig,islope) = 1
261            endif ! co2ice > hmax
262        endif ! iflat
263    enddo !islope
264enddo !ig
265
266END SUBROUTINE
267
268!=======================================================================
269SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
270!-----------------------------------------------------------------------
271! NAME
272!     computeTcondCO2
273!
274! DESCRIPTION
275!     Compute CO2 condensation temperature.
276!
277! AUTHORS & DATE
278!     L. Lange
279!     JB Clement, 2023-2025
280!
281! NOTES
282!
283!-----------------------------------------------------------------------
284
285! DEPENDENCIES
286! ------------
287use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2
288
289! DECLARATION
290! -----------
291implicit none
292
293! ARGUMENTS
294! ---------
295integer,                        intent(in)  :: timelen, ngrid, nslope
296real, dimension(ngrid,timelen), intent(in)  :: vmr_co2_PEM       ! VMR of CO2 in the first layer [mol/mol]
297real, dimension(ngrid,timelen), intent(in)  :: ps_PCM            ! Surface pressure in the PCM [Pa]
298real,                           intent(in)  :: ps_avg_global_ini ! Global averaged surfacepressure in the PCM [Pa]
299real,                           intent(in)  :: ps_avg_global     ! Global averaged surface pressure computed during the PEM iteration
300real, dimension(ngrid,nslope),  intent(out) :: Tcond             ! Condensation temperature of CO2, yearly averaged
301
302! LOCAL VARIABLES
303! ---------------
304integer :: ig, it
305
306! CODE
307! ----
308do ig = 1,ngrid
309    Tcond(ig,:) = 0
310    do it = 1,timelen
311        Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_global_ini/ps_avg_global/100))
312    enddo
313    Tcond(ig,:) = Tcond(ig,:)/timelen
314enddo
315
316END SUBROUTINE computeTcondCO2
317!=======================================================================
318
319END MODULE glaciers
Note: See TracBrowser for help on using the repository browser.