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

Last change on this file since 3984 was 3983, checked in by jbclement, 8 days ago

PEM:

  • Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
  • Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
  • Ice reservoirs representation in the PEM is modernized.

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