source: trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 @ 2856

Last change on this file since 2856 was 2856, checked in by llange, 2 years ago

PEM
CO2 glacier flows added in a module
Maximum thickness of glacier before flow is not anymore hardcoded, but computed with the value of Tcond and the slope angle.
The method is from A.Grau Galfore (LPG), inspired from Nye et al., 2000
Works for high stress regim, will be adapted to consider also low stress regim (Tco2 > 200K);
LL
File size: 10.0 KB
Line 
1      module co2glaciers_mod
2        implicit none
3
4!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5!!!
6!!! Purpose: Compute CO2 glacier flows
7!!!
8!!! Author: LL
9!!!
10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11
12contains
13
14
15subroutine 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_slope,flag_co2flow,flag_co2flow_mesh)
16
17!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
18!!!
19!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
20!!!          the ice transfer
21!!!         
22!!!         
23!!! Author: LL
24!!!
25!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
26
27
28IMPLICIT NONE
29
30! arguments
31! ---------
32
33! Inputs:
34      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! # number of time sample, physical points, subslopes, index of the flat subslope
35      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
36      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
37      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
38      REAL,INTENT(IN) :: global_ave_ps_GCM          ! Global averaged surface pressure from the GCM [Pa]
39      REAL,INTENT(IN) :: global_ave_ps_PEM          ! global averaged surface pressure during the PEM iteration [Pa]
40     
41! Ouputs:
42      REAL,INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
43      REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
44      REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid)  ! same but within the mesh
45
46
47! Local
48      REAL :: Tcond(ngrid) !  Physical field: CO2 condensation temperature [K]
49      REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2  glacier before flow
50
51!-----------------------------
52      call computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
53
54      call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax)
55
56      call  transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
57   RETURN
58end subroutine
59
60
61
62subroutine compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax)
63
64     USE comconst_mod, ONLY: pi
65
66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
67!!!
68!!! Purpose: Compute the maximum thickness of CO2 glaciers given a slope angle
69!!!          before initating flow
70!!!         
71!!! Author: LL, based on theoretical work by A.Grau Galofre (LPG)
72!!!
73!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74
75   IMPLICIT NONE
76
77! arguments
78! --------
79
80! Inputs
81      INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
82      INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
83      REAL,INTENT(IN) :: Tcond(ngrid)     ! Physical field: CO2 condensation temperature [K]
84      REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
85! Outputs
86      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m]
87! Local
88      REAL, PARAMETER :: g = 3.71 ! surface gravity [m/s^2]
89      INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000
90      REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol]
91      REAL,PARAMETER :: Q = 59000. ! Activation Energy [J/mol], Nye et al., 2000
92      DOUBLE PRECISION,PARAMETER :: C = 1.8138e-21 ! Nye et al., 2000 [s/m^n]
93      DOUBLE PRECISION :: Ad = 1.202e11 ! Softness prefactor [MPa^-n] Nye et al., 2000
94      REAL :: Ro,Ho, S,ratio ! gemoetry from Nye et al., 2000
95      DOUBLE PRECISION :: A,Ao ! softness parameter [s/m^n]
96      DOUBLE PRECISION :: C1 ! intermediate variable
97      DOUBLE PRECISION :: t_0 ! relaxation time (assuming radial symetry) [s]
98      DOUBLE PRECISION :: u ! characteristic horizontal deformation rate [m/s]
99      DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the strest that an ice CO2 mass flowing under its weight balanced by viscosity
100      REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3]
101      INTEGER :: ig,islope ! loop variables
102      REAL :: slo_angle
103
104! 0. Geometry parameters
105      Ro = 200e3
106      Ho = 3000.
107      ratio = 2./3.
108      S = Ho/Ro*1/((2+1./n)*(1+1./n))*(1-ratio**(1+1./n))**(1./(2+1./n)-1)*ratio**(1+1./n-1)
109! 1. Flow parameters
110     Ao = 3**(1./(2*n+2))*Ad
111     do ig = 1,n
112     Ao = Ao*1e-6
113     enddo
114! 2. Compute rho_co2
115    DO ig = 1,ngrid
116      rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3  ! Mangan et al. 2017
117    ENDDO
118! 3. Compute max thickness
119    DO ig = 1,ngrid
120       A = Ao*exp(-Q/(Rg*Tcond(ig)))
121       C1 = A*(rho_co2(ig)*g)**(float(n))/float(n+2)
122       t_0 = Ro/(C1*(5*n+3))*(float(2*n+1)/float(n+1))**n
123       u = Ro/t_0
124       tau_d = (u*Ho*(rho_co2(ig)*g)**2*S**2/(2*A))**(1./(n+2))     
125       DO islope = 1,nslope
126          if(islope.eq.iflat) then
127            hmax(ig,islope) = 1.e6
128          else
129            slo_angle = abs(def_slope_mean(islope)*pi/180.)
130            hmax(ig,islope) = tau_d/(rho_co2(ig)*g*slo_angle)
131          endif
132       ENDDO
133    ENDDO
134RETURN
135
136end subroutine
137
138subroutine transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140!!!
141!!! Purpose: Transfer the excess of ice from one subslope to another
142!!!          No transfer between mesh at the time
143!!! Author: LL
144!!!
145!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
146
147      USE comconst_mod, ONLY: pi
148
149
150implicit none
151
152! arguments
153! --------
154
155! Inputs
156      INTEGER, INTENT(IN) :: ngrid,nslope !# of physical points and subslope
157      INTEGER, INTENT(IN) :: iflat ! index of the flat subslope
158      REAL, INTENT(IN) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh
159      REAL, INTENT(IN) :: def_slope_mean(nslope)       ! values of the subgrid slopes
160      REAL, INTENT(IN) :: hmax(ngrid,nslope)           ! maximum height of the CO2 glaciers before initiating flow [m]
161      REAL, INTENT(IN) :: Tcond(ngrid)          ! CO2 condensation temperature [K]
162! Outputs
163      REAL, INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! CO2 in the subslope [kg/m^2]
164      REAL, INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
165      REAL, INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! boolean to check if there is flow in the mesh
166! Local
167      INTEGER ig,islope ! loop
168      REAL rho_co2(ngrid) ! density of CO2, temperature dependant [kg/m^3]
169      INTEGER iaval ! ice will be transfered here
170
171
172
173! 0. Compute rho_co2
174    DO ig = 1,ngrid
175      rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3  ! Mangan et al. 2017
176    ENDDO
177
178! 1. Compute the transfer of ice
179
180       DO ig = 1,ngrid
181        DO islope = 1,nslope
182          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
183! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
184            IF(co2ice_slope(ig,islope).ge.rho_co2(ig)*hmax(ig,islope) * &
185                  cos(pi*def_slope_mean(islope)/180.)) THEN
186! Second: determine the flatest slopes possible:
187                IF(islope.gt.iflat) THEN
188                  iaval=islope-1
189                ELSE
190                 iaval=islope+1
191                ENDIF
192                do while ((iaval.ne.iflat).and.  &
193                    (subslope_dist(ig,iaval).eq.0))
194                  IF(iaval.gt.iflat) THEN
195                     iaval=iaval-1
196                  ELSE
197                     iaval=iaval+1
198                  ENDIF
199                enddo
200              co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
201               (co2ice_slope(ig,islope) - rho_co2(ig)*hmax(ig,islope) *     &
202               cos(pi*def_slope_mean(islope)/180.)) *             &
203               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
204               cos(pi*def_slope_mean(iaval)/180.) /               &
205               cos(pi*def_slope_mean(islope)/180.)
206
207              co2ice_slope(ig,islope)=rho_co2(ig)*hmax(ig,islope) *        &
208               cos(pi*def_slope_mean(islope)/180.)
209
210              flag_co2flow(ig,islope) = 1.
211              flag_co2flow_mesh(ig) = 1.
212            ENDIF ! co2ice > hmax
213          ENDIF ! iflat
214        ENDDO !islope
215       ENDDO !ig
216RETURN
217end subroutine
218
219     subroutine computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
220!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
221!!!
222!!! Purpose: Compute CO2 condensation temperature
223!!!
224!!! Author: LL
225!!!
226!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
227implicit none
228! arguments:
229! ----------
230
231! INPUT
232      INTEGER,INTENT(IN) :: timelen, ngrid ! # of timesample,physical points, subslopes
233      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
234      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
235      REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
236      REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
237! OUTPUT
238      REAL,INTENT(OUT) :: Tcond(ngrid) ! Physical points : condensation temperature of CO2, yearly averaged
239
240! LOCAL
241
242      INTEGER :: ig,it ! for loop
243      REAL :: ave ! intermediate to compute average
244      REAL :: alpha_clap, beta_clap ! Clapeyron law for CO2
245      alpha_clap = 23.3494 ! James et al. 1992
246      beta_clap = 3182.48  ! James et al. 1992
247
248!!!!!!!!!!!!!!!!!!!!!!!!!!!!
249
250
251      DO ig = 1,ngrid
252        ave = 0
253        DO it = 1,timelen
254           ave = ave + beta_clap/(alpha_clap-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100))
255        ENDDO
256        Tcond(ig) = ave/timelen
257       ENDDO
258RETURN
259
260
261      end subroutine
262      end module
Note: See TracBrowser for help on using the repository browser.