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

Last change on this file since 3321 was 3308, checked in by jbclement, 7 months ago

PEM:
Few small corrections to make the PEM work in 3D, in particular concerning the initialization of the planet type and the evolution of ice with slopes.
JBC
File size: 12.9 KB
RevLine 
[3149]1MODULE glaciers_mod
[2995]2
[3149]3implicit none
4
5
[3161]6! Flags for ice management
[3308]7logical :: h2oice_flow  ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
8logical :: co2ice_flow  ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
[3161]9logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def"
10logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def"
11
12! Thresholds for ice management
[3256]13real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
14real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2]
15real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2]
[3161]16
[3149]17!=======================================================================
[2995]18contains
[3149]19!=======================================================================
[2995]20
[3149]21SUBROUTINE 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)
[2995]22
23!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
24!!!
25!!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
26!!!          the ice transfer
27!!!         
28!!!         
29!!! Author: LL
30!!!
31!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32
[3149]33implicit none
[2995]34
35! arguments
36! ---------
37
38! Inputs:
[3149]39      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
40      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
41      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
42      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
43      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
44      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
45      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
[2995]46     
47! Ouputs:
[3149]48      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
49      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
50      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
[2995]51
52! Local
[3149]53      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
54      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
[2995]55
56!-----------------------------
[3149]57write(*,*) "Flow of CO2 glacier"
58   
59call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
60call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
61call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
[2995]62
[3149]63END SUBROUTINE flow_co2glaciers
[2995]64
[3149]65!=======================================================================
[2995]66
[3149]67SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
[2995]68
69!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
70!!!
71!!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
72!!!          the ice transfer
73!!!         
74!!!         
75!!! Author: LL
76!!!
77!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78
[3149]79implicit none
[2995]80
81! arguments
82! ---------
83
84! Inputs:
[3149]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]
[2995]88! Ouputs:
[3149]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
[2995]92! Local
[3149]93      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
[2995]94
95!-----------------------------
[3149]96write(*,*) "Flow of H2O glaciers"
[2995]97
[3149]98call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
99call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
[2995]100
[3149]101END SUBROUTINE flow_h2oglaciers
[2995]102
[3149]103!=======================================================================
[2995]104
[3149]105SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
[2995]106
[3076]107use abort_pem_mod, only: abort_pem
[3082]108#ifndef CPP_STD
109    use comcstfi_h,   only: pi, g
110#else
111    use comcstfi_mod, only: pi, g
112#endif
[2995]113
114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115!!!
116!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
117!!!          before initating flow
118!!!         
119!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
120!!!
121!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122
[3149]123implicit none
[2995]124
125! arguments
126! --------
127
128! Inputs
[3149]129      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
130      integer,intent(in) :: iflat        ! index of the flat subslope
131      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
132      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
133      character(len=3), intent(in) :: name_ice ! Nature of the ice
[2995]134! Outputs
[3149]135      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
[2995]136! Local
137      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
[3149]138      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
139      integer :: ig,islope ! loop variables
140      real :: slo_angle
[2995]141
142! 1. Compute rho
143    if(name_ice.eq."co2") then
144      DO ig = 1,ngrid
145        DO islope = 1,nslope
146        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
147        tau_d = 5.e3
148        ENDDO
149      ENDDO
150    elseif (name_ice.eq."h2o") then
151      DO ig = 1,ngrid
152        DO islope = 1,nslope
153          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
154          tau_d = 1.e5
155        ENDDO
156      ENDDO
157    else
158      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
159    endif
160
161! 3. Compute max thickness
162    DO ig = 1,ngrid
163       DO islope = 1,nslope
164          if(islope.eq.iflat) then
165            hmax(ig,islope) = 1.e8
166          else
167            slo_angle = abs(def_slope_mean(islope)*pi/180.)
168            hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
169          endif
170       ENDDO
171    ENDDO
[3149]172END SUBROUTINE compute_hmaxglaciers
[2995]173
[3149]174!=======================================================================
[2995]175
[3149]176SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
[2995]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
[3076]185use abort_pem_mod, only: abort_pem
[3082]186#ifndef CPP_STD
187    use comcstfi_h,   only: pi
188#else
189    use comcstfi_mod, only: pi
190#endif
[2995]191
192implicit none
193
194! arguments
195! --------
196
197! Inputs
[3149]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
[2995]205
206! Outputs
[3149]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
[2995]210! Local
[3149]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
[2995]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
[3149]270END SUBROUTINE
[2995]271
[3149]272!=======================================================================
273
274SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
[2995]275!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
276!!!
277!!! Purpose: Compute CO2 condensation temperature
278!!!
279!!! Author: LL
280!!!
281!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282 
283use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
284
285implicit none
286
287! arguments:
288! ----------
289
290! INPUT
[3149]291integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
292real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
293real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
294real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
295real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
296
[2995]297! OUTPUT
[3149]298real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
[2995]299
300! LOCAL
[3149]301integer :: ig, it ! For loop
302real    :: ave    ! Intermediate to compute average
[2995]303
304!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3149]305do ig = 1,ngrid
306    ave = 0
307    do it = 1,timelen
308        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))
309    enddo
310    Tcond(ig,:) = ave/timelen
311enddo
[2995]312
[3149]313END SUBROUTINE computeTcondCO2
[2995]314
[3149]315END MODULE glaciers_mod
Note: See TracBrowser for help on using the repository browser.