Context Navigation

surf_landice_mod.F90 @ 5136

Last change on this file since 5136 was 5128, checked in by abarral, 8 weeks ago
Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments. Merge the three different versions of abort_gcm into one Fix seq, para 3D compilation broken from r5107 onwards (lint) usual + Remove uneeded fixed-form continuations
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format