Context Navigation

surf_landice_mod.F90 @ 5969

Last change on this file since 5969 was 5969, checked in by yann meurdesoif, 3 weeks ago

Bug Fix + gpu port of surf_landice

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: 28.7 KB

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

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

Download in other formats:

Original Format