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

Last change on this file since 3552 was 3532, checked in by jbclement, 3 weeks ago

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