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

Last change on this file since 2995 was 2995, checked in by llange, 16 months ago

Mars PEM
*Implementation of the H2O glacier flow laws
*The algorithm for glacier flow is now more generic and not specific for co2 ice
*Principle: if ice thickness > ice mass, computed from models (cf note attached in the wiki), then the excess of ice is transfered
LL
File size: 11.5 KB
Line 
1module glaciers_mod
2       
3 implicit none
4 LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
5 LOGICAL  h2oglaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
6!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7!!!
8!!! Purpose: Compute CO2 glacier flows
9!!!
10!!! Author: LL
11!!!
12!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13
14contains
15
16
17subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
18
19!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20!!!
21!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
22!!!          the ice transfer
23!!!         
24!!!         
25!!! Author: LL
26!!!
27!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
28
29
30IMPLICIT NONE
31
32! arguments
33! ---------
34
35! Inputs:
36      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
37      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
38      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
39      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
40      REAL,INTENT(IN) :: global_ave_ps_GCM          ! Global averaged surface pressure from the GCM [Pa]
41      REAL,INTENT(IN) :: global_ave_ps_PEM          ! global averaged surface pressure during the PEM iteration [Pa]
42     
43! Ouputs:
44      REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
45      REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
46      REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid)  ! same but within the mesh
47
48
49! Local
50      REAL :: Tcond(ngrid,nslope) !  Physical field: CO2 condensation temperature [K]
51      REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
52
53!-----------------------------
54      call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
55
56      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
57
58      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
59   RETURN
60end subroutine
61
62
63
64
65
66subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
67
68!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69!!!
70!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
71!!!          the ice transfer
72!!!         
73!!!         
74!!! Author: LL
75!!!
76!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77
78
79IMPLICIT NONE
80
81! arguments
82! ---------
83
84! Inputs:
85      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
86      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
87      REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K]
88! Ouputs:
89      REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
90      REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
91      REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
92! Local
93      REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
94
95!-----------------------------
96
97      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
98
99      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
100
101   RETURN
102end subroutine
103
104
105
106subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
107
108     USE comconst_mod, ONLY: pi,g
109
110!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111!!!
112!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
113!!!          before initating flow
114!!!         
115!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
116!!!
117!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118
119   IMPLICIT NONE
120
121! arguments
122! --------
123
124! Inputs
125      INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
126      INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
127      REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
128      REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
129      character(len=30), INTENT(IN) :: name_ice ! Nature of the ice
130! Outputs
131      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
132! Local
133      DOUBLE PRECISION :: 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
134      REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
135      INTEGER :: ig,islope ! loop variables
136      REAL :: slo_angle
137
138
139! 1. Compute rho
140    if(name_ice.eq."co2") then
141      DO ig = 1,ngrid
142        DO islope = 1,nslope
143        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
144        tau_d = 5.e3
145        ENDDO
146      ENDDO
147    elseif (name_ice.eq."h2o") then
148      DO ig = 1,ngrid
149        DO islope = 1,nslope
150          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
151          tau_d = 1.e5
152        ENDDO
153      ENDDO
154    else
155      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
156    endif
157
158! 3. Compute max thickness
159    DO ig = 1,ngrid
160       DO islope = 1,nslope
161          if(islope.eq.iflat) then
162            hmax(ig,islope) = 1.e8
163          else
164            slo_angle = abs(def_slope_mean(islope)*pi/180.)
165            hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
166          endif
167       ENDDO
168    ENDDO
169RETURN
170
171end subroutine
172
173
174
175
176subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
177!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178!!!
179!!! Purpose: Transfer the excess of ice from one subslope to another
180!!!          No transfer between mesh at the time
181!!! Author: LL
182!!!
183!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184
185      USE comconst_mod, ONLY: pi
186
187
188implicit none
189
190! arguments
191! --------
192
193! Inputs
194      INTEGER, INTENT(IN) :: ngrid,nslope               !# of physical points and subslope
195      INTEGER, INTENT(IN) :: iflat                      ! index of the flat subslope
196      REAL, INTENT(IN) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
197      REAL, INTENT(IN) :: def_slope_mean(nslope)        ! values of the subgrid slopes
198      REAL, INTENT(IN) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
199      REAL, INTENT(IN) :: Tice(ngrid,nslope)            ! Ice temperature[K]
200      character(len=30), INTENT(IN) :: name_ice              ! Nature of the ice
201
202! Outputs
203      REAL, INTENT(INOUT) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
204      REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
205      REAL, INTENT(INOUT) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
206! Local
207      INTEGER ig,islope ! loop
208      REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
209      INTEGER iaval ! ice will be transfered here
210
211
212
213! 0. Compute rho
214    if(name_ice.eq."co2") then
215      DO ig = 1,ngrid
216        DO islope = 1,nslope
217        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
218        ENDDO
219      ENDDO
220    elseif (name_ice.eq."h2o") then
221      DO ig = 1,ngrid
222        DO islope = 1,nslope
223          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
224        ENDDO
225      ENDDO
226    else
227      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
228    endif
229
230! 1. Compute the transfer of ice
231
232       DO ig = 1,ngrid
233        DO islope = 1,nslope
234          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
235! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
236            IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
237                  cos(pi*def_slope_mean(islope)/180.)) THEN
238! Second: determine the flatest slopes possible:
239                IF(islope.gt.iflat) THEN
240                  iaval=islope-1
241                ELSE
242                 iaval=islope+1
243                ENDIF
244                do while ((iaval.ne.iflat).and.  &
245                    (subslope_dist(ig,iaval).eq.0))
246                  IF(iaval.gt.iflat) THEN
247                     iaval=iaval-1
248                  ELSE
249                     iaval=iaval+1
250                  ENDIF
251                enddo
252              qice(ig,iaval) = qice(ig,iaval) + &
253               (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
254               cos(pi*def_slope_mean(islope)/180.)) *             &
255               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
256               cos(pi*def_slope_mean(iaval)/180.) /               &
257               cos(pi*def_slope_mean(islope)/180.)
258
259              qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
260               cos(pi*def_slope_mean(islope)/180.)
261
262              flag_flow(ig,islope) = 1.
263              flag_flowmesh(ig) = 1.
264            ENDIF ! co2ice > hmax
265          ENDIF ! iflat
266        ENDDO !islope
267       ENDDO !ig
268RETURN
269end subroutine
270
271     subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
272!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273!!!
274!!! Purpose: Compute CO2 condensation temperature
275!!!
276!!! Author: LL
277!!!
278!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
279 
280use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
281
282implicit none
283
284! arguments:
285! ----------
286
287! INPUT
288      INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
289      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
290      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
291      REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
292      REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
293! OUTPUT
294      REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged
295
296! LOCAL
297
298      INTEGER :: ig,it,islope ! for loop
299      REAL :: ave ! intermediate to compute average
300
301!!!!!!!!!!!!!!!!!!!!!!!!!!!!
302
303
304      DO ig = 1,ngrid
305        ave = 0
306        DO it = 1,timelen
307           ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100))
308        ENDDO
309        DO islope = 1,nslope
310          Tcond(ig,islope) = ave/timelen
311        ENDDO
312      ENDDO
313RETURN
314
315
316end subroutine
317end module
Note: See TracBrowser for help on using the repository browser.