1 | MODULE surf_landice_mod |
---|
2 | |
---|
3 | IMPLICIT NONE |
---|
4 | |
---|
5 | CONTAINS |
---|
6 | |
---|
7 | !**************************************************************************************** |
---|
8 | |
---|
9 | SUBROUTINE surf_landice(itime, dtime, knon, knindex, & |
---|
10 | rlon, rlat, debut, lafin, & |
---|
11 | rmu0, lwdownm, albedo, pphi1, & |
---|
12 | swnet, lwnet, tsurf, p1lay, & |
---|
13 | cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, & |
---|
14 | AcoefH, AcoefQ, BcoefH, BcoefQ, & |
---|
15 | AcoefU, AcoefV, BcoefU, BcoefV, & |
---|
16 | AcoefQBS, BcoefQBS, & |
---|
17 | ps, u1, v1, gustiness, rugoro, pctsrf, & |
---|
18 | snow, qsurf, qsol, qbs1, agesno, & |
---|
19 | tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, & |
---|
20 | tsurf_new, dflux_s, dflux_l, & |
---|
21 | alt, slope, cloudf, & |
---|
22 | snowhgt, qsnow, to_ice, sissnow, & |
---|
23 | alb3, runoff, & |
---|
24 | flux_u1, flux_v1 & |
---|
25 | #ifdef ISO |
---|
26 | ,xtprecip_rain, xtprecip_snow,xtspechum,Rland_ice & |
---|
27 | ,xtsnow,xtsol,xtevap & |
---|
28 | #endif |
---|
29 | ) |
---|
30 | |
---|
31 | USE dimphy |
---|
32 | USE geometry_mod, ONLY: longitude, latitude |
---|
33 | USE surface_data, ONLY: type_ocean, calice, calsno, landice_opt, iflag_albcalc |
---|
34 | USE fonte_neige_mod, ONLY: fonte_neige, run_off_lic, fqcalving_global, ffonte_global, fqfonte_global, runofflic_global |
---|
35 | USE cpl_mod, ONLY: cpl_send_landice_fields |
---|
36 | USE calcul_fluxs_mod |
---|
37 | USE phys_local_var_mod, ONLY: zxrhoslic, zxustartlic, zxqsaltlic |
---|
38 | USE phys_output_var_mod, ONLY: snow_o, zfra_o |
---|
39 | #ifdef ISO |
---|
40 | USE fonte_neige_mod, ONLY: xtrun_off_lic |
---|
41 | USE infotrac_phy, ONLY: ntiso,niso |
---|
42 | USE isotopes_routines_mod, ONLY: calcul_iso_surf_lic_vectall |
---|
43 | #ifdef ISOVERIF |
---|
44 | USE isotopes_mod, ONLY: iso_eau,ridicule |
---|
45 | USE isotopes_verif_mod |
---|
46 | #endif |
---|
47 | #endif |
---|
48 | |
---|
49 | !FC |
---|
50 | USE ioipsl_getin_p_mod, ONLY: getin_p |
---|
51 | USE lmdz_blowing_snow_ini, ONLY: c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs |
---|
52 | USE lmdz_blowing_snow_ini, ONLY: rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs |
---|
53 | USE surf_inlandsis_mod, ONLY: surf_inlandsis |
---|
54 | USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INLANDSIS |
---|
55 | |
---|
56 | USE indice_sol_mod |
---|
57 | |
---|
58 | ! INCLUDE "indicesol.h" |
---|
59 | INCLUDE "dimsoil.h" |
---|
60 | INCLUDE "YOMCST.h" |
---|
61 | INCLUDE "clesphys.h" |
---|
62 | |
---|
63 | ! Input variables |
---|
64 | !**************************************************************************************** |
---|
65 | INTEGER, INTENT(IN) :: itime, knon |
---|
66 | INTEGER, DIMENSION(klon), INTENT(in) :: knindex |
---|
67 | REAL, INTENT(in) :: dtime |
---|
68 | REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiance |
---|
69 | REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiance |
---|
70 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf |
---|
71 | REAL, DIMENSION(klon), INTENT(IN) :: p1lay |
---|
72 | REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm |
---|
73 | REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs |
---|
74 | REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum |
---|
75 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ |
---|
76 | REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ |
---|
77 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV |
---|
78 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS |
---|
79 | REAL, DIMENSION(klon), INTENT(IN) :: ps |
---|
80 | REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1 |
---|
81 | REAL, DIMENSION(klon), INTENT(IN) :: rugoro |
---|
82 | REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf |
---|
83 | #ifdef ISO |
---|
84 | REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow |
---|
85 | REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum |
---|
86 | #endif |
---|
87 | |
---|
88 | |
---|
89 | LOGICAL, INTENT(IN) :: debut !true if first step |
---|
90 | LOGICAL, INTENT(IN) :: lafin !true if last step |
---|
91 | REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat |
---|
92 | REAL, DIMENSION(klon), INTENT(IN) :: rmu0 |
---|
93 | REAL, DIMENSION(klon), INTENT(IN) :: lwdownm !ylwdown |
---|
94 | REAL, DIMENSION(klon), INTENT(IN) :: albedo !mean albedo |
---|
95 | REAL, DIMENSION(klon), INTENT(IN) :: pphi1 |
---|
96 | REAL, DIMENSION(klon), INTENT(IN) :: alt !mean altitude of the grid box |
---|
97 | REAL, DIMENSION(klon), INTENT(IN) :: slope !mean slope in grid box |
---|
98 | REAL, DIMENSION(klon), INTENT(IN) :: cloudf !total cloud fraction |
---|
99 | |
---|
100 | ! In/Output variables |
---|
101 | !**************************************************************************************** |
---|
102 | REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol |
---|
103 | REAL, DIMENSION(klon), INTENT(INOUT) :: agesno |
---|
104 | REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil |
---|
105 | #ifdef ISO |
---|
106 | REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow, xtsol |
---|
107 | REAL, DIMENSION(niso,klon), INTENT(INOUT) :: Rland_ice |
---|
108 | #endif |
---|
109 | |
---|
110 | |
---|
111 | ! Output variables |
---|
112 | !**************************************************************************************** |
---|
113 | REAL, DIMENSION(klon), INTENT(OUT) :: qsurf |
---|
114 | REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h |
---|
115 | !albedo SB >>> |
---|
116 | ! REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! new albedo in visible SW interval |
---|
117 | ! REAL, DIMENSION(klon), INTENT(OUT) :: alb2 ! new albedo in near IR interval |
---|
118 | REAL, DIMENSION(6), INTENT(IN) :: SFRWL |
---|
119 | REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir, alb_dif |
---|
120 | !albedo SB <<< |
---|
121 | REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat |
---|
122 | REAL, DIMENSION(klon), INTENT(OUT) :: fluxbs |
---|
123 | REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new |
---|
124 | REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l |
---|
125 | REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 |
---|
126 | |
---|
127 | REAL, DIMENSION(klon), INTENT(OUT) :: alb3 |
---|
128 | REAL, DIMENSION(klon), INTENT(OUT) :: qsnow !column water in snow [kg/m2] |
---|
129 | REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt !Snow height (m) |
---|
130 | REAL, DIMENSION(klon), INTENT(OUT) :: to_ice |
---|
131 | REAL, DIMENSION(klon), INTENT(OUT) :: sissnow |
---|
132 | REAL, DIMENSION(klon), INTENT(OUT) :: runoff !Land ice runoff |
---|
133 | #ifdef ISO |
---|
134 | REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap |
---|
135 | ! real, DIMENSION(niso,klon) :: xtrun_off_lic_0_diag ! est une variable globale de |
---|
136 | ! fonte_neige |
---|
137 | #endif |
---|
138 | |
---|
139 | |
---|
140 | ! Local variables |
---|
141 | !**************************************************************************************** |
---|
142 | REAL, DIMENSION(klon) :: soilcap, soilflux |
---|
143 | REAL, DIMENSION(klon) :: cal, beta, dif_grnd |
---|
144 | REAL, DIMENSION(klon) :: zfra, alb_neig |
---|
145 | REAL, DIMENSION(klon) :: radsol |
---|
146 | REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar |
---|
147 | INTEGER :: i, j, nt |
---|
148 | REAL, DIMENSION(klon) :: fqfonte, ffonte |
---|
149 | REAL, DIMENSION(klon) :: run_off_lic_frac |
---|
150 | #ifdef ISO |
---|
151 | REAL, PARAMETER :: t_coup = 273.15 |
---|
152 | REAL, DIMENSION(klon) :: fqfonte_diag |
---|
153 | REAL, DIMENSION(klon) :: fq_fonte_diag |
---|
154 | REAL, DIMENSION(klon) :: snow_evap_diag |
---|
155 | REAL, DIMENSION(klon) :: fqcalving_diag |
---|
156 | REAL max_eau_sol_diag |
---|
157 | REAL, DIMENSION(klon) :: runoff_diag |
---|
158 | REAL, DIMENSION(klon) :: run_off_lic_diag |
---|
159 | REAL :: coeff_rel_diag |
---|
160 | INTEGER :: ixt |
---|
161 | REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec |
---|
162 | REAL, DIMENSION(klon) :: snow_prec,qsol_prec |
---|
163 | ! real, DIMENSION(klon) :: run_off_lic_0_diag |
---|
164 | #endif |
---|
165 | |
---|
166 | |
---|
167 | REAL, DIMENSION(klon) :: emis_new !Emissivity |
---|
168 | REAL, DIMENSION(klon) :: swdown, lwdown |
---|
169 | REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis) |
---|
170 | REAL, DIMENSION(klon) :: erod !erosion of surface snow (flux, kg/m2/s like evap) |
---|
171 | REAL, DIMENSION(klon) :: zsl_height, wind_velo !surface layer height, wind spd |
---|
172 | REAL, DIMENSION(klon) :: dens_air, snow_cont_air !air density; snow content air |
---|
173 | REAL, DIMENSION(klon) :: alb_soil !albedo of underlying ice |
---|
174 | REAL, DIMENSION(klon) :: pexner !Exner potential |
---|
175 | REAL :: pref |
---|
176 | REAL, DIMENSION(klon, nsoilmx) :: tsoil0 !modif |
---|
177 | REAL :: dtis ! subtimestep |
---|
178 | LOGICAL :: debut_is, lafin_is ! debut and lafin for inlandsis |
---|
179 | |
---|
180 | CHARACTER (len = 20) :: modname = 'surf_landice' |
---|
181 | CHARACTER (len = 80) :: abort_message |
---|
182 | |
---|
183 | REAL, DIMENSION(klon) :: alb1, alb2 |
---|
184 | REAL, DIMENSION(klon) :: precip_totsnow, evap_totsnow |
---|
185 | REAL, DIMENSION (klon, 6) :: alb6 |
---|
186 | REAL :: esalt |
---|
187 | REAL :: lambdasalt, fluxsalt, csalt, nunu, aa, bb, cc |
---|
188 | REAL :: tau_dens, maxerosion |
---|
189 | REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt |
---|
190 | REAL, DIMENSION(klon) :: fluxbs_1, fluxbs_2, bsweight_fresh |
---|
191 | LOGICAL, DIMENSION(klon) :: ok_remaining_freshsnow |
---|
192 | REAL :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd |
---|
193 | |
---|
194 | |
---|
195 | ! End definition |
---|
196 | !**************************************************************************************** |
---|
197 | !FC |
---|
198 | !FC |
---|
199 | REAL, SAVE :: alb_vis_sno_lic |
---|
200 | !$OMP THREADPRIVATE(alb_vis_sno_lic) |
---|
201 | REAL, SAVE :: alb_nir_sno_lic |
---|
202 | !$OMP THREADPRIVATE(alb_nir_sno_lic) |
---|
203 | LOGICAL, SAVE :: firstcall = .TRUE. |
---|
204 | !$OMP THREADPRIVATE(firstcall) |
---|
205 | |
---|
206 | |
---|
207 | !FC firtscall initializations |
---|
208 | !****************************************************************************************** |
---|
209 | #ifdef ISO |
---|
210 | #ifdef ISOVERIF |
---|
211 | ! write(*,*) 'surf_land_ice 1499' |
---|
212 | DO i=1,knon |
---|
213 | IF (iso_eau > 0) THEN |
---|
214 | CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), & |
---|
215 | 'surf_land_ice 126',errmax,errmaxrel) |
---|
216 | ENDIF !IF (iso_eau > 0) THEN |
---|
217 | ENDDO !DO i=1,knon |
---|
218 | #endif |
---|
219 | #endif |
---|
220 | |
---|
221 | IF (firstcall) THEN |
---|
222 | alb_vis_sno_lic = 0.77 |
---|
223 | CALL getin_p('alb_vis_sno_lic', alb_vis_sno_lic) |
---|
224 | PRINT*, 'alb_vis_sno_lic', alb_vis_sno_lic |
---|
225 | alb_nir_sno_lic = 0.77 |
---|
226 | CALL getin_p('alb_nir_sno_lic', alb_nir_sno_lic) |
---|
227 | PRINT*, 'alb_nir_sno_lic', alb_nir_sno_lic |
---|
228 | |
---|
229 | firstcall = .false. |
---|
230 | ENDIF |
---|
231 | !****************************************************************************************** |
---|
232 | |
---|
233 | ! Initialize output variables |
---|
234 | alb3(:) = 999999. |
---|
235 | alb2(:) = 999999. |
---|
236 | alb1(:) = 999999. |
---|
237 | fluxbs(:) = 0. |
---|
238 | runoff(:) = 0. |
---|
239 | !**************************************************************************************** |
---|
240 | ! Calculate total absorbed radiance at surface |
---|
241 | |
---|
242 | !**************************************************************************************** |
---|
243 | radsol(:) = 0.0 |
---|
244 | radsol(1:knon) = swnet(1:knon) + lwnet(1:knon) |
---|
245 | |
---|
246 | !**************************************************************************************** |
---|
247 | |
---|
248 | !**************************************************************************************** |
---|
249 | ! landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... |
---|
250 | ! landice_opt = 1 : prepare and call INterace Lmdz SISvat (INLANDSIS) |
---|
251 | !**************************************************************************************** |
---|
252 | |
---|
253 | IF (landice_opt == 1) THEN |
---|
254 | |
---|
255 | !**************************************************************************************** |
---|
256 | ! CALL to INLANDSIS interface |
---|
257 | !**************************************************************************************** |
---|
258 | IF (CPPKEY_INLANDSIS) THEN |
---|
259 | |
---|
260 | #ifdef ISO |
---|
261 | CALL abort_gcm('surf_landice 235','isotopes pas dans INLANDSIS',1) |
---|
262 | #endif |
---|
263 | |
---|
264 | debut_is = debut |
---|
265 | lafin_is = .false. |
---|
266 | ! Suppose zero surface speed |
---|
267 | u0(:) = 0.0 |
---|
268 | v0(:) = 0.0 |
---|
269 | |
---|
270 | CALL calcul_flux_wind(knon, dtime, & |
---|
271 | u0, v0, u1, v1, gustiness, cdragm, & |
---|
272 | AcoefU, AcoefV, BcoefU, BcoefV, & |
---|
273 | p1lay, temp_air, & |
---|
274 | flux_u1, flux_v1) |
---|
275 | |
---|
276 | |
---|
277 | ! Set constants and compute some input for SISVAT |
---|
278 | ! = 1000 hPa |
---|
279 | ! and calculate incoming flux for SW and LW interval: swdown, lwdown |
---|
280 | swdown(:) = 0.0 |
---|
281 | lwdown(:) = 0.0 |
---|
282 | snow_cont_air(:) = 0. ! the snow content in air is not a prognostic variable of the model |
---|
283 | alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set |
---|
284 | ustar(:) = 0. |
---|
285 | pref = 100000. |
---|
286 | DO i = 1, knon |
---|
287 | swdown(i) = swnet(i) / (1 - albedo(i)) |
---|
288 | lwdown(i) = lwdownm(i) |
---|
289 | wind_velo(i) = u1(i)**2 + v1(i)**2 |
---|
290 | wind_velo(i) = wind_velo(i)**0.5 |
---|
291 | pexner(i) = (p1lay(i) / pref)**(RD / RCPD) |
---|
292 | dens_air(i) = p1lay(i) / RD / temp_air(i) ! dry air density |
---|
293 | zsl_height(i) = pphi1(i) / RG |
---|
294 | tsoil0(i, :) = tsoil(i, :) |
---|
295 | ustar(i) = (cdragm(i) * (wind_velo(i)**2))**0.5 |
---|
296 | END DO |
---|
297 | |
---|
298 | dtis = dtime |
---|
299 | |
---|
300 | IF (lafin) THEN |
---|
301 | lafin_is = .true. |
---|
302 | END IF |
---|
303 | |
---|
304 | CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, & |
---|
305 | rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow, & |
---|
306 | zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & |
---|
307 | rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, & |
---|
308 | radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & |
---|
309 | AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & |
---|
310 | run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s, dflux_l, & |
---|
311 | tsurf_new, alb1, alb2, alb3, alb6, & |
---|
312 | emis_new, z0m, z0h, qsurf) |
---|
313 | |
---|
314 | debut_is = .false. |
---|
315 | |
---|
316 | |
---|
317 | ! Treatment of snow melting and calving |
---|
318 | |
---|
319 | ! for consistency with standard LMDZ, add calving to run_off_lic |
---|
320 | run_off_lic(:) = run_off_lic(:) + to_ice(:) |
---|
321 | |
---|
322 | DO i = 1, knon |
---|
323 | ffonte_global(knindex(i), is_lic) = ffonte(i) |
---|
324 | fqfonte_global(knindex(i), is_lic) = fqfonte(i)! net melting= melting - refreezing |
---|
325 | fqcalving_global(knindex(i), is_lic) = to_ice(i) ! flux |
---|
326 | runofflic_global(knindex(i)) = run_off_lic(i) |
---|
327 | ENDDO |
---|
328 | ! Here, we assume that the calving term is equal to the to_ice term |
---|
329 | ! (no ice accumulation) |
---|
330 | |
---|
331 | ELSE |
---|
332 | abort_message = 'Pb de coherence: landice_opt = 1 mais CPP_INLANDSIS = .false.' |
---|
333 | CALL abort_physic(modname, abort_message, 1) |
---|
334 | END IF |
---|
335 | |
---|
336 | ELSE |
---|
337 | |
---|
338 | !**************************************************************************************** |
---|
339 | ! Soil calculations |
---|
340 | |
---|
341 | !**************************************************************************************** |
---|
342 | |
---|
343 | ! EV: use calbeta |
---|
344 | CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd) |
---|
345 | |
---|
346 | |
---|
347 | ! use soil model and recalculate properly cal |
---|
348 | IF (soil_model) THEN |
---|
349 | CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, & |
---|
350 | longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux) |
---|
351 | cal(1:knon) = RCPD / soilcap(1:knon) |
---|
352 | radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) |
---|
353 | ELSE |
---|
354 | cal = RCPD * calice |
---|
355 | WHERE (snow > 0.0) cal = RCPD * calsno |
---|
356 | ENDIF |
---|
357 | |
---|
358 | |
---|
359 | !**************************************************************************************** |
---|
360 | ! Calulate fluxes |
---|
361 | |
---|
362 | !**************************************************************************************** |
---|
363 | ! beta(:) = 1.0 |
---|
364 | ! dif_grnd(:) = 0.0 |
---|
365 | |
---|
366 | ! Suppose zero surface speed |
---|
367 | u0(:) = 0.0 |
---|
368 | v0(:) = 0.0 |
---|
369 | u1_lay(:) = u1(:) - u0(:) |
---|
370 | v1_lay(:) = v1(:) - v0(:) |
---|
371 | |
---|
372 | CALL calcul_fluxs(knon, is_lic, dtime, & |
---|
373 | tsurf, p1lay, cal, beta, cdragh, cdragh, ps, & |
---|
374 | precip_rain, precip_snow, snow, qsurf, & |
---|
375 | radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
---|
376 | 1., AcoefH, AcoefQ, BcoefH, BcoefQ, & |
---|
377 | tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
---|
378 | |
---|
379 | #ifdef ISO |
---|
380 | #ifdef ISOVERIF |
---|
381 | !write(*,*) 'surf_land_ice 1499' |
---|
382 | DO i=1,knon |
---|
383 | IF (iso_eau > 0) THEN |
---|
384 | IF (snow(i) > ridicule) THEN |
---|
385 | CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), & |
---|
386 | 'surf_land_ice 1151',errmax,errmaxrel) |
---|
387 | ENDIF !IF ((snow(i) > ridicule)) THEN |
---|
388 | ENDIF !IF (iso_eau > 0) THEN |
---|
389 | ENDDO !DO i=1,knon |
---|
390 | #endif |
---|
391 | |
---|
392 | DO i=1,knon |
---|
393 | snow_prec(i)=snow(i) |
---|
394 | DO ixt=1,niso |
---|
395 | xtsnow_prec(ixt,i)=xtsnow(ixt,i) |
---|
396 | ENDDO !DO ixt=1,niso |
---|
397 | ! initialisation: |
---|
398 | fq_fonte_diag(i)=0.0 |
---|
399 | fqfonte_diag(i)=0.0 |
---|
400 | snow_evap_diag(i)=0.0 |
---|
401 | ENDDO !DO i=1,knon |
---|
402 | #endif |
---|
403 | |
---|
404 | CALL calcul_flux_wind(knon, dtime, & |
---|
405 | u0, v0, u1, v1, gustiness, cdragm, & |
---|
406 | AcoefU, AcoefV, BcoefU, BcoefV, & |
---|
407 | p1lay, temp_air, & |
---|
408 | flux_u1, flux_v1) |
---|
409 | |
---|
410 | |
---|
411 | !**************************************************************************************** |
---|
412 | ! Calculate albedo |
---|
413 | |
---|
414 | !**************************************************************************************** |
---|
415 | |
---|
416 | !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux" |
---|
417 | ! alb1(1 : knon) = 0.6 !IM cf FH/GK |
---|
418 | ! alb1(1 : knon) = 0.82 |
---|
419 | ! alb1(1 : knon) = 0.77 !211003 Ksta0.77 |
---|
420 | ! alb1(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5 |
---|
421 | !IM: KstaTER0.77 & LMD_ARMIP6 |
---|
422 | |
---|
423 | ! Attantion: alb1 and alb2 are not the same! |
---|
424 | alb1(1:knon) = alb_vis_sno_lic |
---|
425 | alb2(1:knon) = alb_nir_sno_lic |
---|
426 | |
---|
427 | |
---|
428 | !**************************************************************************************** |
---|
429 | ! Rugosity |
---|
430 | |
---|
431 | !**************************************************************************************** |
---|
432 | |
---|
433 | if (z0m_landice > 0.) then |
---|
434 | z0m(1:knon) = z0m_landice |
---|
435 | z0h(1:knon) = z0h_landice |
---|
436 | else |
---|
437 | ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2018 |
---|
438 | coefa = 0.1658 !0.1862 !Ant |
---|
439 | coefb = -50.3869 !-55.7718 !Ant |
---|
440 | ta1 = 253.15 !255. Ant |
---|
441 | ta2 = 273.15 |
---|
442 | ta3 = 273.15 + 3 |
---|
443 | z01 = exp(coefa * ta1 + coefb) !~0.2 ! ~0.25 mm |
---|
444 | z02 = exp(coefa * ta2 + coefb) !~6 !~7 mm |
---|
445 | z03 = z01 |
---|
446 | coefc = log(z03 / z02) / (ta3 - ta2) |
---|
447 | coefd = log(z03) - coefc * ta3 |
---|
448 | do j = 1, knon |
---|
449 | if (temp_air(j) < ta1) then |
---|
450 | z0m(j) = z01 |
---|
451 | else if (temp_air(j)>=ta1 .and. temp_air(j)<ta2) then |
---|
452 | z0m(j) = exp(coefa * temp_air(j) + coefb) |
---|
453 | else if (temp_air(j)>=ta2 .and. temp_air(j)<ta3) then |
---|
454 | ! if st > 0, melting induce smooth surface |
---|
455 | z0m(j) = exp(coefc * temp_air(j) + coefd) |
---|
456 | else |
---|
457 | z0m(j) = z03 |
---|
458 | endif |
---|
459 | z0h(j) = z0m(j) |
---|
460 | enddo |
---|
461 | |
---|
462 | endif |
---|
463 | |
---|
464 | |
---|
465 | !**************************************************************************************** |
---|
466 | ! Simple blowing snow param |
---|
467 | !**************************************************************************************** |
---|
468 | ! we proceed in 2 steps: |
---|
469 | ! first we erode - if possible -the accumulated snow during the time step |
---|
470 | ! then we update the density of the underlying layer and see if we can also erode |
---|
471 | ! this layer |
---|
472 | |
---|
473 | if (ok_bs) then |
---|
474 | fluxbs(:) = 0. |
---|
475 | do j = 1, knon |
---|
476 | ws1(j) = (u1(j)**2 + v1(j)**2)**0.5 |
---|
477 | ustar(j) = (cdragm(j) * (u1(j)**2 + v1(j)**2))**0.5 |
---|
478 | rhod(j) = p1lay(j) / RD / temp_air(j) |
---|
479 | ustart0(j) = (log(2.868) - log(1.625)) / 0.085 * sqrt(cdragm(j)) |
---|
480 | enddo |
---|
481 | |
---|
482 | ! 1st step: erosion of fresh snow accumulated during the time step |
---|
483 | do j = 1, knon |
---|
484 | if (precip_snow(j) > 0.) then |
---|
485 | rhos(j) = rhofresh_bs |
---|
486 | ! blowing snow flux formula used in MAR |
---|
487 | ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs)) |
---|
488 | ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs |
---|
489 | ! computation of qbs at the top of the saltation layer |
---|
490 | ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) |
---|
491 | esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j))) |
---|
492 | hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27) |
---|
493 | qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt |
---|
494 | ! calculation of erosion (flux positive towards the surface here) |
---|
495 | ! consistent with implicit resolution of turbulent mixing equation |
---|
496 | ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s |
---|
497 | ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) |
---|
498 | ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer |
---|
499 | ! (rho*qsalt*hsalt) |
---|
500 | ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step |
---|
501 | maxerosion = min(precip_snow(j), hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs) |
---|
502 | |
---|
503 | fluxbs_1(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) & |
---|
504 | / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime) |
---|
505 | fluxbs_1(j) = max(-maxerosion, fluxbs_1(j)) |
---|
506 | |
---|
507 | if (precip_snow(j) > abs(fluxbs_1(j))) then |
---|
508 | ok_remaining_freshsnow(j) = .true. |
---|
509 | bsweight_fresh(j) = 1. |
---|
510 | else |
---|
511 | ok_remaining_freshsnow(j) = .false. |
---|
512 | bsweight_fresh(j) = exp(-(abs(fluxbs_1(j)) - precip_snow(j)) / precip_snow(j)) |
---|
513 | endif |
---|
514 | else |
---|
515 | ok_remaining_freshsnow(j) = .false. |
---|
516 | fluxbs_1(j) = 0. |
---|
517 | bsweight_fresh(j) = 0. |
---|
518 | endif |
---|
519 | enddo |
---|
520 | |
---|
521 | |
---|
522 | ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step) |
---|
523 | ! this is done through the routine albsno |
---|
524 | CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:) + fluxbs_1(:)) |
---|
525 | |
---|
526 | ! 2nd step: |
---|
527 | ! computation of threshold friction velocity |
---|
528 | ! which depends on surface snow density |
---|
529 | do j = 1, knon |
---|
530 | if (ok_remaining_freshsnow(j)) then |
---|
531 | fluxbs_2(j) = 0. |
---|
532 | else |
---|
533 | ! we start eroding the underlying layer |
---|
534 | ! estimation of snow density |
---|
535 | ! snow density increases with snow age and |
---|
536 | ! increases even faster in case of sedimentation of blowing snow or rain |
---|
537 | tau_dens = max(tau_densmin_bs, tau_dens0_bs * exp(-abs(precip_bs(j)) / pbst_bs - & |
---|
538 | abs(precip_rain(j)) / prt_bs) * exp(-max(tsurf(j) - RTT, 0.))) |
---|
539 | rhos(j) = rhofresh_bs + (rhohard_bs - rhofresh_bs) * (1. - exp(-agesno(j) * 86400.0 / tau_dens)) |
---|
540 | ! blowing snow flux formula used in MAR |
---|
541 | ustart(j) = ustart0(j) * exp(max(rhoice_bs / rhofresh_bs - rhoice_bs / rhos(j), 0.)) * exp(max(0., rhos(j) - rhohard_bs)) |
---|
542 | ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs |
---|
543 | ! computation of qbs at the top of the saltation layer |
---|
544 | ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) |
---|
545 | esalt = 1. / (c_esalt_bs * max(1.e-6, ustar(j))) |
---|
546 | hsalt(j) = 0.08436 * (max(1.e-6, ustar(j))**1.27) |
---|
547 | qsalt(j) = (max(ustar(j)**2 - ustart(j)**2, 0.)) / (RG * hsalt(j)) * esalt |
---|
548 | ! calculation of erosion (flux positive towards the surface here) |
---|
549 | ! consistent with implicit resolution of turbulent mixing equation |
---|
550 | ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s |
---|
551 | ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) |
---|
552 | ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer |
---|
553 | ! (rho*qsalt*hsalt) |
---|
554 | maxerosion = hsalt(j) * qsalt(j) * rhod(j) / tau_eqsalt_bs |
---|
555 | fluxbs_2(j) = rhod(j) * ws1(j) * cdragh(j) * zeta_bs * (AcoefQBS(j) - qsalt(j)) & |
---|
556 | / (1. - rhod(j) * ws1(j) * cdragh(j) * zeta_bs * BcoefQBS(j) * dtime) |
---|
557 | fluxbs_2(j) = max(-maxerosion, fluxbs_2(j)) |
---|
558 | endif |
---|
559 | enddo |
---|
560 | |
---|
561 | |
---|
562 | |
---|
563 | |
---|
564 | ! final flux and outputs |
---|
565 | do j = 1, knon |
---|
566 | ! total flux is the erosion of fresh snow + |
---|
567 | ! a fraction of the underlying snow (if all the fresh snow has been eroded) |
---|
568 | ! the calculation of the fraction is quite delicate since we do not know |
---|
569 | ! how much time was needed to erode the fresh snow. We assume that this time |
---|
570 | ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh |
---|
571 | |
---|
572 | fluxbs(j) = fluxbs_1(j) + fluxbs_2(j) * (1. - bsweight_fresh(j)) |
---|
573 | i = knindex(j) |
---|
574 | zxustartlic(i) = ustart(j) |
---|
575 | zxrhoslic(i) = rhos(j) |
---|
576 | zxqsaltlic(i) = qsalt(j) |
---|
577 | enddo |
---|
578 | |
---|
579 | else ! not ok_bs |
---|
580 | ! those lines are useful to calculate the snow age |
---|
581 | CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) |
---|
582 | |
---|
583 | endif ! if ok_bs |
---|
584 | |
---|
585 | |
---|
586 | |
---|
587 | !**************************************************************************************** |
---|
588 | ! Calculate snow amount |
---|
589 | |
---|
590 | !**************************************************************************************** |
---|
591 | IF (ok_bs) THEN |
---|
592 | precip_totsnow(:) = precip_snow(:) + precip_bs(:) |
---|
593 | evap_totsnow(:) = evap(:) - fluxbs(:) ! flux bs is positive towards the surface (snow erosion) |
---|
594 | ELSE |
---|
595 | precip_totsnow(:) = precip_snow(:) |
---|
596 | evap_totsnow(:) = evap(:) |
---|
597 | ENDIF |
---|
598 | |
---|
599 | CALL fonte_neige(knon, is_lic, knindex, dtime, & |
---|
600 | tsurf, precip_rain, precip_totsnow, & |
---|
601 | snow, qsol, tsurf_new, evap_totsnow & |
---|
602 | #ifdef ISO |
---|
603 | ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag & |
---|
604 | ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag & |
---|
605 | #endif |
---|
606 | ) |
---|
607 | |
---|
608 | |
---|
609 | #ifdef ISO |
---|
610 | #ifdef ISOVERIF |
---|
611 | DO i=1,knon |
---|
612 | IF (iso_eau > 0) THEN |
---|
613 | CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, & |
---|
614 | 'surf_landice_mod 217',errmax,errmaxrel) |
---|
615 | ENDIF !IF (iso_eau > 0) THEN |
---|
616 | ENDDO !DO i=1,knon |
---|
617 | #endif |
---|
618 | |
---|
619 | CALL calcul_iso_surf_lic_vectall(klon,knon, & |
---|
620 | evap,snow_evap_diag,Tsurf_new,snow, & |
---|
621 | fq_fonte_diag,fqfonte_diag,dtime,t_coup, & |
---|
622 | precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, & |
---|
623 | xtspechum,spechum,ps,Rland_ice, & |
---|
624 | xtevap,xtsnow,fqcalving_diag, & |
---|
625 | knindex,is_lic,run_off_lic_diag,coeff_rel_diag & |
---|
626 | ) |
---|
627 | |
---|
628 | ! call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag) |
---|
629 | |
---|
630 | #endif |
---|
631 | |
---|
632 | WHERE (snow(1:knon) < 0.0001) agesno(1:knon) = 0. |
---|
633 | zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon) / (snow(1:knon) + 10.0))) |
---|
634 | |
---|
635 | END IF ! landice_opt |
---|
636 | |
---|
637 | |
---|
638 | !**************************************************************************************** |
---|
639 | ! Send run-off on land-ice to coupler if coupled ocean. |
---|
640 | ! run_off_lic has been calculated in fonte_neige or surf_inlandsis |
---|
641 | ! If landice_opt>=2, corresponding call is done from surf_land_orchidee |
---|
642 | !**************************************************************************************** |
---|
643 | IF (type_ocean=='couple' .AND. landice_opt < 2) THEN |
---|
644 | ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic)) |
---|
645 | run_off_lic_frac(:) = 0.0 |
---|
646 | DO j = 1, knon |
---|
647 | i = knindex(j) |
---|
648 | run_off_lic_frac(j) = pctsrf(i, is_lic) |
---|
649 | ENDDO |
---|
650 | |
---|
651 | CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac) |
---|
652 | ENDIF |
---|
653 | |
---|
654 | ! transfer runoff rate [kg/m2/s](!) to physiq for output |
---|
655 | runoff(1:knon) = run_off_lic(1:knon) / dtime |
---|
656 | |
---|
657 | snow_o = 0. |
---|
658 | zfra_o = 0. |
---|
659 | DO j = 1, knon |
---|
660 | i = knindex(j) |
---|
661 | snow_o(i) = snow(j) |
---|
662 | zfra_o(i) = zfra(j) |
---|
663 | ENDDO |
---|
664 | |
---|
665 | |
---|
666 | !albedo SB >>> |
---|
667 | select case(NSW) |
---|
668 | case(2) |
---|
669 | alb_dir(1:knon, 1) = alb1(1:knon) |
---|
670 | alb_dir(1:knon, 2) = alb2(1:knon) |
---|
671 | case(4) |
---|
672 | alb_dir(1:knon, 1) = alb1(1:knon) |
---|
673 | alb_dir(1:knon, 2) = alb2(1:knon) |
---|
674 | alb_dir(1:knon, 3) = alb2(1:knon) |
---|
675 | alb_dir(1:knon, 4) = alb2(1:knon) |
---|
676 | case(6) |
---|
677 | alb_dir(1:knon, 1) = alb1(1:knon) |
---|
678 | alb_dir(1:knon, 2) = alb1(1:knon) |
---|
679 | alb_dir(1:knon, 3) = alb1(1:knon) |
---|
680 | alb_dir(1:knon, 4) = alb2(1:knon) |
---|
681 | alb_dir(1:knon, 5) = alb2(1:knon) |
---|
682 | alb_dir(1:knon, 6) = alb2(1:knon) |
---|
683 | |
---|
684 | IF ((landice_opt == 1) .AND. (iflag_albcalc == 2)) THEN |
---|
685 | alb_dir(1:knon, 1) = alb6(1:knon, 1) |
---|
686 | alb_dir(1:knon, 2) = alb6(1:knon, 2) |
---|
687 | alb_dir(1:knon, 3) = alb6(1:knon, 3) |
---|
688 | alb_dir(1:knon, 4) = alb6(1:knon, 4) |
---|
689 | alb_dir(1:knon, 5) = alb6(1:knon, 5) |
---|
690 | alb_dir(1:knon, 6) = alb6(1:knon, 6) |
---|
691 | ENDIF |
---|
692 | |
---|
693 | end select |
---|
694 | alb_dif = alb_dir |
---|
695 | !albedo SB <<< |
---|
696 | |
---|
697 | END SUBROUTINE surf_landice |
---|
698 | |
---|
699 | !**************************************************************************************** |
---|
700 | |
---|
701 | END MODULE surf_landice_mod |
---|
702 | |
---|
703 | |
---|
704 | |
---|