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

Last change on this file since 3969 was 3961, checked in by jbclement, 3 months ago

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