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

Last change on this file since 3149 was 3149, checked in by jbclement, 12 months ago

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC
File size: 12.3 KB
Line 
1MODULE glaciers_mod
2
3implicit none
4
5logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
6logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
7
8!=======================================================================
9contains
10!=======================================================================
11
12SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
13
14!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15!!!
16!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
17!!!          the ice transfer
18!!!         
19!!!         
20!!! Author: LL
21!!!
22!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24implicit none
25
26! arguments
27! ---------
28
29! Inputs:
30      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
31      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
32      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
33      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
34      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
35      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
36      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
37     
38! Ouputs:
39      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
40      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
41      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
42
43! Local
44      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
45      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
46
47!-----------------------------
48write(*,*) "Flow of CO2 glacier"
49   
50call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
51call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
52call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
53
54END SUBROUTINE flow_co2glaciers
55
56!=======================================================================
57
58SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
59
60!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
61!!!
62!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
63!!!          the ice transfer
64!!!         
65!!!         
66!!! Author: LL
67!!!
68!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69
70implicit none
71
72! arguments
73! ---------
74
75! Inputs:
76      integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
77      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
78      real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
79! Ouputs:
80      real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
81      real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
82      real,intent(inout) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
83! Local
84      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
85
86!-----------------------------
87write(*,*) "Flow of H2O glaciers"
88
89call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
90call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
91
92END SUBROUTINE flow_h2oglaciers
93
94!=======================================================================
95
96SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
97
98use abort_pem_mod, only: abort_pem
99#ifndef CPP_STD
100    use comcstfi_h,   only: pi, g
101#else
102    use comcstfi_mod, only: pi, g
103#endif
104
105!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106!!!
107!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
108!!!          before initating flow
109!!!         
110!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
111!!!
112!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
113
114implicit none
115
116! arguments
117! --------
118
119! Inputs
120      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
121      integer,intent(in) :: iflat        ! index of the flat subslope
122      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
123      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
124      character(len=3), intent(in) :: name_ice ! Nature of the ice
125! Outputs
126      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
127! Local
128      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
129      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
130      integer :: ig,islope ! loop variables
131      real :: slo_angle
132
133! 1. Compute rho
134    if(name_ice.eq."co2") then
135      DO ig = 1,ngrid
136        DO islope = 1,nslope
137        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
138        tau_d = 5.e3
139        ENDDO
140      ENDDO
141    elseif (name_ice.eq."h2o") then
142      DO ig = 1,ngrid
143        DO islope = 1,nslope
144          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
145          tau_d = 1.e5
146        ENDDO
147      ENDDO
148    else
149      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
150    endif
151
152! 3. Compute max thickness
153    DO ig = 1,ngrid
154       DO islope = 1,nslope
155          if(islope.eq.iflat) then
156            hmax(ig,islope) = 1.e8
157          else
158            slo_angle = abs(def_slope_mean(islope)*pi/180.)
159            hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
160          endif
161       ENDDO
162    ENDDO
163END SUBROUTINE compute_hmaxglaciers
164
165!=======================================================================
166
167SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
169!!!
170!!! Purpose: Transfer the excess of ice from one subslope to another
171!!!          No transfer between mesh at the time
172!!! Author: LL
173!!!
174!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175
176use abort_pem_mod, only: abort_pem
177#ifndef CPP_STD
178    use comcstfi_h,   only: pi
179#else
180    use comcstfi_mod, only: pi
181#endif
182
183implicit none
184
185! arguments
186! --------
187
188! Inputs
189      integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
190      integer, intent(in) :: iflat                      ! index of the flat subslope
191      real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
192      real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
193      real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
194      real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
195      character(len=3), intent(in) :: name_ice              ! Nature of the ice
196
197! Outputs
198      real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
199      real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
200      real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
201! Local
202      integer ig,islope ! loop
203      real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
204      integer iaval ! ice will be transfered here
205
206! 0. Compute rho
207    if(name_ice.eq."co2") then
208      DO ig = 1,ngrid
209        DO islope = 1,nslope
210        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
211        ENDDO
212      ENDDO
213    elseif (name_ice.eq."h2o") then
214      DO ig = 1,ngrid
215        DO islope = 1,nslope
216          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
217        ENDDO
218      ENDDO
219    else
220      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
221    endif
222
223! 1. Compute the transfer of ice
224
225       DO ig = 1,ngrid
226        DO islope = 1,nslope
227          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
228! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
229            IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
230                  cos(pi*def_slope_mean(islope)/180.)) THEN
231! Second: determine the flatest slopes possible:
232                IF(islope.gt.iflat) THEN
233                  iaval=islope-1
234                ELSE
235                 iaval=islope+1
236                ENDIF
237                do while ((iaval.ne.iflat).and.  &
238                    (subslope_dist(ig,iaval).eq.0))
239                  IF(iaval.gt.iflat) THEN
240                     iaval=iaval-1
241                  ELSE
242                     iaval=iaval+1
243                  ENDIF
244                enddo
245              qice(ig,iaval) = qice(ig,iaval) + &
246               (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
247               cos(pi*def_slope_mean(islope)/180.)) *             &
248               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
249               cos(pi*def_slope_mean(iaval)/180.) /               &
250               cos(pi*def_slope_mean(islope)/180.)
251
252              qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
253               cos(pi*def_slope_mean(islope)/180.)
254
255              flag_flow(ig,islope) = 1.
256              flag_flowmesh(ig) = 1.
257            ENDIF ! co2ice > hmax
258          ENDIF ! iflat
259        ENDDO !islope
260       ENDDO !ig
261END SUBROUTINE
262
263!=======================================================================
264
265SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
266!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
267!!!
268!!! Purpose: Compute CO2 condensation temperature
269!!!
270!!! Author: LL
271!!!
272!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273 
274use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
275
276implicit none
277
278! arguments:
279! ----------
280
281! INPUT
282integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
283real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
284real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
285real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
286real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
287
288! OUTPUT
289real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
290
291! LOCAL
292integer :: ig, it ! For loop
293real    :: ave    ! Intermediate to compute average
294
295!!!!!!!!!!!!!!!!!!!!!!!!!!!!
296do ig = 1,ngrid
297    ave = 0
298    do it = 1,timelen
299        ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100))
300    enddo
301    Tcond(ig,:) = ave/timelen
302enddo
303
304END SUBROUTINE computeTcondCO2
305
306END MODULE glaciers_mod
Note: See TracBrowser for help on using the repository browser.