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

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

PEM:
Big cleaning/improvements of the PEM:

  • Conversion of "abort_pem.F" and "soil_settings_PEM.F" into Fortran 90;
  • Transformation of every PEM subroutines into module;
  • Rewriting of many subroutines with modern Fortran syntax;
  • Correction of a bug in "pem.F90" when calling 'recomp_tend_co2_slope'. The arguments were given in disorder and emissivity was missing;
  • Update of "launch_pem.sh" in deftank.

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