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

Last change on this file since 2918 was 2918, checked in by llange, 20 months ago

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