Context Navigation

surf_landice_mod.F90 @ 4835

Last change on this file since 4835 was 4835, checked in by evignon, 3 months ago
commission pour la suite du travail sur la mise a jour de la param de neige soufflee
File size: 25.9 KB

Rev	Line
[3927]	1	!
	2	MODULE surf_landice_mod
	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, &
[4523]	14	cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
[3927]	15	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	16	AcoefU, AcoefV, BcoefU, BcoefV, &
[4529]	17	AcoefQBS, BcoefQBS, &
[3927]	18	ps, u1, v1, gustiness, rugoro, pctsrf, &
[4523]	19	snow, qsurf, qsol, qbs1, agesno, &
	20	tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
[3927]	21	tsurf_new, dflux_s, dflux_l, &
[3940]	22	alt, slope, cloudf, &
[3927]	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
[3940]	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
[3927]	35	USE cpl_mod, ONLY : cpl_send_landice_fields
	36	USE calcul_fluxs_mod
[4835]	37	USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
[4523]	38	USE phys_output_var_mod, ONLY : snow_o,zfra_o
[3927]	39	#ifdef ISO
	40	USE fonte_neige_mod, ONLY : xtrun_off_lic
[4143]	41	USE infotrac_phy, ONLY : ntiso,niso
[3927]	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
[3975]	48	USE geometry_mod, ONLY : longitude,latitude
	49
[3927]	50	!FC
	51	USE ioipsl_getin_p_mod, ONLY : getin_p
[4835]	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
[3940]	54	#ifdef CPP_INLANDSIS
	55	USE surf_inlandsis_mod, ONLY : surf_inlandsis
[3927]	56	#endif
[3940]	57
[3927]	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
[4523]	75	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs
[3927]	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
[4529]	80	REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS
[3927]	81	REAL, DIMENSION(klon), INTENT(IN) :: ps
[4523]	82	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1
[3927]	83	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
	84	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
	85	#ifdef ISO
[4143]	86	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
	87	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
[3927]	88	#endif
	89
	90	LOGICAL, INTENT(IN) :: debut !true if first step
	91	LOGICAL, INTENT(IN) :: lafin !true if last step
	92	REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
	93	REAL, DIMENSION(klon), INTENT(IN) :: rmu0
	94	REAL, DIMENSION(klon), INTENT(IN) :: lwdownm !ylwdown
	95	REAL, DIMENSION(klon), INTENT(IN) :: albedo !mean albedo
	96	REAL, DIMENSION(klon), INTENT(IN) :: pphi1
[3940]	97	REAL, DIMENSION(klon), INTENT(IN) :: alt !mean altitude of the grid box
[3927]	98	REAL, DIMENSION(klon), INTENT(IN) :: slope !mean slope in grid box
	99	REAL, DIMENSION(klon), INTENT(IN) :: cloudf !total cloud fraction
	100
	101	! In/Output variables
	102	!****************************************************************************************
	103	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
	104	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
	105	REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
	106	#ifdef ISO
	107	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow, xtsol
	108	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: Rland_ice
	109	#endif
	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
[3940]	118	REAL, DIMENSION(6), INTENT(IN) :: SFRWL
	119	REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir,alb_dif
[3927]	120	!albedo SB <<<
	121	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
[4523]	122	REAL, DIMENSION(klon), INTENT(OUT) :: fluxbs
[3927]	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
[4143]	134	REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
[3927]	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
[3940]	146	REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar
	147	INTEGER :: i,j,nt
	148	REAL, DIMENSION(klon) :: fqfonte,ffonte
[4285]	149	REAL, DIMENSION(klon) :: run_off_lic_frac
[3927]	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	REAL, DIMENSION(klon) :: emis_new !Emissivity
	167	REAL, DIMENSION(klon) :: swdown,lwdown
[3940]	168	REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
	169	REAL, DIMENSION(klon) :: erod !erosion of surface snow (flux, kg/m2/s like evap)
	170	REAL, DIMENSION(klon) :: zsl_height, wind_velo !surface layer height, wind spd
[3927]	171	REAL, DIMENSION(klon) :: dens_air, snow_cont_air !air density; snow content air
	172	REAL, DIMENSION(klon) :: alb_soil !albedo of underlying ice
	173	REAL, DIMENSION(klon) :: pexner !Exner potential
	174	REAL :: pref
[3940]	175	REAL, DIMENSION(klon,nsoilmx) :: tsoil0 !modif
	176	REAL :: dtis ! subtimestep
	177	LOGICAL :: debut_is, lafin_is ! debut and lafin for inlandsis
[3927]	178
	179	CHARACTER (len = 20) :: modname = 'surf_landice'
	180	CHARACTER (len = 80) :: abort_message
	181
	182
[3940]	183	REAL,DIMENSION(klon) :: alb1,alb2
[4523]	184	REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
[3940]	185	REAL, DIMENSION (klon,6) :: alb6
[4835]	186	REAL :: esalt
[4529]	187	REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
[4835]	188	REAL :: tau_dens, maxerosion, fluxbs_2
	189	REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
[4672]	190
[3927]	191	! End definition
	192	!****************************************************************************************
	193	!FC
	194	!FC
	195	REAL,SAVE :: alb_vis_sno_lic
	196	!$OMP THREADPRIVATE(alb_vis_sno_lic)
	197	REAL,SAVE :: alb_nir_sno_lic
	198	!$OMP THREADPRIVATE(alb_nir_sno_lic)
	199	LOGICAL, SAVE :: firstcall = .TRUE.
	200	!$OMP THREADPRIVATE(firstcall)
	201
	202
[3940]	203	!FC firtscall initializations
	204	!******************************************************************************************
[3927]	205	#ifdef ISO
	206	#ifdef ISOVERIF
	207	! write(,) 'surf_land_ice 1499'
	208	do i=1,knon
	209	if (iso_eau.gt.0) then
	210	call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
	211	& 'surf_land_ice 126',errmax,errmaxrel)
	212	endif !if (iso_eau.gt.0) then
	213	enddo !do i=1,knon
	214	#endif
	215	#endif
	216
	217	IF (firstcall) THEN
	218	alb_vis_sno_lic=0.77
	219	CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
	220	PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
	221	alb_nir_sno_lic=0.77
	222	CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
	223	PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
[3940]	224
[3927]	225	firstcall=.false.
	226	ENDIF
[3940]	227	!******************************************************************************************
[4523]	228
[3927]	229	! Initialize output variables
	230	alb3(:) = 999999.
	231	alb2(:) = 999999.
	232	alb1(:) = 999999.
[4523]	233	fluxbs(:)=0.
[3927]	234	runoff(:) = 0.
	235	!****************************************************************************************
	236	! Calculate total absorbed radiance at surface
	237	!
	238	!****************************************************************************************
	239	radsol(:) = 0.0
	240	radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
	241
	242	!****************************************************************************************
[4523]	243
	244	!****************************************************************************************
[3940]	245	! landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...
	246	! landice_opt = 1 : prepare and call INterace Lmdz SISvat (INLANDSIS)
[3927]	247	!****************************************************************************************
[3940]	248
	249
	250	IF (landice_opt .EQ. 1) THEN
	251
	252	!****************************************************************************************
	253	! CALL to INLANDSIS interface
	254	!****************************************************************************************
	255	#ifdef CPP_INLANDSIS
	256
	257	#ifdef ISO
	258	CALL abort_gcm('surf_landice 235','isotopes pas dans INLANDSIS',1)
	259	#endif
	260
	261	debut_is=debut
	262	lafin_is=.false.
	263	! Suppose zero surface speed
	264	u0(:) = 0.0
	265	v0(:) = 0.0
	266
	267
	268	CALL calcul_flux_wind(knon, dtime, &
	269	u0, v0, u1, v1, gustiness, cdragm, &
	270	AcoefU, AcoefV, BcoefU, BcoefV, &
	271	p1lay, temp_air, &
	272	flux_u1, flux_v1)
	273
[3927]	274
[3940]	275	! Set constants and compute some input for SISVAT
	276	! = 1000 hPa
	277	! and calculate incoming flux for SW and LW interval: swdown, lwdown
[3927]	278	swdown(:) = 0.0
	279	lwdown(:) = 0.0
[3940]	280	snow_cont_air(:) = 0. ! the snow content in air is not a prognostic variable of the model
	281	alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
	282	ustar(:) = 0.
	283	pref = 100000.
[3927]	284	DO i = 1, knon
	285	swdown(i) = swnet(i)/(1-albedo(i))
	286	lwdown(i) = lwdownm(i)
	287	wind_velo(i) = u1(i)2 + v1(i)2
	288	wind_velo(i) = wind_velo(i)**0.5
	289	pexner(i) = (p1lay(i)/pref)**(RD/RCPD)
	290	dens_air(i) = p1lay(i)/RD/temp_air(i) ! dry air density
[3940]	291	zsl_height(i) = pphi1(i)/RG
	292	tsoil0(i,:) = tsoil(i,:)
	293	ustar(i)= (cdragm(i)(wind_velo(i)2))*0.5
[3927]	294	END DO
[3940]	295
[3927]	296
[3940]	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
[3927]	332	#else
[3940]	333	abort_message='Pb de coherence: landice_opt = 1 mais CPP_INLANDSIS = .false.'
[3927]	334	CALL abort_physic(modname,abort_message,1)
	335	#endif
	336
[3940]	337
	338	ELSE
	339
[3927]	340	!****************************************************************************************
	341	! Soil calculations
	342	!
	343	!****************************************************************************************
[3940]	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
[3927]	350	IF (soil_model) THEN
[3975]	351	CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
[4523]	352	& longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
[3927]	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	!****************************************************************************************
[3940]	365	! beta(:) = 1.0
	366	! dif_grnd(:) = 0.0
[3927]	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
	382	#ifdef ISO
	383	! verif
	384	#ifdef ISOVERIF
[4033]	385	!write(,) 'surf_land_ice 1499'
[3927]	386	do i=1,knon
	387	if (iso_eau.gt.0) then
	388	if (snow(i).gt.ridicule) then
	389	call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
	390	& 'surf_land_ice 1151',errmax,errmaxrel)
	391	endif !if ((snow(i).gt.ridicule)) then
	392	endif !if (iso_eau.gt.0) then
	393	enddo !do i=1,knon
	394	#endif
	395	!#ifdef ISOVERIF
	396
	397	do i=1,knon
	398	snow_prec(i)=snow(i)
	399	do ixt=1,niso
	400	xtsnow_prec(ixt,i)=xtsnow(ixt,i)
	401	enddo !do ixt=1,niso
	402	! initialisation:
	403	fq_fonte_diag(i)=0.0
	404	fqfonte_diag(i)=0.0
	405	snow_evap_diag(i)=0.0
	406	enddo !do i=1,knon
	407	#endif
	408	!#ifdef ISO
	409	CALL calcul_flux_wind(knon, dtime, &
	410	u0, v0, u1, v1, gustiness, cdragm, &
	411	AcoefU, AcoefV, BcoefU, BcoefV, &
	412	p1lay, temp_air, &
	413	flux_u1, flux_v1)
	414
[4523]	415
[3927]	416	!****************************************************************************************
[4523]	417	! Calculate albedo
	418	!
	419	!****************************************************************************************
	420
	421	!
	422	!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
	423	! alb1(1 : knon) = 0.6 !IM cf FH/GK
	424	! alb1(1 : knon) = 0.82
	425	! alb1(1 : knon) = 0.77 !211003 Ksta0.77
	426	! alb1(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5
	427	!IM: KstaTER0.77 & LMD_ARMIP6
	428
	429	! Attantion: alb1 and alb2 are not the same!
	430	alb1(1:knon) = alb_vis_sno_lic
	431	alb2(1:knon) = alb_nir_sno_lic
	432
	433
	434	!****************************************************************************************
	435	! Rugosity
	436	!
	437	!****************************************************************************************
	438	z0m = z0m_landice
	439	z0h = z0h_landice
	440	!z0m = SQRT(z0m2+rugoro2)
	441
	442
[4835]	443	!****************************************************************************************
	444	! Simple blowing snow param
	445	!****************************************************************************************
	446	! we proceed in 2 steps:
	447	! first we erode - if possible -the accumulated snow during the time step
	448	! then we update the density of the underlying layer and see if we can also erode
	449	! this layer
[4523]	450
[4529]	451
[4835]	452	if (ok_bs) then
	453	fluxbs(:)=0.
	454	do j=1,klon
	455	ws1(j)=(u1(j)2+v1(j)2)**0.5
	456	ustar(j)=(cdragm(j)(u1(j)2+v1(j)2))*0.5
	457	rhod(j)=p1lay(j)/RD/temp_air(j)
	458	ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
	459	enddo
	460
	461	! 1st step: erosion of fresh snow accumulated during the time step
	462	do j=1, knon
	463
	464	rhos(j)=rhofresh_bs
	465	! blowing snow flux formula used in MAR
	466	ustart(j)=ustart0(j)exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))exp(max(0.,rhos(j)-rhohard_bs))
	467	! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
	468	! computation of qbs at the top of the saltation layer
	469	! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
	470	esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
	471	hsalt(j)=0.08436(max(1.e-6,ustar(j))*1.27)
	472	qsalt(j)=(max(ustar(j)2-ustart(j)2,0.))/(RGhsalt(j))esalt
	473	! calculation of erosion (flux positive towards the surface here)
	474	! consistent with implicit resolution of turbulent mixing equation
	475	! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
	476	! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
	477	! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
	478	! (rhoqsalthsalt)
	479	! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
	480	maxerosion=min(precip_snow(j),hsalt(j)qsalt(j)rhod(j)/tau_eqsalt_bs)
	481
	482	fluxbs(j)=rhod(j)ws1(j)cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
	483	/ (1.-rhod(j)ws1(j)cdragm(j)BcoefQBS(j)dtime)
	484	!fluxbs(j)= rhodws1(j)cdragm(j)*(qbs1(j)-qsalt(j))
	485	fluxbs(j)=max(-maxerosion,fluxbs(j))
	486	enddo
	487
	488
	489	! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
	490	! this is done through the routine albsno
	491	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:))
	492
	493	! 2nd step:
[4529]	494	! computation of threshold friction velocity
	495	! which depends on surface snow density
[4672]	496	do j = 1, knon
[4523]	497	! estimation of snow density
	498	! snow density increases with snow age and
[4672]	499	! increases even faster in case of sedimentation of blowing snow or rain
[4835]	500	tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - &
	501	abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
	502	rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)(1.-exp(-agesno(j)86400.0/tau_dens))
[4523]	503	! blowing snow flux formula used in MAR
[4835]	504	ustart(j)=ustart0(j)exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))exp(max(0.,rhos(j)-rhohard_bs))
	505	! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
	506	! computation of qbs at the top of the saltation layer
	507	! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
	508	esalt=c_esalt_bs*ustar(j)
	509	hsalt(j)=0.08436(max(1.e-6,ustar(j))*1.27)
	510	qsalt(j)=(max(ustar(j)2-ustart(j)2,0.))/(RGhsalt(j))esalt
	511	! calculation of erosion (flux positive towards the surface here)
	512	! consistent with implicit resolution of turbulent mixing equation
	513	! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
	514	! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
	515	! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
	516	! (rhoqsalthsalt)
	517	maxerosion=hsalt(j)qsalt(j)rhod(j)/tau_eqsalt_bs
	518	fluxbs_2=rhod(j)ws1(j)cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
	519	/ (1.-rhod(j)ws1(j)cdragm(j)BcoefQBS(j)dtime)
	520	!fluxbs_2= rhodws1(j)cdragm(j)*(qbs1(j)-qsalt(j))
	521	fluxbs_2=max(-maxerosion,fluxbs_2)
	522	fluxbs(j)=fluxbs(j)+fluxbs_2
[4523]	523	enddo
	524
[4835]	525	! for outputs
	526	do j=1, knon
	527	i = knindex(j)
	528	zxustartlic(i) = ustart(j)
	529	zxrhoslic(i) = rhos(j)
	530	zxqsaltlic(i)=qsalt(j)
	531	enddo
[4529]	532
	533
[4835]	534	else
	535	! those lines are useful to calculate the snow age
	536	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
[4529]	537
[4835]	538	endif ! if ok_bs
[4672]	539
[4523]	540
	541
	542
	543
[4835]	544
[4523]	545	!****************************************************************************************
[4835]	546	! Calculate snow amount
[3927]	547	!
	548	!****************************************************************************************
[4523]	549	IF (ok_bs) THEN
[4528]	550	precip_totsnow(:)=precip_snow(:)+precip_bs(:)
	551	evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
[4523]	552	ELSE
[4528]	553	precip_totsnow(:)=precip_snow(:)
	554	evap_totsnow(:)=evap(:)
[4523]	555	ENDIF
	556
[3940]	557	CALL fonte_neige(knon, is_lic, knindex, dtime, &
[4523]	558	tsurf, precip_rain, precip_totsnow, &
	559	snow, qsol, tsurf_new, evap_totsnow &
[3927]	560	#ifdef ISO
	561	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
	562	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
	563	#endif
	564	& )
	565
	566
	567	#ifdef ISO
	568	#ifdef ISOVERIF
	569	do i=1,knon
	570	if (iso_eau.gt.0) then
	571	call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
	572	& 'surf_landice_mod 217',errmax,errmaxrel)
	573	endif !if (iso_eau.gt.0) then
	574	enddo !do i=1,knon
	575	#endif
	576	!#ifdef ISOVERIF
	577
	578	call calcul_iso_surf_lic_vectall(klon,knon, &
	579	& evap,snow_evap_diag,Tsurf_new,snow, &
	580	& fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
	581	& precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, &
	582	& xtspechum,spechum,ps,Rland_ice, &
	583	& xtevap,xtsnow,fqcalving_diag, &
	584	& knindex,is_lic,run_off_lic_diag,coeff_rel_diag &
	585	& )
	586
	587	! call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
	588
	589	#endif
	590	!#ifdef ISO
	591
	592
[4523]	593	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
	594	zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[3940]	595
[3927]	596
[3940]	597	END IF ! landice_opt
	598
	599
[3927]	600	!****************************************************************************************
	601	! Send run-off on land-ice to coupler if coupled ocean.
[3940]	602	! run_off_lic has been calculated in fonte_neige or surf_inlandsis
[4285]	603	! If landice_opt>=2, corresponding call is done from surf_land_orchidee
[3927]	604	!****************************************************************************************
[4285]	605	IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
	606	! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
	607	run_off_lic_frac(:)=0.0
	608	DO j = 1, knon
	609	i = knindex(j)
	610	run_off_lic_frac(j) = pctsrf(i,is_lic)
	611	ENDDO
[4523]	612
[4285]	613	CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
[3927]	614	ENDIF
	615
	616	! transfer runoff rate [kg/m2/s](!) to physiq for output
	617	runoff(1:knon)=run_off_lic(1:knon)/dtime
	618
	619	snow_o=0.
	620	zfra_o = 0.
	621	DO j = 1, knon
	622	i = knindex(j)
	623	snow_o(i) = snow(j)
	624	zfra_o(i) = zfra(j)
	625	ENDDO
	626
	627
	628	!albedo SB >>>
	629	select case(NSW)
	630	case(2)
	631	alb_dir(1:knon,1)=alb1(1:knon)
	632	alb_dir(1:knon,2)=alb2(1:knon)
	633	case(4)
	634	alb_dir(1:knon,1)=alb1(1:knon)
	635	alb_dir(1:knon,2)=alb2(1:knon)
	636	alb_dir(1:knon,3)=alb2(1:knon)
	637	alb_dir(1:knon,4)=alb2(1:knon)
	638	case(6)
	639	alb_dir(1:knon,1)=alb1(1:knon)
	640	alb_dir(1:knon,2)=alb1(1:knon)
	641	alb_dir(1:knon,3)=alb1(1:knon)
	642	alb_dir(1:knon,4)=alb2(1:knon)
	643	alb_dir(1:knon,5)=alb2(1:knon)
	644	alb_dir(1:knon,6)=alb2(1:knon)
[3940]	645
	646	IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
	647	alb_dir(1:knon,1)=alb6(1:knon,1)
	648	alb_dir(1:knon,2)=alb6(1:knon,2)
	649	alb_dir(1:knon,3)=alb6(1:knon,3)
	650	alb_dir(1:knon,4)=alb6(1:knon,4)
	651	alb_dir(1:knon,5)=alb6(1:knon,5)
	652	alb_dir(1:knon,6)=alb6(1:knon,6)
	653	ENDIF
	654
[3927]	655	end select
	656	alb_dif=alb_dir
	657	!albedo SB <<<
	658
	659
	660	END SUBROUTINE surf_landice
	661	!
	662	!****************************************************************************************
	663	!
	664	END MODULE surf_landice_mod
	665
	666
	667

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

Download in other formats:

Original Format