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

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

MARS PEM

  • Fix a bug in conf_pem ("Flux_geo" was expected in the .def while it should be 'fluxgeo', corrected)
  • Fix a bug in glaciers modules (the length for 'name_ice' was too large)
  • Fix a bug when writing tsoil in the PEM (ngrid x nsoil_PEM x nslope was given, while it was expecting a ngrid x nsoil_GCM x nslope)

PEM runs correctly now
LL
File size: 11.5 KB
RevLine 
[2995]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)
[3002]59   RETURN   
[2995]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      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
99
100   RETURN
101end subroutine
102
103
104
105subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
106
107     USE comconst_mod, ONLY: pi,g
108
109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
110!!!
111!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
112!!!          before initating flow
113!!!         
114!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
115!!!
116!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117
118   IMPLICIT NONE
119
120! arguments
121! --------
122
123! Inputs
124      INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
125      INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
126      REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
127      REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
[3002]128      character(len=3), INTENT(IN) :: name_ice ! Nature of the ice
[2995]129! Outputs
130      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
131! Local
132      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
133      REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
134      INTEGER :: ig,islope ! loop variables
135      REAL :: slo_angle
136
137! 1. Compute rho
138    if(name_ice.eq."co2") then
139      DO ig = 1,ngrid
140        DO islope = 1,nslope
141        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
142        tau_d = 5.e3
143        ENDDO
144      ENDDO
145    elseif (name_ice.eq."h2o") then
146      DO ig = 1,ngrid
147        DO islope = 1,nslope
148          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
149          tau_d = 1.e5
150        ENDDO
151      ENDDO
152    else
153      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
154    endif
155
156! 3. Compute max thickness
157    DO ig = 1,ngrid
158       DO islope = 1,nslope
159          if(islope.eq.iflat) then
160            hmax(ig,islope) = 1.e8
161          else
162            slo_angle = abs(def_slope_mean(islope)*pi/180.)
163            hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
164          endif
165       ENDDO
166    ENDDO
167RETURN
168
169end subroutine
170
171
172
173
174subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
175!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176!!!
177!!! Purpose: Transfer the excess of ice from one subslope to another
178!!!          No transfer between mesh at the time
179!!! Author: LL
180!!!
181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182
183      USE comconst_mod, ONLY: pi
184
185
186implicit none
187
188! arguments
189! --------
190
191! Inputs
192      INTEGER, INTENT(IN) :: ngrid,nslope               !# of physical points and subslope
193      INTEGER, INTENT(IN) :: iflat                      ! index of the flat subslope
194      REAL, INTENT(IN) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
195      REAL, INTENT(IN) :: def_slope_mean(nslope)        ! values of the subgrid slopes
196      REAL, INTENT(IN) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
197      REAL, INTENT(IN) :: Tice(ngrid,nslope)            ! Ice temperature[K]
[3002]198      character(len=3), INTENT(IN) :: name_ice              ! Nature of the ice
[2995]199
200! Outputs
201      REAL, INTENT(INOUT) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
202      REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
203      REAL, INTENT(INOUT) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
204! Local
205      INTEGER ig,islope ! loop
206      REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
207      INTEGER iaval ! ice will be transfered here
208
209! 0. Compute rho
210    if(name_ice.eq."co2") then
211      DO ig = 1,ngrid
212        DO islope = 1,nslope
213        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
214        ENDDO
215      ENDDO
216    elseif (name_ice.eq."h2o") then
217      DO ig = 1,ngrid
218        DO islope = 1,nslope
219          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
220        ENDDO
221      ENDDO
222    else
223      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
224    endif
225
226! 1. Compute the transfer of ice
227
228       DO ig = 1,ngrid
229        DO islope = 1,nslope
230          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
231! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
232            IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
233                  cos(pi*def_slope_mean(islope)/180.)) THEN
234! Second: determine the flatest slopes possible:
235                IF(islope.gt.iflat) THEN
236                  iaval=islope-1
237                ELSE
238                 iaval=islope+1
239                ENDIF
240                do while ((iaval.ne.iflat).and.  &
241                    (subslope_dist(ig,iaval).eq.0))
242                  IF(iaval.gt.iflat) THEN
243                     iaval=iaval-1
244                  ELSE
245                     iaval=iaval+1
246                  ENDIF
247                enddo
248              qice(ig,iaval) = qice(ig,iaval) + &
249               (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
250               cos(pi*def_slope_mean(islope)/180.)) *             &
251               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
252               cos(pi*def_slope_mean(iaval)/180.) /               &
253               cos(pi*def_slope_mean(islope)/180.)
254
255              qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
256               cos(pi*def_slope_mean(islope)/180.)
257
258              flag_flow(ig,islope) = 1.
259              flag_flowmesh(ig) = 1.
260            ENDIF ! co2ice > hmax
261          ENDIF ! iflat
262        ENDDO !islope
263       ENDDO !ig
264RETURN
265end subroutine
266
267     subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
268!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
269!!!
270!!! Purpose: Compute CO2 condensation temperature
271!!!
272!!! Author: LL
273!!!
274!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275 
276use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
277
278implicit none
279
280! arguments:
281! ----------
282
283! INPUT
284      INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
285      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
286      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
287      REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
288      REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
289! OUTPUT
290      REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged
291
292! LOCAL
293
294      INTEGER :: ig,it,islope ! for loop
295      REAL :: ave ! intermediate to compute average
296
297!!!!!!!!!!!!!!!!!!!!!!!!!!!!
298
299
300      DO ig = 1,ngrid
301        ave = 0
302        DO it = 1,timelen
303           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))
304        ENDDO
305        DO islope = 1,nslope
306          Tcond(ig,islope) = ave/timelen
307        ENDDO
308      ENDDO
309RETURN
310
311
312end subroutine
313end module
Note: See TracBrowser for help on using the repository browser.