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

Last change on this file since 3082 was 3082, checked in by jbclement, 14 months ago

PEM:

  • Correction of a bug in the initialization of constants. The correct modules are now used: 'comcstfi_h' (and no longer 'comconst_mod'!) in the general case and 'comcstfi_mod' in the case of generic model;
  • Addition of the variable 'ecritpem' in "run_PEM.def" to set the frequency of outputs in the "diagfi.nc". By default, 'ecritpem = 1' which means there is one output at each PEM year.

JBC
File size: 11.6 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
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,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
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(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,nslope) !  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 computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
53
54      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
55
56      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
57   RETURN   
58end subroutine
59
60
61
62
63
64subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
65
66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
67!!!
68!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
69!!!          the ice transfer
70!!!         
71!!!         
72!!! Author: LL
73!!!
74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75
76
77IMPLICIT NONE
78
79! arguments
80! ---------
81
82! Inputs:
83      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
84      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
85      REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K]
86! Ouputs:
87      REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
88      REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
89      REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
90! Local
91      REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
92
93!-----------------------------
94
95      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
96      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
97
98   RETURN
99end subroutine
100
101
102subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
103
104use abort_pem_mod, only: abort_pem
105#ifndef CPP_STD
106    use comcstfi_h,   only: pi, g
107#else
108    use comcstfi_mod, only: pi, g
109#endif
110
111!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
112!!!
113!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
114!!!          before initating flow
115!!!         
116!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
117!!!
118!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119
120   IMPLICIT NONE
121
122! arguments
123! --------
124
125! Inputs
126      INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
127      INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
128      REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
129      REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
130      character(len=3), INTENT(IN) :: name_ice ! Nature of the ice
131! Outputs
132      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
133! Local
134      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
135      REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
136      INTEGER :: ig,islope ! loop variables
137      REAL :: slo_angle
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
185use abort_pem_mod, only: abort_pem
186#ifndef CPP_STD
187    use comcstfi_h,   only: pi
188#else
189    use comcstfi_mod, only: pi
190#endif
191
192implicit none
193
194! arguments
195! --------
196
197! Inputs
198      INTEGER, INTENT(IN) :: ngrid,nslope               !# of physical points and subslope
199      INTEGER, INTENT(IN) :: iflat                      ! index of the flat subslope
200      REAL, INTENT(IN) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
201      REAL, INTENT(IN) :: def_slope_mean(nslope)        ! values of the subgrid slopes
202      REAL, INTENT(IN) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
203      REAL, INTENT(IN) :: Tice(ngrid,nslope)            ! Ice temperature[K]
204      character(len=3), INTENT(IN) :: name_ice              ! Nature of the ice
205
206! Outputs
207      REAL, INTENT(INOUT) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
208      REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
209      REAL, INTENT(INOUT) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
210! Local
211      INTEGER ig,islope ! loop
212      REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
213      INTEGER iaval ! ice will be transfered here
214
215! 0. Compute rho
216    if(name_ice.eq."co2") then
217      DO ig = 1,ngrid
218        DO islope = 1,nslope
219        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
220        ENDDO
221      ENDDO
222    elseif (name_ice.eq."h2o") then
223      DO ig = 1,ngrid
224        DO islope = 1,nslope
225          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
226        ENDDO
227      ENDDO
228    else
229      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
230    endif
231
232! 1. Compute the transfer of ice
233
234       DO ig = 1,ngrid
235        DO islope = 1,nslope
236          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
237! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
238            IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
239                  cos(pi*def_slope_mean(islope)/180.)) THEN
240! Second: determine the flatest slopes possible:
241                IF(islope.gt.iflat) THEN
242                  iaval=islope-1
243                ELSE
244                 iaval=islope+1
245                ENDIF
246                do while ((iaval.ne.iflat).and.  &
247                    (subslope_dist(ig,iaval).eq.0))
248                  IF(iaval.gt.iflat) THEN
249                     iaval=iaval-1
250                  ELSE
251                     iaval=iaval+1
252                  ENDIF
253                enddo
254              qice(ig,iaval) = qice(ig,iaval) + &
255               (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
256               cos(pi*def_slope_mean(islope)/180.)) *             &
257               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
258               cos(pi*def_slope_mean(iaval)/180.) /               &
259               cos(pi*def_slope_mean(islope)/180.)
260
261              qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
262               cos(pi*def_slope_mean(islope)/180.)
263
264              flag_flow(ig,islope) = 1.
265              flag_flowmesh(ig) = 1.
266            ENDIF ! co2ice > hmax
267          ENDIF ! iflat
268        ENDDO !islope
269       ENDDO !ig
270RETURN
271end subroutine
272
273     subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
274!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275!!!
276!!! Purpose: Compute CO2 condensation temperature
277!!!
278!!! Author: LL
279!!!
280!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281 
282use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
283
284implicit none
285
286! arguments:
287! ----------
288
289! INPUT
290      INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
291      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
292      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
293      REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
294      REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
295! OUTPUT
296      REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged
297
298! LOCAL
299
300      INTEGER :: ig,it,islope ! for loop
301      REAL :: ave ! intermediate to compute average
302
303!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304
305
306      DO ig = 1,ngrid
307        ave = 0
308        DO it = 1,timelen
309           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))
310        ENDDO
311        DO islope = 1,nslope
312          Tcond(ig,islope) = ave/timelen
313        ENDDO
314      ENDDO
315RETURN
316
317
318end subroutine
319end module
Note: See TracBrowser for help on using the repository browser.