Context Navigation

surf_landice_mod.F90 @ 4582

Last change on this file since 4582 was 4538, checked in by evignon, 18 months ago
fix bug in surf_landice_mod
File size: 24.1 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
[4523]	37	USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
	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
[4529]	52	USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
[3927]	53
[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
[3975]	60
[3927]	61	! INCLUDE "indicesol.h"
	62	INCLUDE "dimsoil.h"
	63	INCLUDE "YOMCST.h"
	64	INCLUDE "clesphys.h"
	65
	66	! Input variables
	67	!****************************************************************************************
	68	INTEGER, INTENT(IN) :: itime, knon
	69	INTEGER, DIMENSION(klon), INTENT(in) :: knindex
	70	REAL, INTENT(in) :: dtime
	71	REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiance
	72	REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiance
	73	REAL, DIMENSION(klon), INTENT(IN) :: tsurf
	74	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
	75	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
[4523]	76	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs
[3927]	77	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
	78	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
	79	REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
	80	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
[4529]	81	REAL, DIMENSION(klon), INTENT(IN) :: AcoefQBS, BcoefQBS
[3927]	82	REAL, DIMENSION(klon), INTENT(IN) :: ps
[4523]	83	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness, qbs1
[3927]	84	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
	85	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
	86	#ifdef ISO
[4143]	87	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
	88	REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
[3927]	89	#endif
	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
[3940]	98	REAL, DIMENSION(klon), INTENT(IN) :: alt !mean altitude of the grid box
[3927]	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	! Output variables
	113	!****************************************************************************************
	114	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
	115	REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h
	116	!albedo SB >>>
	117	! REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! new albedo in visible SW interval
	118	! REAL, DIMENSION(klon), INTENT(OUT) :: alb2 ! new albedo in near IR interval
[3940]	119	REAL, DIMENSION(6), INTENT(IN) :: SFRWL
	120	REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir,alb_dif
[3927]	121	!albedo SB <<<
	122	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
[4523]	123	REAL, DIMENSION(klon), INTENT(OUT) :: fluxbs
[3927]	124	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
	125	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
	126	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
	127
	128	REAL, DIMENSION(klon), INTENT(OUT) :: alb3
	129	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow !column water in snow [kg/m2]
	130	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt !Snow height (m)
	131	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice
	132	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow
	133	REAL, DIMENSION(klon), INTENT(OUT) :: runoff !Land ice runoff
	134	#ifdef ISO
[4143]	135	REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
[3927]	136	! real, DIMENSION(niso,klon) :: xtrun_off_lic_0_diag ! est une variable globale de
	137	! fonte_neige
	138	#endif
	139
	140
	141	! Local variables
	142	!****************************************************************************************
	143	REAL, DIMENSION(klon) :: soilcap, soilflux
	144	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
	145	REAL, DIMENSION(klon) :: zfra, alb_neig
	146	REAL, DIMENSION(klon) :: radsol
[3940]	147	REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar
	148	INTEGER :: i,j,nt
	149	REAL, DIMENSION(klon) :: fqfonte,ffonte
[4285]	150	REAL, DIMENSION(klon) :: run_off_lic_frac
[3927]	151	#ifdef ISO
	152	real, parameter :: t_coup = 273.15
	153	real, dimension(klon) :: fqfonte_diag
	154	real, dimension(klon) :: fq_fonte_diag
	155	real, dimension(klon) :: snow_evap_diag
	156	real, dimension(klon) :: fqcalving_diag
	157	real max_eau_sol_diag
	158	real, dimension(klon) :: runoff_diag
	159	real, dimension(klon) :: run_off_lic_diag
	160	real :: coeff_rel_diag
	161	integer ixt
	162	REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
	163	REAL, DIMENSION(klon) :: snow_prec,qsol_prec
	164	! real, DIMENSION(klon) :: run_off_lic_0_diag
	165	#endif
	166
	167	REAL, DIMENSION(klon) :: emis_new !Emissivity
	168	REAL, DIMENSION(klon) :: swdown,lwdown
[3940]	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
[3927]	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
[3940]	176	REAL, DIMENSION(klon,nsoilmx) :: tsoil0 !modif
	177	REAL :: dtis ! subtimestep
	178	LOGICAL :: debut_is, lafin_is ! debut and lafin for inlandsis
[3927]	179
	180	CHARACTER (len = 20) :: modname = 'surf_landice'
	181	CHARACTER (len = 80) :: abort_message
	182
	183
[3940]	184	REAL,DIMENSION(klon) :: alb1,alb2
[4523]	185	REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
[3940]	186	REAL, DIMENSION (klon,6) :: alb6
[4529]	187	REAL :: rho0, rhoice, ustart0, hsalt, esalt, rhod
	188	REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
[4523]	189	REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
[4529]	190	REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt
[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	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
	422
	423
	424	! EV: following lines are obsolete since we set alb1 and alb2 to constant values
	425	! I therefore comment them
	426	! alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
	427	! 0.6 * (1.0-zfra(1:knon))
	428	!
	429	!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
	430	! alb1(1 : knon) = 0.6 !IM cf FH/GK
	431	! alb1(1 : knon) = 0.82
	432	! alb1(1 : knon) = 0.77 !211003 Ksta0.77
	433	! alb1(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5
	434	!IM: KstaTER0.77 & LMD_ARMIP6
	435
	436	! Attantion: alb1 and alb2 are not the same!
	437	alb1(1:knon) = alb_vis_sno_lic
	438	alb2(1:knon) = alb_nir_sno_lic
	439
	440
	441	!****************************************************************************************
	442	! Rugosity
	443	!
	444	!****************************************************************************************
	445	z0m = z0m_landice
	446	z0h = z0h_landice
	447	!z0m = SQRT(z0m2+rugoro2)
	448
	449
	450
	451	! Simple blowing snow param
	452	if (ok_bs) then
	453	ustart0 = 0.211
	454	rhoice = 920.0
	455	rho0 = 200.0
	456	rhomax=450.0
	457	rhohard=400.0
	458	tau_dens0=86400.0*10. ! 10 days by default, in s
	459	tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
[4529]	460
	461	! computation of threshold friction velocity
	462	! which depends on surface snow density
[4523]	463	do i = 1, knon
	464	! estimation of snow density
	465	! snow density increases with snow age and
	466	! increases even faster in case of sedimentation of blowing snow
	467	tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
	468	rhos(i)=rho0+(rhohard-rho0)(1.-exp(-agesno(i)86400.0/tau_dens))
	469	! blowing snow flux formula used in MAR
	470	ws1(i)=(u1(i)2+v1(i)2)**0.5
	471	ustar(i)=(cdragm(i)(u1(i)2+v1(i)2))*0.5
	472	ustart(i)=ustart0exp(max(rhoice/rho0-rhoice/rhos(i),0.))exp(max(0.,rhos(i)-rhomax))
	473	! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
	474	! rhohard<450)
	475	enddo
[4529]	476
	477	! computation of qbs at the top of the saltation layer
	478	! two formulations possible
	479	! pay attention that qbs is a mixing ratio and has to be converted
	480	! to specific content
	481
	482	if (iflag_saltation_bs .eq. 1) then
	483	! expression from CRYOWRF (Sharma et al. 2022)
	484	aa=2.6
	485	bb=2.5
	486	cc=2.0
	487	lambdasalt=0.45
	488	do i =1, knon
	489	rhod=p1lay(i)/RD/temp_air(i)
	490	nunu=max(ustar(i)/ustart(i),1.e-3)
	491	fluxsalt=rhod/RG(ustar(i)3)(1.-nunu*(-2)) &
	492	(aa+bbnunu(-2)+ccnunu**(-1))
	493	csalt=fluxsalt/(2.8*ustart(i))
	494	hsalt=0.08436ustar(i)*1.27
	495	qsalt(i)=1./rhodcsaltlambdasaltRG/(max(ustar(i)*2,1E-6)) &
	496	* exp(-lambdasaltRGhsalt/max(ustar(i)**2,1E-6))
	497	qsalt(i)=max(qsalt(i),0.)
	498	enddo
[4523]	499
[4529]	500
	501	else
	502	! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
	503	do i=1, knon
	504	esalt=1./(3.25*max(ustar(i),0.001))
	505	hsalt=0.08436ustar(i)*1.27
	506	qsalt(i)=(max(ustar(i)2-ustart(i)2,0.))/(RGhsalt)esalt
	507	!ep=qsaltcdragm(i)sqrt(u1(i)2+v1(i)2)
	508	enddo
	509	endif
	510
	511	! calculation of erosion (emission flux towards the first atmospheric level)
	512	! consistent with implicit resolution of turbulent mixing equation
	513	do i=1, knon
	514	rhod=p1lay(i)/RD/temp_air(i)
	515	fluxbs(i)=rhodws1(i)cdragm(i)zeta_bs(AcoefQBS(i)-qsalt(i)) &
[4538]	516	/ (1.-rhodws1(i)zeta_bscdragm(i)BcoefQBS(i)*dtime)
[4529]	517	!fluxbs(i)= zeta_bsrhodws1(i)cdragm(i)(qbs1(i)-qsalt(i))
	518	enddo
	519
[4523]	520	! for outputs
	521	do j = 1, knon
	522	i = knindex(j)
	523	zxustartlic(i) = ustart(j)
	524	zxrhoslic(i) = rhos(j)
	525	enddo
	526
	527	endif
	528
	529
	530
	531	!****************************************************************************************
	532	! Calculate surface snow amount
[3927]	533	!
	534	!****************************************************************************************
[4523]	535	IF (ok_bs) THEN
[4528]	536	precip_totsnow(:)=precip_snow(:)+precip_bs(:)
	537	evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
[4523]	538	ELSE
[4528]	539	precip_totsnow(:)=precip_snow(:)
	540	evap_totsnow(:)=evap(:)
[4523]	541	ENDIF
	542
[3940]	543	CALL fonte_neige(knon, is_lic, knindex, dtime, &
[4523]	544	tsurf, precip_rain, precip_totsnow, &
	545	snow, qsol, tsurf_new, evap_totsnow &
[3927]	546	#ifdef ISO
	547	& ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag &
	548	& ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
	549	#endif
	550	& )
	551
	552
	553	#ifdef ISO
	554	#ifdef ISOVERIF
	555	do i=1,knon
	556	if (iso_eau.gt.0) then
	557	call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
	558	& 'surf_landice_mod 217',errmax,errmaxrel)
	559	endif !if (iso_eau.gt.0) then
	560	enddo !do i=1,knon
	561	#endif
	562	!#ifdef ISOVERIF
	563
	564	call calcul_iso_surf_lic_vectall(klon,knon, &
	565	& evap,snow_evap_diag,Tsurf_new,snow, &
	566	& fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
	567	& precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, &
	568	& xtspechum,spechum,ps,Rland_ice, &
	569	& xtevap,xtsnow,fqcalving_diag, &
	570	& knindex,is_lic,run_off_lic_diag,coeff_rel_diag &
	571	& )
	572
	573	! call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
	574
	575	#endif
	576	!#ifdef ISO
	577
	578
[4523]	579	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
	580	zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[3940]	581
[3927]	582
[3940]	583	END IF ! landice_opt
	584
	585
[3927]	586	!****************************************************************************************
	587	! Send run-off on land-ice to coupler if coupled ocean.
[3940]	588	! run_off_lic has been calculated in fonte_neige or surf_inlandsis
[4285]	589	! If landice_opt>=2, corresponding call is done from surf_land_orchidee
[3927]	590	!****************************************************************************************
[4285]	591	IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
	592	! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
	593	run_off_lic_frac(:)=0.0
	594	DO j = 1, knon
	595	i = knindex(j)
	596	run_off_lic_frac(j) = pctsrf(i,is_lic)
	597	ENDDO
[4523]	598
[4285]	599	CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
[3927]	600	ENDIF
	601
	602	! transfer runoff rate [kg/m2/s](!) to physiq for output
	603	runoff(1:knon)=run_off_lic(1:knon)/dtime
	604
	605	snow_o=0.
	606	zfra_o = 0.
	607	DO j = 1, knon
	608	i = knindex(j)
	609	snow_o(i) = snow(j)
	610	zfra_o(i) = zfra(j)
	611	ENDDO
	612
	613
	614	!albedo SB >>>
	615	select case(NSW)
	616	case(2)
	617	alb_dir(1:knon,1)=alb1(1:knon)
	618	alb_dir(1:knon,2)=alb2(1:knon)
	619	case(4)
	620	alb_dir(1:knon,1)=alb1(1:knon)
	621	alb_dir(1:knon,2)=alb2(1:knon)
	622	alb_dir(1:knon,3)=alb2(1:knon)
	623	alb_dir(1:knon,4)=alb2(1:knon)
	624	case(6)
	625	alb_dir(1:knon,1)=alb1(1:knon)
	626	alb_dir(1:knon,2)=alb1(1:knon)
	627	alb_dir(1:knon,3)=alb1(1:knon)
	628	alb_dir(1:knon,4)=alb2(1:knon)
	629	alb_dir(1:knon,5)=alb2(1:knon)
	630	alb_dir(1:knon,6)=alb2(1:knon)
[3940]	631
	632	IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
	633	alb_dir(1:knon,1)=alb6(1:knon,1)
	634	alb_dir(1:knon,2)=alb6(1:knon,2)
	635	alb_dir(1:knon,3)=alb6(1:knon,3)
	636	alb_dir(1:knon,4)=alb6(1:knon,4)
	637	alb_dir(1:knon,5)=alb6(1:knon,5)
	638	alb_dir(1:knon,6)=alb6(1:knon,6)
	639	ENDIF
	640
[3927]	641	end select
	642	alb_dif=alb_dir
	643	!albedo SB <<<
	644
	645
	646	END SUBROUTINE surf_landice
	647	!
	648	!****************************************************************************************
	649	!
	650	END MODULE surf_landice_mod
	651
	652
	653

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

Download in other formats:

Original Format