source: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90 @ 3985

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

PEM:
Addition of a module "phys_constants" to read and store physical parameter of the planet properly, i.e. without going through the module "comcstfi_h" and/or "comconst_mod".
JBC
File size: 12.9 KB
Line 
1MODULE glaciers_mod
2
3implicit none
4
5! Different types of ice
6type :: ice
7    real :: h2o
8    real :: co2
9end type ice
10
11! Perennial ice manage by the PEM
12type(ice), dimension(:,:), allocatable :: perice
13
14! Flags for ice flow
15logical :: h2oice_flow  ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
16logical :: co2ice_flow  ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
17
18! Threshold to consider H2O ice as watercap
19real :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
20
21! Ice properties
22real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
23real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
24
25!=======================================================================
26contains
27!=======================================================================
28
29SUBROUTINE set_perice4PCM(ngrid,nslope,PCMfrost,is_perh2oice,PCMh2oice,PCMco2ice)
30
31use metamorphism,   only: iPCM_h2ofrost
32use comslope_mod,   only: subslope_dist, def_slope_mean
33use phys_constants, only: pi
34
35
36implicit none
37
38! Arguments
39!----------
40integer, intent(in) :: ngrid, nslope
41real, dimension(:,:,:), intent(inout) :: PCMfrost
42logical, dimension(ngrid),        intent(out) :: is_perh2oice
43real,    dimension(ngrid,nslope), intent(out) :: PCMco2ice, PCMh2oice
44
45! Local variables
46!----------------
47integer :: i
48
49! Code
50!-----
51write(*,*) '> Reconstructing perennial ic for the PCM'
52PCMco2ice(:,:) = perice(:,:)%co2
53PCMh2oice(:,:) = 0. ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
54do i = 1,ngrid
55    ! Is H2O ice still considered as an infinite reservoir for the PCM?
56    if (sum(perice(i,:)%h2o*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
57        ! There is enough ice to be considered as an infinite reservoir
58        is_perh2oice(i) = .true.
59    else
60        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
61        is_perh2oice(i) = .false.
62        PCMfrost(i,iPCM_h2ofrost,:) = PCMfrost(i,iPCM_h2ofrost,:) + perice(i,:)%h2o
63        perice(i,:)%h2o = 0.
64    endif
65enddo
66
67END SUBROUTINE set_perice4PCM
68!=======================================================================
69
70SUBROUTINE ini_ice(ngrid,nslope)
71
72implicit none
73
74! Arguments
75!----------
76integer, intent(in) :: ngrid, nslope
77
78! Local variables
79!----------------
80
81! Code
82!-----
83if (.not. allocated(perice)) allocate(perice(ngrid,nslope))
84
85END SUBROUTINE ini_ice
86!=======================================================================
87
88SUBROUTINE end_ice()
89
90implicit none
91
92! Arguments
93!----------
94
95! Local variables
96!----------------
97
98! Code
99!-----
100if (allocated(perice)) deallocate(perice)
101
102END SUBROUTINE end_ice
103
104!=======================================================================
105
106SUBROUTINE 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)
107
108!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
109!!!
110!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
111!!!          the ice transfer
112!!!
113!!!
114!!! Author: LL
115!!!
116!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117
118implicit none
119
120! Inputs
121!-------
122integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
123real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
124real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Grid points: values of the sub grid slope angles
125real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Grid points x Time field : VMR of co2 in the first layer [mol/mol]
126real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
127real,                              intent(in) :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
128real,                              intent(in) :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
129! Ouputs
130!-------
131real, dimension(ngrid,nslope), intent(inout) :: co2ice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
132integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_co2flow ! Flag to see if there is flow on the subgrid slopes
133! Local
134!------
135real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
136real, dimension(ngrid,nslope) :: hmax  ! Grid points x Slope field: maximum thickness for co2  glacier before flow
137
138write(*,*) "> Flow of CO2 glaciers"
139call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
140call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
141call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow)
142
143END SUBROUTINE flow_co2glaciers
144
145!=======================================================================
146
147SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow)
148
149!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150!!!
151!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
152!!!          the ice transfer
153!!!
154!!!
155!!! Author: LL
156!!!
157!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158
159implicit none
160
161! arguments
162! ---------
163
164! Inputs
165! ------
166integer,                       intent(in) :: ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
167real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Grid points x Slopes : Distribution of the subgrid slopes
168real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
169real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
170! Outputs
171!--------
172real, dimension(ngrid,nslope), intent(inout) :: h2oice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
173integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_h2oflow ! Flag to see if there is flow on the subgrid slopes
174! Local
175! -----
176real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2  glacier before flow
177
178write(*,*) "> Flow of H2O glaciers"
179call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
180call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow)
181
182END SUBROUTINE flow_h2oglaciers
183
184!=======================================================================
185
186SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
187
188use ice_table_mod,  only: rho_ice
189use abort_pem_mod,  only: abort_pem
190use phys_constants, only: pi, g
191
192!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193!!!
194!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
195!!!
196!!! Author: LL, based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
197!!!
198!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199
200implicit none
201
202! Inputs
203! ------
204integer,                       intent(in) :: ngrid, nslope  ! # of grid points and subslopes
205integer,                       intent(in) :: iflat          ! index of the flat subslope
206real, dimension(nslope),       intent(in) :: def_slope_mean ! Slope field: Values of the subgrid slope angles [deg]
207real, dimension(ngrid,nslope), intent(in) :: Tice           ! Physical field:  ice temperature [K]
208character(3),                  intent(in) :: name_ice       ! Nature of ice
209! Outputs
210! -------
211real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum  thickness before flaw [m]
212! Local
213! -----
214real    :: 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
215integer :: ig, islope ! loop variables
216real    :: slo_angle
217
218select case (trim(adjustl(name_ice)))
219    case('h2o')
220        tau_d = 1.e5
221    case('co2')
222        tau_d = 5.e3
223    case default
224        call abort_pem("compute_hmaxglaciers","Type of ice not known!",1)
225end select
226
227do ig = 1,ngrid
228    do islope = 1,nslope
229        if (islope == iflat) then
230            hmax(ig,islope) = 1.e8
231        else
232            slo_angle = abs(def_slope_mean(islope)*pi/180.)
233            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
234        endif
235    enddo
236enddo
237
238END SUBROUTINE compute_hmaxglaciers
239
240!=======================================================================
241
242SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow)
243!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
244!!!
245!!! Purpose: Transfer the excess of ice from one subslope to another
246!!!          No transfer between mesh at the time
247!!! Author: LL
248!!!
249!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
250
251use ice_table_mod,  only: rho_ice
252use abort_pem_mod,  only: abort_pem
253use phys_constants, only: pi
254
255implicit none
256
257! Inputs
258!-------
259integer,                       intent(in) :: ngrid, nslope  ! # of physical points and subslope
260integer,                       intent(in) :: iflat          ! index of the flat subslope
261real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Distribution of the subgrid slopes within the mesh
262real, dimension(nslope),       intent(in) :: def_slope_mean ! values of the subgrid slopes
263real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
264real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
265! Outputs
266!--------
267real, dimension(ngrid,nslope), intent(inout) :: qice ! CO2 in the subslope [kg/m^2]
268integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_flow ! Flag to see if there is flow on the subgrid slopes
269! Local
270!------
271integer :: ig, islope ! loop
272integer :: iaval      ! ice will be transfered here
273
274flag_flow = 0
275
276do ig = 1,ngrid
277    do islope = 1,nslope
278        if (islope /= iflat) then ! ice can be infinite on flat ground
279! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
280            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
281! Second: determine the flatest slopes possible:
282                if (islope > iflat) then
283                    iaval=islope-1
284                else
285                    iaval = islope + 1
286                endif
287                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
288                    if (iaval > iflat) then
289                        iaval = iaval - 1
290                    else
291                        iaval = iaval + 1
292                    endif
293                enddo
294                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.)) &
295                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
296
297                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
298                flag_flow(ig,islope) = 1
299            endif ! co2ice > hmax
300        endif ! iflat
301    enddo !islope
302enddo !ig
303
304END SUBROUTINE
305
306!=======================================================================
307
308SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310!!!
311!!! Purpose: Compute CO2 condensation temperature
312!!!
313!!! Author: LL
314!!!
315!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
316
317use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2
318
319implicit none
320
321! arguments:
322! ----------
323
324! INPUT
325integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
326real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Grid points x times field: VMR of CO2 in the first layer [mol/mol]
327real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Grid points x times field: surface pressure in the PCM [Pa]
328real,                           intent(in) :: ps_avg_global_ini      ! Global averaged surfacepressure in the PCM [Pa]
329real,                           intent(in) :: ps_avg_global          ! Global averaged surface pressure computed during the PEM iteration
330
331! OUTPUT
332real, dimension(ngrid,nslope), intent(out) :: Tcond ! Grid points: condensation temperature of CO2, yearly averaged
333
334! LOCAL
335integer :: ig, it ! For loop
336
337do ig = 1,ngrid
338    Tcond(ig,:) = 0
339    do it = 1,timelen
340        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))
341    enddo
342    Tcond(ig,:) = Tcond(ig,:)/timelen
343enddo
344
345END SUBROUTINE computeTcondCO2
346
347END MODULE glaciers_mod
Note: See TracBrowser for help on using the repository browser.