Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/explicitsnow.f90 @ 7485

Last change on this file since 7485 was 7266, checked in by agnes.ducharne, 3 years ago
Integrated snow related patches of the trunk (r6402 and r6523). Weak changes in sowy areas with a 5d test.
Property svn:keywords set to `Date Revision HeadURL`
File size: 73.0 KB

Rev	Line
[4282]	1	! ================================================================================================================================
	2	! MODULE : explicitsnow
	3	!
[4470]	4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
[4282]	5	!
	6	! LICENCE : IPSL (2006)
	7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
	8	!
	9	!>\BRIEF Computes hydrologic snow processes on continental points.
[2223]	10	!!
[4282]	11	!!\n DESCRIPTION: This module computes hydrologic snow processes on continental points.
[2223]	12	!!
[4282]	13	!! RECENT CHANGES:
	14	!!
	15	!! REFERENCES : None
	16	!!
	17	!! SVN :
	18	!! $HeadURL$
	19	!! $Date$
	20	!! $Revision$
	21	!! \n
	22	!_ ================================================================================================================================
	23
[2223]	24	MODULE explicitsnow
[2650]	25	USE ioipsl_para
[2223]	26	USE constantes_soil
	27	USE constantes
[4646]	28	USE time, ONLY : one_day, dt_sechiba
[2223]	29	USE pft_parameters
	30	USE qsat_moisture
[4281]	31	USE sechiba_io_p
[3410]	32	USE xios_orchidee
[2223]	33
	34	IMPLICIT NONE
	35
	36	! Public routines :
	37	PRIVATE
[2650]	38	PUBLIC :: explicitsnow_main, explicitsnow_initialize, explicitsnow_finalize
[2223]	39
	40	CONTAINS
	41
[2650]	42	!================================================================================================================================
	43	!! SUBROUTINE : explicitsnow_initialize
[2223]	44	!!
[2650]	45	!>\BRIEF Read variables for explictsnow module from restart file
	46	!!
[4282]	47	!! DESCRIPTION : Read variables for explictsnow module from restart file
[2223]	48	!!
[2650]	49	!! \n
	50	!_
	51	!================================================================================================================================
	52	SUBROUTINE explicitsnow_initialize( kjit, kjpindex, rest_id, snowrho, &
[3059]	53	snowtemp, snowdz, snowheat, snowgrain)
[2650]	54
	55	!! 0.1 Input variables
	56	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
	57	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
	58	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
	59
	60	!! 0.2 Output variables
[3410]	61	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho !! Snow density (Kg/m^3)
[2650]	62	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp !! Snow temperature
	63	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz !! Snow layer thickness
[3410]	64	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat !! Snow heat content (J/m^2)
	65	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain !! Snow grainsize (m)
[2650]	66
	67
	68	!! 1. Read from restart file
	69	CALL restget_p (rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit,.TRUE.,snowrho, "gather", nbp_glo, index_g)
	70	CALL setvar_p (snowrho, val_exp, 'Snow Density profile', xrhosmin)
	71
	72	CALL restget_p (rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit,.TRUE.,snowtemp, "gather", nbp_glo, index_g)
	73	CALL setvar_p (snowtemp, val_exp, 'Snow Temperature profile', tp_00)
	74
	75	CALL restget_p (rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit,.TRUE.,snowdz, "gather", nbp_glo, index_g)
	76	CALL setvar_p (snowdz, val_exp, 'Snow depth profile', 0.0)
	77
	78	CALL restget_p (rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit,.TRUE.,snowheat, "gather", nbp_glo, index_g)
	79	CALL setvar_p (snowheat, val_exp, 'Snow Heat profile', 0.0)
	80
	81	CALL restget_p (rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit,.TRUE.,snowgrain, "gather", nbp_glo, index_g)
	82	CALL setvar_p (snowgrain, val_exp, 'Snow grain profile', 0.0)
	83
	84	END SUBROUTINE explicitsnow_initialize
	85
	86
	87	!================================================================================================================================
	88	!! SUBROUTINE : explicitsnow_main
[2223]	89	!!
[4282]	90	!>\BRIEF Main call for snow calculations
[2650]	91	!!
[4282]	92	!! DESCRIPTION : Main routine for calculation of the snow processes with explictsnow module.
[2223]	93	!!
[2650]	94	!! RECENT CHANGE(S) : None
	95	!!
	96	!! MAIN OUTPUT VARIABLE(S): None
	97	!!
	98	!! REFERENCE(S) :
	99	!!
	100	!! FLOWCHART : None
[2223]	101	!! \n
[2650]	102	!_
	103	!================================================================================================================================
	104	SUBROUTINE explicitsnow_main(kjpindex, precip_rain, precip_snow, temp_air, pb, & ! in
	105	u, v, temp_sol_new, soilcap, pgflux, & ! in
[3269]	106	frac_nobio, totfrac_nobio,gtemp, & ! in
[5091]	107	lambda_snow, cgrnd_snow, dgrnd_snow, contfrac, & ! in
[2650]	108	vevapsno, snow_age, snow_nobio_age,snow_nobio, snowrho, & ! inout
	109	snowgrain, snowdz, snowtemp, snowheat, snow, & ! inout
[3269]	110	temp_sol_add, & ! inout
[2650]	111	snowliq, subsnownobio, grndflux, snowmelt, tot_melt, & ! output
[3269]	112	subsinksoil ) ! output
[2650]	113
[2223]	114
	115	!! 0. Variable and parameter declaration
	116
	117	!! 0.1 Input variables
	118	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
	119	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
	120	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall
	121	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
[2650]	122	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
[2223]	123	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: u,v !! Horizontal wind speed
[2650]	124	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature
[2223]	125	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
	126	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pgflux !! Net energy into snowpack
	127	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
	128	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
	129	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: gtemp !! First soil layer temperature
[3269]	130	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	131	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
	132	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
[5091]	133	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
[2223]	134
	135	!! 0.2 Output fields
[2650]	136	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowliq !! Snow liquid content (m)
	137	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
	138	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: grndflux !! Net flux into soil [W/m2]
[3687]	139	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowmelt !! Snow melt [mm/dt_sechiba]
[2223]	140	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt from ice and snow
[2650]	141	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil
[2223]	142
	143	!! 0.3 Modified fields
	144	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
	145	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
	146	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
	147	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
[3410]	148	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
	149	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain !! Snow grainsize (m)
	150	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow layer thickness [m]
[2650]	151	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
[3410]	152	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat !! Snow heat content (J/m^2)
[2223]	153	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
[3269]	154	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
[2223]	155
[2650]	156
[2223]	157	!! 0.4 Local declaration
	158	INTEGER(i_std) :: ji, iv, jj,m,jv
	159	REAL(r_std),DIMENSION (kjpindex) :: snow_depth_tmp
	160	REAL(r_std),DIMENSION (kjpindex) :: snowmelt_from_maxmass
	161	REAL(r_std) :: snowdzm1
	162	REAL(r_std), DIMENSION (kjpindex) :: thrufal !! Water leaving snow pack [kg/m2/s]
	163	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
	164	REAL(r_std), DIMENSION (kjpindex) :: xx !! Temporary
	165	REAL(r_std), DIMENSION (kjpindex) :: snowmelt_tmp,snowmelt_ice,icemelt,temp_sol_new_old
	166	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowdz_old
	167	REAL(r_std), DIMENSION (kjpindex) :: ZLIQHEATXS
[3521]	168	REAL(r_std), DIMENSION (kjpindex) :: ZSNOWEVAPS, ZSNOWDZ,subsnowveg,ZSNOWMELT_XS
	169	REAL(r_std) :: maxmass_snowdepth, snow_d1k
[2650]	170	REAL(r_std), DIMENSION (kjpindex,nsnow) :: WSNOWDZ,SMASS
[3521]	171	REAL(r_std), DIMENSION (kjpindex) :: SMASSC,snowacc,snow_remove
[2650]	172	INTEGER(i_std) :: locjj
	173	REAL(r_std) :: grndflux_tmp
	174	REAL(r_std), DIMENSION (nsnow) :: snowtemp_tmp
	175	REAL(r_std) :: s2flux_tmp,fromsoilflux
	176	REAL(r_std), DIMENSION (kjpindex,nsnow) :: pcapa_snow
	177	REAL(r_std), DIMENSION (kjpindex) :: psnowhmass
	178	REAL(r_std), PARAMETER :: XP00 = 1.E5
[4544]	179	REAL(r_std), DIMENSION (kjpindex) :: snowheattot_begin,snowheattot_end,del_snowheattot !![J/m2]
[5091]	180	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowliq_diag ! Snow liquid content per snow layer [kg/m2]
	181	REAL(r_std), DIMENSION (kjpindex) :: snowliqtot_diag ! Snow liquid content integrated over snow depth [kg/m2]
[4881]	182	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowrho_diag
	183	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowheat_diag
	184	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowgrain_diag
[2223]	185
	186	!! 1. Initialization
[2650]	187
[2223]	188	temp_sol_new_old = temp_sol_new
[3521]	189	snowmelt_ice(:) = zero
	190	icemelt(:) = zero
	191	tot_melt(:) = zero
	192	snowmelt(:) = zero
	193	snowmelt_from_maxmass(:) = zero
	194
[2650]	195
	196	!! 2. on Vegetation
	197	! 2.1 Snow fall
[4544]	198	snowheattot_begin(:)=SUM(snowheat(:,:),2)
[2650]	199	CALL explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
	200	& snowheat,snowgrain,snowtemp,psnowhmass)
	201
	202	! 2.2 calculate the new snow discretization
	203	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
	204
	205	snowdz_old = snowdz
	206
	207	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
	208
	209	! 2.3 Snow heat redistribution
	210	CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
	211
	212	! 2.4 Diagonize water portion of the snow from snow heat content:
	213	DO ji=1, kjpindex
	214	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
	215	snowtemp(ji,:) = snow3ltemp_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:))
	216	snowliq(ji,:) = snow3lliq_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
	217	ELSE
	218	snowliq(ji,:) = zero
	219	snowtemp(ji,:) = tp_00
	220	ENDIF
[2223]	221	END DO
[2650]	222
	223	! 2.5 snow compaction
	224	CALL explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
	225	! Update snow heat
	226	DO ji = 1, kjpindex
	227	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
	228	ENDDO
[3269]	229
[4575]	230	! 2.6 Calculate the snow temperature profile based on heat diffusion
[3269]	231	CALL explicitsnow_profile (kjpindex,cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
[2650]	232
[4575]	233	! 2.7 Test whether snow is existed on the ground or not
[2223]	234	grndflux(:)=0.0
[2650]	235	CALL explicitsnow_gone(kjpindex,pgflux,&
	236	& snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
	237
[4575]	238	! 2.8 Calculate snow melt/refreezing processes
[2650]	239	CALL explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
	240	& snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
	241
	242
[3269]	243	! 2.9 Snow sublimation changing snow thickness
[2223]	244	snow(:) = 0.0
[2650]	245	DO ji=1,kjpindex !domain size
[2223]	246	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
[2650]	247	ENDDO
	248
	249	subsnownobio(:,:) = zero
[2223]	250	subsinksoil(:) = zero
[2650]	251
	252	DO ji=1, kjpindex ! domain size
[2223]	253	IF ( snow(ji) > snowcri ) THEN
[4575]	254	! There is enough snow to split the sublimation demand between the nobio and vegetation parts.
	255	IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
	256	! There is enough snow on the ice fraction to sustain the sublimation demand.
	257	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
	258	ELSE
	259	! There is not enough snow on the ice fraction to sustain the full sublimation demand.
	260	! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
	261	! We do this because the nobio snow cannot deal with a too high sublimation demand, whereas the
	262	! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into
	263	! account in the water budget (see TWBR in hydrol.f90).
	264	subsnownobio(ji,iice) = snow_nobio(ji,iice)
	265	ENDIF
[2223]	266	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
	267	ELSE
[4575]	268	! There is a small amount of snow.
[2223]	269	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
[4575]	270	! The ice fraction is not too small. In this case, we make the hypothesis that the snow is mainly on the ice fraction.
	271	IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
	272	! There is enough snow on the ice fraction to sustain the sublimation demand.
	273	subsnownobio(ji,iice) = vevapsno(ji)
	274	subsnowveg(ji) = zero
	275	ELSE
	276	! There is not enough snow on the ice fraction to sustain the sublimation demand.
	277	! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
	278	! We do this beacuse the nobio snow cannot deal with a too high sublimation demand, whereas the
	279	! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into
	280	! account in the water budget (see TWBR in hydrol.f90).
	281	subsnownobio(ji,iice) = snow_nobio(ji,iice)
	282	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
	283	ENDIF
[2223]	284	ELSE
[4575]	285	! The nobio snow fraction is too small, the vegetation fraction deals with the whole sublimation demand.
[2223]	286	subsnownobio(ji,iice) = zero
	287	subsnowveg(ji) = vevapsno(ji)
[2650]	288	ENDIF
	289	ENDIF
	290
[3059]	291	!! 2.8.1 Check that sublimation on the vegetated fraction is possible.
[2223]	292	IF (subsnowveg(ji) .GT. snow(ji)) THEN
	293	IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN
	294	subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji))
	295	END IF
	296	! Sublimation is thus limited to what is available
	297	subsnowveg(ji) = snow(ji)
	298	snow(ji) = zero
	299	snowdz(ji,:) = 0
	300	snowliq(ji,:) = 0
	301	snowtemp(ji,:) = tp_00
	302	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
	303	ELSE
[2650]	304	! Calculating the snow accumulation
	305	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:)
	306	SMASSC (ji)= 0.0
	307	DO jj=1,nsnow
	308	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
	309	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
	310	ENDDO
	311	! Finding the layer
[3521]	312	locjj=1
[2650]	313	DO jj=1,nsnow-1
	314	IF ((SMASS(ji,jj) .LE. subsnowveg(ji)) .AND. (SMASS(ji,jj+1) .GE. subsnowveg(ji)) ) THEN
	315	locjj=jj+1
	316	ENDIF
	317	ENDDO
	318
	319	! Calculating the removal of snow depth
	320	IF (locjj .EQ. 1) THEN
	321	ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1)
	322	ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
	323	snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
	324	ELSE IF (locjj .GT. 1) THEN
	325	snowacc(ji)=0
	326	DO jj=1,locjj-1
	327	snowacc(ji)=snowacc(ji)+snowdz(ji,jj)*snowrho(ji,jj)
	328	snowdz(ji,jj)=0
	329	ENDDO
	330	ZSNOWEVAPS(ji) = (subsnowveg(ji)-snowacc(ji))/snowrho(ji,locjj)
	331	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWEVAPS(ji)
	332	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
[3521]	333	! ELSE
	334	! ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1)
	335	! ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
	336	! snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
[2650]	337	ENDIF
	338
	339	ENDIF
	340	ENDDO
	341
[3521]	342	! 2.9.1 Snowmelt_from_maxmass
	343
	344	DO ji=1,kjpindex !domain size
	345	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
	346	IF(snow(ji).GT.maxmass_snow)THEN
[3574]	347	IF (snow(ji).GT.(1.2*maxmass_snow)) THEN
	348	! melt faster for points where accumulation is too high
	349	snow_d1k = trois * soilcap(ji) / chalfu0
	350	ELSE
	351	! standard melting
	352	snow_d1k = un * soilcap(ji) / chalfu0
	353	ENDIF
[3521]	354	snowmelt_from_maxmass(ji) = MIN((snow(ji) - maxmass_snow),snow_d1k)
	355	! Calculating the snow accumulation
	356	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) ! WSNOWDZ en kg/m2
	357	SMASSC (ji)= 0.0
	358	DO jj=nsnow,1,-1
	359	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
	360	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
	361	ENDDO
	362
	363	! Finding the layer
	364	locjj = nsnow
	365	DO jj=nsnow,2,-1
	366	IF ((SMASS(ji,jj) .LE. snowmelt_from_maxmass(ji)) .AND. (SMASS(ji,jj-1) .GE. snowmelt_from_maxmass(ji)) ) THEN
	367	locjj=jj-1
	368	ENDIF
	369	ENDDO
	370
	371	! Calculating the removal of snow depth
	372	IF (locjj .EQ. 3) THEN
	373	! ZSNOWMELT_XS : Epaisseur de la couche qui a disparu partiellement
	374	ZSNOWMELT_XS(ji) = snowmelt_from_maxmass(ji)/snowrho(ji,3)
	375	ZSNOWDZ(ji) = snowdz(ji,3) - ZSNOWMELT_XS(ji)
	376	snowdz(ji,3) = MAX(0.0, ZSNOWDZ(ji))
	377	ELSE IF (locjj .LT. 3) THEN
	378	snow_remove(ji)=0 ! masse de neige des couches qui disparaissent totalement
	379	DO jj=3,locjj+1,-1
	380	snow_remove(ji)=snow_remove(ji)+snowdz(ji,jj)*snowrho(ji,jj)
	381	snowdz(ji,jj)=0
	382	ENDDO
	383	ZSNOWMELT_XS(ji) = (snowmelt_from_maxmass(ji) - snow_remove(ji))/snowrho(ji,locjj)
	384	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWMELT_XS(ji)
	385	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
	386	ENDIF
	387	ENDIF
	388
	389	ENDDO
	390
[3269]	391	!2.10 Calculate snow grain size using the updated thermal gradient
[2650]	392
[2591]	393	CALL explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
[2650]	394
[3269]	395	!2.11 Update snow heat
[2223]	396	! Update the heat content (variable stored each time step)
	397	! using current snow temperature and liquid water content:
	398	!
	399	! First, make check to make sure heat content not too large
	400	! (this can result due to signifigant heating of thin snowpacks):
	401	! add any excess heat to ground flux:
	402	!
	403	DO ji=1,kjpindex
	404	DO jj=1,nsnow
[2591]	405	ZLIQHEATXS(ji) = MAX(0.0, snowliq(ji,jj)ph2o - 0.10snowdz(ji,jj)snowrho(ji,jj))chalfu0/dt_sechiba
	406	snowliq(ji,jj) = snowliq(ji,jj) - ZLIQHEATXS(ji)dt_sechiba/(ph2ochalfu0)
[2223]	407	snowliq(ji,jj) = MAX(0.0, snowliq(ji,jj))
	408	grndflux(ji) = grndflux(ji) + ZLIQHEATXS(ji)
[2650]	409	ENDDO
	410	ENDDO
	411
[2223]	412	snow(:) = 0.0
[2650]	413	DO ji=1,kjpindex !domain size
[2223]	414	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
[2650]	415	ENDDO
	416
[2223]	417	DO ji = 1, kjpindex
	418	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
[2650]	419	ENDDO
[4544]	420	snowheattot_end(:)=SUM(snowheat(:,:),2)
	421	del_snowheattot(:)=snowheattot_end(:)-snowheattot_begin(:)
[2650]	422
	423	!3. on land ice (using the default ORCHIDEE snow scheme)
	424
	425	DO ji = 1,kjpindex
	426
[2223]	427	!! 3.1. It is snowing
[2650]	428
[2223]	429	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
	430	& frac_nobio(ji,iice)*precip_rain(ji)
[2650]	431
[4575]	432	!! 3.2. Sublimation - was calculated before, it cannot give us negative snow_nobio (see 2.9).
[2650]	433
[2223]	434	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
[2650]	435
[2223]	436	!! 3.3. ice melt only for continental ice fraction
[2650]	437
[2223]	438	snowmelt_tmp(ji) = zero
	439	IF (temp_sol_new_old(ji) .GT. tp_00) THEN
[2650]	440
[2223]	441	!! 3.3.1 If there is snow on the ice-fraction it can melt
[2650]	442
[2223]	443	snowmelt_tmp(ji) = frac_nobio(ji,iice)(temp_sol_new_old(ji) - tp_00) soilcap(ji) / chalfu0
[2650]	444
[2223]	445	IF ( snowmelt_tmp(ji) .GT. snow_nobio(ji,iice) ) THEN
	446	snowmelt_tmp(ji) = MAX( 0., snow_nobio(ji,iice))
[2650]	447	ENDIF
[2223]	448	snowmelt_ice(ji) = snowmelt_ice(ji) + snowmelt_tmp(ji)
	449	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp(ji)
[2650]	450
	451	ENDIF
[2223]	452
	453	!! Ice melt only if there is more than a given mass : maxmass_snow, i.e. only weight melts glaciers !
[2650]	454
[2223]	455	IF ( snow_nobio(ji,iice) .GE. maxmass_snow ) THEN
	456	icemelt(ji) = snow_nobio(ji,iice) - maxmass_snow
	457	snow_nobio(ji,iice) = maxmass_snow
[2650]	458	ENDIF
	459
[2223]	460	END DO
[2650]	461
	462
[2223]	463	!! 4. On other surface types - not done yet
[2650]	464
[2223]	465	IF ( nnobio .GT. 1 ) THEN
[3006]	466	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
	467	WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
	468	CALL ipslerr_p(3,'explicitsnow_main','Cannot treat snow on these surface types.','','')
[2650]	469	ENDIF
	470
[2223]	471	!! 5. Computes snow age on land and land ice (for albedo)
[2650]	472
[2223]	473	DO ji = 1, kjpindex
[2650]	474
[2223]	475	!! 5.1. Snow age on land
[2650]	476
[2223]	477	IF (snow(ji) .LE. zero) THEN
	478	snow_age(ji) = zero
	479	ELSE
[2591]	480	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
[2223]	481	& * EXP(-precip_snow(ji) / snow_trans)
[2650]	482	ENDIF
	483
[2223]	484	!! 5.2. Snow age on land ice
[2650]	485
[2223]	486	!! Age of snow on ice: a little bit different because in cold regions, we really
	487	!! cannot negect the effect of cold temperatures on snow metamorphism any more.
[2650]	488
[2223]	489	IF (snow_nobio(ji,iice) .LE. zero) THEN
	490	snow_nobio_age(ji,iice) = zero
	491	ELSE
[2650]	492
[2223]	493	d_age(ji) = ( snow_nobio_age(ji,iice) + &
[2591]	494	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
[2223]	495	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
	496	IF (d_age(ji) .GT. 0. ) THEN
	497	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
	498	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
	499	d_age(ji) = d_age(ji) / (un+xx(ji))
[2650]	500	ENDIF
[2223]	501	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
[2650]	502
	503	ENDIF
	504
	505	ENDDO
	506
	507
[2223]	508	!! 6. Check the snow on land
	509	DO ji=1,kjpindex
	510	IF (snow(ji) .EQ. 0) THEN
	511	snowrho(ji,:)=50.0
	512	snowgrain(ji,:)=0.0
	513	snowdz(ji,:)=0.0
	514	snowliq(ji,:)=0.0
[2650]	515	ENDIF
	516	ENDDO
	517
	518
[2223]	519	! Snow melt only if there is more than a given mass : maxmass_snow
	520	! Here I suggest to remove the snow based on a certain threshold of snow
	521	! depth instead of snow mass because it is quite difficult for
	522	! explictsnow module to calculate other snow properties following the
	523	! removal of snow mass
	524	! to define the threshold of snow depth based on old snow density (330
	525	! kg/m3)
	526	! maxmass_snowdepth=maxmass_snow/sn_dens
[3521]	527	! snowmelt_from_maxmass(:) = 0.0
[2650]	528
[2223]	529	!! 7. compute total melt
[2650]	530
[2223]	531	DO ji=1,kjpindex
	532	tot_melt(ji) = icemelt(ji) + snowmelt(ji) + snowmelt_ice(ji) + snowmelt_from_maxmass(ji)
[2650]	533	ENDDO
[2348]	534	IF (printlev>=3) WRITE(numout,*) 'explicitsnow_main done'
[7266]	535
	536
	537	!! 8. Recalculate the new snow discretization
	538	! This is done to make sure that all levels have snow
	539	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
	540	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
[2650]	541
[4881]	542
[3410]	543	! Write output variables
[4881]	544
	545	! Add XIOS default value where no snow
	546	DO ji=1,kjpindex
	547	IF (snow(ji) .GT. zero) THEN
[5091]	548	! Output for snowliq and snowliqtot change unit from m to kg/m2 and are multiplied by contfrac
	549	! to take into acount the portion form land fraction only
	550	snowliq_diag(ji,:) = snowliq(ji,:) * 1000 * contfrac(ji)
	551	snowliqtot_diag(ji) = SUM(snowliq_diag(ji,:)) * 1000 * contfrac(ji)
[4881]	552	snowrho_diag(ji,:) = snowrho(ji,:)
	553	snowheat_diag(ji,:) = snowheat(ji,:)
	554	snowgrain_diag(ji,:) = snowgrain(ji,:)
	555	ELSE
	556	snowliq_diag(ji,:) = xios_default_val
	557	snowliqtot_diag(ji) = xios_default_val
	558	snowrho_diag(ji,:) = xios_default_val
	559	snowheat_diag(ji,:) = xios_default_val
	560	snowgrain_diag(ji,:) = xios_default_val
	561	END IF
	562	END DO
	563
	564	CALL xios_orchidee_send_field("snowliq",snowliq_diag)
	565	CALL xios_orchidee_send_field("snowliqtot", snowliqtot_diag)
	566	CALL xios_orchidee_send_field("snowrho",snowrho_diag)
	567	CALL xios_orchidee_send_field("snowheat",snowheat_diag)
	568	CALL xios_orchidee_send_field("snowgrain",snowgrain_diag)
[3521]	569	CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass)
	570	CALL xios_orchidee_send_field("soilcap",soilcap)
[4544]	571	CALL xios_orchidee_send_field("del_snowheattot",del_snowheattot)
[3410]	572
[2223]	573	END SUBROUTINE explicitsnow_main
[2650]	574
[2223]	575
[2650]	576	!================================================================================================================================
	577	!! SUBROUTINE : explicitsnow_finalize
	578	!!
	579	!>\BRIEF Write variables for explictsnow module to restart file
	580	!!
	581	!! DESCRIPTION : Write variables for explictsnow module to restart file
	582	!!
	583	!! \n
	584	!_
	585	!================================================================================================================================
	586	SUBROUTINE explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, &
[3059]	587	snowtemp, snowdz, snowheat, snowgrain)
[2650]	588
[2223]	589	!! 0.1 Input variables
[2650]	590	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
	591	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
	592	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
[3410]	593	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
[2650]	594	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
	595	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness
[3410]	596	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content (J/m^2)
	597	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grainsize (m)
[2650]	598
	599
	600	!! 1. Write to restart file
	601	CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g)
	602	CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g)
	603	CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g)
	604	CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g)
	605	CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g)
	606
	607	END SUBROUTINE explicitsnow_finalize
	608
[2223]	609
[2650]	610	!================================================================================================================================
	611	!! SUBROUTINE : explicitsnow_grain
	612	!!
	613	!>\BRIEF Compute evolution of snow grain size
	614	!!
	615	!! DESCRIPTION :
	616	!!
	617	!! RECENT CHANGE(S) : None
	618	!!
	619	!! MAIN OUTPUT VARIABLE(S): None
	620	!!
	621	!! REFERENCE(S) : R. Jordan (1991)
	622	!!
	623	!! FLOWCHART : None
	624	!! \n
	625	!_
	626	!================================================================================================================================
[2223]	627
	628
[2650]	629	SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
[2223]	630
[2650]	631	!! 0.1 Input variables
	632	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
	633	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content
	634	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth (m)
	635	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
	636	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
	637	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pb !! Surface pressure (hpa)
[2223]	638
[2650]	639	!! 0.2 Output variables
[2223]	640
[2650]	641	!! 0.3 Modified variables
[2223]	642
[3410]	643	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
[2223]	644
[2650]	645	!! 0.4 Local variables
	646	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowdz,zdz,ztheta
	647	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: ztemp,zdiff,ztgrad,zthetaa,zfrac,&
	648	zexpo,zckt_liq,zckt_ice,zckt
	649	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrhomin,zgrainmin
	650	INTEGER(i_std) :: ji,jj
[2223]	651
[2650]	652	!! 0.5 Local parameters
	653	REAL(r_std), PARAMETER :: ztheta_crit = 0.02 !! m3 m-3
	654	REAL(r_std), PARAMETER :: zc1_ice = 8.047E+9 !! kg m-3 K
	655	REAL(r_std), PARAMETER :: zc1_liq = 5.726E+8 !! kg m-3 K
	656	REAL(r_std), PARAMETER :: zdeos = 0.92E-4 !! effective diffusion
	657	!! coef for water vapor in snow
	658	!! at 0C and 1000 mb (m2 s-1)
	659	REAL(r_std), PARAMETER :: zg1 = 5.0E-7 !! m4 kg-1
	660	REAL(r_std), PARAMETER :: zg2 = 4.0E-12 !! m2 s-1
	661	REAL(r_std), PARAMETER :: ztheta_w = 0.05 !! m3 m-3
	662	REAL(r_std), PARAMETER :: ztheta_crit_w = 0.14 !! m3 m-3
	663	REAL(r_std), PARAMETER :: zdzmin = 0.01 !! m : minimum thickness
	664	!! for thermal gradient evaluation:
	665	!! to prevent excessive gradients
	666	!! for vanishingly thin snowpacks.
	667	REAL(r_std), PARAMETER :: xp00=1.E5
	668	!! 1. initialize
	669
	670	DO ji=1,kjpindex
[2223]	671
	672
[2650]	673	zsnowdz(ji,:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
[2223]	674
[2650]	675	DO jj=1,nsnow-1
	676	zdz(ji,jj) = zsnowdz(ji,jj) + zsnowdz(ji,jj+1)
	677	ENDDO
	678	zdz(ji,nsnow) = zsnowdz(ji,nsnow)
[2223]	679
[2650]	680	! compute interface average volumetric water content (m3 m-3):
	681	! first, layer avg VWC:
	682	!
	683	ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:))
[2223]	684
[2650]	685	! at interfaces:
	686	zthetaa(ji,0) = ztheta(ji,1)
	687	DO jj=1,nsnow-1
	688	zthetaa(ji,jj) = (zsnowdz(ji,jj) *ztheta(ji,jj) + &
	689	zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj)
	690	ENDDO
	691	zthetaa(ji,nsnow) = ztheta(ji,nsnow)
	692	! compute interface average temperatures (K):
	693	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	694	!
	695	ztemp(ji,0) = snowtemp(ji,1)
	696	DO jj=1,nsnow-1
	697	ztemp(ji,jj) = (zsnowdz(ji,jj) *snowtemp(ji,jj) + &
	698	zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj)
	699	ENDDO
	700	ztemp(ji,nsnow) = snowtemp(ji,nsnow)
	701	!
	702	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	703	! compute variation of saturation vapor pressure with temperature
	704	! for solid and liquid phases:
	705	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	706	zexpo(ji,:) = chalsu0/(xrv*ztemp(ji,:))
	707	zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
	708	!
	709	zexpo(ji,:) = chalev0/(xrv*ztemp(ji,:))
	710	zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
	711	!
	712	! compute the weighted ice/liquid total variation (N m-2 K):
	713	!
	714	zfrac(ji,:) = MIN(1.0, zthetaa(ji,:)/ztheta_crit)
	715	zckt(ji,:) = zfrac(ji,:)zckt_liq(ji,:) + (1.0 - zfrac(ji,:))zckt_ice(ji,:)
[2223]	716
[2650]	717	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	718	! Compute effective diffusion coefficient (m2 s-1):
	719	! -diffudivity relative to a reference diffusion at 1000 mb and freezing point
	720	! multiplied by phase energy coefficient
	721	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	722	!
	723	DO jj=0,nsnow
	724	zdiff(ji,jj) = zdeos(xp00/(pb(ji)100.))((ztemp(ji,jj)/tp_00)6)zckt(ji,jj)
	725	ENDDO
[2223]	726
[2650]	727	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	728	! Temperature gradient (K m-1):
[2223]	729
[2650]	730	ztgrad(ji,0) = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal
	731	DO jj=1,nsnow-1
	732	ztgrad(ji,jj) = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj))
	733	ENDDO
	734	!
	735	! assume at base of snow, temperature is in equilibrium with soil
	736	! (but obviously must be at or below freezing point!)
	737	!
	738	ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow))
	739	! prognostic grain size (m) equation:
	740	!-------------------------------------------------------------------
	741	! first compute the minimum grain size (m):
	742	!
	743	zrhomin(ji,:) = xrhosmin
	744	zgrainmin(ji,:) = snow3lgrain_1d(zrhomin(ji,:))
[2223]	745
[2650]	746	! dry snow:
	747	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	748	!
	749	DO jj=1,nsnow
[2223]	750
[2650]	751	IF(ztheta(ji,jj) == 0.0) THEN
[2223]	752
[2650]	753	! only allow growth due to vapor flux INTO layer:
	754	! aab add sublimation(only condensation) as upper BC...?
[2223]	755
[2650]	756	snowgrain(ji,jj) = snowgrain(ji,jj) + &
	757	(dt_sechibazg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
	758	( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) - &
	759	zdiff(ji,jj) *MIN(0.0,ztgrad(ji,jj)) )
	760	ELSE
	761
	762	! wet snow
	763	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	764	!
	765	snowgrain(ji,jj) = snowgrain(ji,jj) + &
	766	(dt_sechibazg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
	767	MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w)
	768	END IF
[2223]	769
[2650]	770	ENDDO
[2223]	771
	772
[2650]	773	ENDDO
[2223]	774
	775
[2650]	776	END SUBROUTINE explicitsnow_grain
[2223]	777
[2650]	778	!!
	779	!================================================================================================================================
	780	!! SUBROUTINE : explicitsnow_compactn
	781	!!
	782	!>\BRIEF Compute Compaction/Settling
	783	!!
	784	!! DESCRIPTION :
	785	!! Snow compaction due to overburden and settling.
	786	!! Mass is unchanged: layer thickness is reduced
	787	!! in proportion to density increases. Method
	788	!! of Anderson (1976): see Loth and Graf, 1993,
	789	!! J. of Geophys. Res., 98, 10,451-10,464.
	790	!!
	791	!! RECENT CHANGE(S) : None
	792	!!
	793	!! MAIN OUTPUT VARIABLE(S): None
	794	!!
	795	!! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976)
	796	!!
	797	!! FLOWCHART : None
	798	!! \n
	799	!_
	800	!================================================================================================================================
[2223]	801
	802
[2650]	803	SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
[2223]	804
[2650]	805	!! 0.1 Input variables
[2223]	806
[2650]	807	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
	808	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
[2223]	809
[2650]	810	!! 0.2 Output variables
[2223]	811
[2650]	812	!! 0.3 Modified variables
[2223]	813
[3410]	814	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
[2650]	815	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth
[2223]	816
[2650]	817	!! 0.4 Local variables
[2223]	818
[2650]	819	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle
	820	REAL(r_std),DIMENSION(kjpindex) :: zsmassc !! cummulative snow mass (kg/m2)
	821	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
	822	INTEGER(i_std) :: ji,jj
[2223]	823
[2650]	824	!! 1. initialize
[2223]	825
[2650]	826	zsnowrho2 = snowrho
	827	zsettle(:,:) = ZSNOWCMPCT_ACM
	828	zviscocity(:,:) = ZSNOWCMPCT_V0
[2223]	829
[2650]	830	!! 2. Calculating Cumulative snow mass (kg/m2):
[2223]	831
[2650]	832	DO ji=1, kjpindex
[2223]	833
	834
[2650]	835	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
[2223]	836
[2650]	837	zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:)
[2223]	838
[2650]	839	zsmassc (ji)= 0.0
[2223]	840
[2650]	841	DO jj=1,nsnow
	842	zsmass(ji,jj) = zsmassc(ji) + zwsnowdz(ji,jj)
	843	zsmassc(ji) = zsmassc(ji) + zwsnowdz(ji,jj)
	844	ENDDO
[2223]	845
	846
[2650]	847	!! 3. Computing compaction/Settling
	848	! ----------------------
	849	! Compaction/settling if density below upper limit
	850	! (compaction is generally quite small above ~ 500 kg m-3):
	851	!
	852	DO jj=1,nsnow
	853	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
	854
	855	!
	856	! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile)
	857	!
	858
	859	zsettle(ji,jj) = ZSNOWCMPCT_ACM*EXP( &
	860	-ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj))) &
	861	-ZSNOWCMPCT_CCM*MAX(0.0, &
	862	snowrho(ji,jj)-ZSNOWCMPCT_RHOD))
	863	!
	864	! Snow viscocity:
	865	!
	866	zviscocity(ji,jj) = ZSNOWCMPCT_V0EXP( ZSNOWCMPCT_VT(tp_00-MIN(tp_00,snowtemp(ji,jj))) + &
	867	ZSNOWCMPCT_VR*snowrho(ji,jj) )
[2223]	868
[2650]	869	! Calculate snow density: compaction from weight/over-burden
	870	! Anderson 1976 method:
	871	zsnowrho2(ji,jj) = snowrho(ji,jj) + snowrho(ji,jj)dt_sechiba( &
	872	(cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) &
	873	+ zsettle(ji,jj) )
	874	! Conserve mass by decreasing grid thicknesses in response
	875	! to density increases
	876	!
	877	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
[2223]	878
[2650]	879	ENDIF
	880	ENDDO
[2223]	881
[2650]	882	! Update density (kg m-3):
	883	snowrho(ji,:) = zsnowrho2(ji,:)
[2223]	884
[2650]	885	ENDIF
[2223]	886
[2650]	887	ENDDO
[2223]	888
[2650]	889	END SUBROUTINE explicitsnow_compactn
[2223]	890
[2650]	891	!!
	892	!================================================================================================================================
	893	!! SUBROUTINE : explicitsnow_transf
	894	!!
	895	!>\BRIEF Computing snow mass and heat redistribution due to grid thickness configuration resetting
	896	!!
	897	!! DESCRIPTION : Snow mass and heat redistibution due to grid thickness
	898	!! configuration resetting. Total mass and heat content
	899	!! of the overall snowpack unchanged/conserved within this routine.
	900	!! RECENT CHANGE(S) : None
	901	!!
	902	!! MAIN OUTPUT VARIABLE(S): None
	903	!!
	904	!! REFERENCE(S) :
	905	!!
	906	!! FLOWCHART : None
	907	!! \n
	908	!_
	909	!================================================================================================================================
[2223]	910
[2650]	911	SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
[2223]	912
[2650]	913	!! 0.1 Input variables
[2223]	914
[2650]	915	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
	916	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz_old !! Snow depth at the previous time step
[2223]	917
[2650]	918	!! 0.2 Output variables
[2223]	919
[2650]	920	!! 0.3 Modified variables
[2223]	921
[3410]	922	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
	923	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
[2650]	924	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth (m)
	925	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
	926
	927	!! 0.4 Local varibles
	928
[4305]	929	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzo
	930	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzn
[2650]	931	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowddz
	932	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zdelta
	933	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrhon,zsnowheatn,zsnowgrainn
	934	REAL(r_std),DIMENSION(kjpindex) :: zsnowmix_delta
	935	REAL(r_std),DIMENSION(kjpindex) :: zsumheat, zsumswe, zsumgrain
	936	INTEGER(i_std),DIMENSION(nsnow,2) :: locflag
	937	REAL(r_std) :: psnow
	938	INTEGER(i_std) :: ji,jj,jjj
[2223]	939
[2650]	940	! Initialization
	941	zsumheat(:) = 0.0
	942	zsumswe(:) = 0.0
	943	zsumgrain(:) = 0.0
	944	zsnowmix_delta(:) = 0.0
	945	locflag(:,:) = 0
[2223]	946
[2650]	947	DO ji=1, kjpindex
	948
[2223]	949
[2650]	950	psnow = SUM(snowdz(ji,:))
[2223]	951
[4305]	952	IF (psnow .GE. xsnowcritd .AND. snowdz_old(ji,1) .NE. 0 .AND. snowdz_old(ji,2) .NE. 0 .AND. snowdz_old(ji,3) .NE. 0) THEN
[2223]	953	!
[4305]	954	zsnowzo(ji,0) = 0.
	955	zsnowzn(ji,0) = 0.
[2223]	956	zsnowzo(ji,1) = snowdz_old(ji,1)
	957	zsnowzn(ji,1) = snowdz(ji,1)
	958	!
[3085]	959	DO jj=2,nsnow
[2223]	960	zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj)
	961	zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj)
[2650]	962	ENDDO
[2223]	963
[4305]	964	DO jj=1,nsnow
[2223]	965
[4305]	966	IF (jj .EQ. 1) THEN
	967	locflag(jj,1)=1
	968	ELSE
	969	DO jjj=nsnow,1,-1
[2223]	970
[4305]	971	!upper bound of the snow layer
	972	IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN
	973	locflag(jj,1) = jjj
	974	ENDIF
	975	ENDDO
	976	ENDIF
[2223]	977
[4305]	978	IF (jj .EQ. nsnow) THEN
	979	locflag(jj,2)=nsnow
	980	ELSE
	981	DO jjj=nsnow,1,-1
[2223]	982
[4305]	983	!lower bound of the snow layer
	984	IF (zsnowzn(ji,jj) .LE. zsnowzo(ji,jjj)) THEN
	985	locflag(jj,2) = jjj
	986	ENDIF
	987	ENDDO
	988	ENDIF
[2223]	989
	990	!to interpolate
	991	! when heavy snow occurred
	992	IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN
	993
	994	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1))
	995
	996	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1))
	997
	998	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1))
	999
	1000	ELSE
	1001
	1002	!snow density
	1003	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * &
[2650]	1004	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
	1005	snowrho(ji,locflag(jj,2)) * &
	1006	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
[2223]	1007
	1008	DO jjj=locflag(jj,1),locflag(jj,2)-1
	1009	zsnowrhon(ji,jj)=zsnowrhon(ji,jj) + &
[2650]	1010	(jjj-locflag(jj,1))snowrho(ji,jjj)snowdz_old(ji,jjj)
	1011	ENDDO
[2223]	1012
	1013	zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj)
	1014
	1015
	1016	!snow heat
[4305]	1017	IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN
[2223]	1018
[4305]	1019	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * &
[2650]	1020	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + &
	1021	snowheat(ji,locflag(jj,2)) * &
	1022	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2))
[4305]	1023	ELSE
[4316]	1024	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2)) * &
	1025	(zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1)))
[4305]	1026	ENDIF
[2223]	1027
	1028	DO jjj=locflag(jj,1),locflag(jj,2)-1
	1029	zsnowheatn(ji,jj)=zsnowheatn(ji,jj) + &
[2650]	1030	(jjj-locflag(jj,1))*snowheat(ji,jjj)
	1031	ENDDO
[2223]	1032
	1033
	1034	!snow grain
	1035	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * &
[2650]	1036	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
	1037	snowgrain(ji,locflag(jj,2)) * &
	1038	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
[2223]	1039
	1040	DO jjj=locflag(jj,1),locflag(jj,2)-1
	1041	zsnowgrainn(ji,jj)=zsnowgrainn(ji,jj) + &
[2650]	1042	(jjj-locflag(jj,1))snowgrain(ji,jjj)snowdz_old(ji,jjj)
	1043	ENDDO
[2223]	1044
	1045	zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj)
	1046
[2650]	1047	ENDIF
	1048
	1049	ENDDO
[4305]	1050
[2223]	1051	snowrho(ji,:) = zsnowrhon(ji,:)
	1052	snowheat(ji,:) = zsnowheatn(ji,:)
	1053	snowgrain(ji,:) = zsnowgrainn(ji,:)
	1054
[4305]	1055	ENDIF
[2223]	1056
[2650]	1057	! Vanishing or very thin snowpack check:
	1058	! -----------------------------------------
	1059	!
	1060	! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
	1061	! this avoids problems related to heat and mass exchange for
	1062	! thin layers during heavy snowfall or signifigant melt: one
	1063	! new/old layer can exceed the thickness of several old/new layers.
	1064	! Therefore, mix (conservative):
	1065	!
	1066	! modified by Tao Wang
[4305]	1067	IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) &
[2223]	1068	& .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN
[2650]	1069	zsumheat(ji) = SUM(snowheat(ji,:))
	1070	zsumswe(ji) = SUM(snowrho(ji,:)*snowdz_old(ji,:))
	1071	zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:))
	1072	zsnowmix_delta(ji) = 1.0
	1073	DO jj=1,nsnow
	1074	zsnowheatn(ji,jj) = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow)
	1075	snowdz(ji,jj) = zsnowmix_delta(ji)*(psnow/nsnow)
	1076	zsnowrhon(ji,jj) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow)
	1077	zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow)
	1078	ENDDO
[2223]	1079	! Update mass (density and thickness), heat and grain size:
	1080	! ------------------------------------------------------------
	1081	!
	1082	snowrho(ji,:) = zsnowrhon(ji,:)
	1083	snowheat(ji,:) = zsnowheatn(ji,:)
	1084	snowgrain(ji,:) = zsnowgrainn(ji,:)
	1085
[2650]	1086	ENDIF
[2223]	1087
[2650]	1088	ENDDO
[2223]	1089
	1090
[2650]	1091	END SUBROUTINE explicitsnow_transf
[2223]	1092
[2650]	1093
	1094	!!
	1095	!================================================================================================================================
	1096	!! SUBROUTINE : explicitsnow_fall
	1097	!!
	1098	!>\BRIEF Computes snowfall
	1099	!!
[4282]	1100	!! DESCRIPTION : Computes snowfall
[2650]	1101	!routine.
	1102	!! RECENT CHANGE(S) : None
	1103	!!
	1104	!! MAIN OUTPUT VARIABLE(S): None
	1105	!!
	1106	!! REFERENCE(S) :
	1107	!!
	1108	!! FLOWCHART : None
	1109	!! \n
	1110	!_
	1111	!================================================================================================================================
	1112
	1113	SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
	1114	& snowheat,snowgrain,snowtemp,psnowhmass)
[2223]	1115
[2650]	1116	!! 0.1 Input variables
[2223]	1117
[2650]	1118	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
	1119	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba)
	1120	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
	1121	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! Horizontal wind speed
	1122	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
	1123	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio
	1124	!! 0.2 Output variables
[2223]	1125
[2650]	1126	REAL(r_std), DIMENSION(kjpindex),INTENT(out) :: psnowhmass !! Heat content of snowfall (J/m2)
[2223]	1127
[2650]	1128	!! 0.3 Modified variables
	1129
	1130	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density profile (kg/m3)
	1131	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow layer thickness profile (m)
	1132	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
[3410]	1133	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
[2223]	1134
[2650]	1135	!! 0.4 Local variables
[2223]	1136
[2650]	1137	REAL(r_std), DIMENSION(kjpindex) :: rhosnew !! Snowfall density
	1138	REAL(r_std), DIMENSION(kjpindex) :: dsnowfall !! Snowfall thickness (m)
	1139	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowdz_old
	1140	REAL(r_std), DIMENSION(kjpindex) :: snow_depth,snow_depth_old,newgrain
	1141	REAL(r_std) :: snowfall_delta,speed
	1142	INTEGER(i_std) :: ji,jj
[2223]	1143
[2650]	1144	!! 1. initialize the variables
[2223]	1145
[2650]	1146	snowdz_old = snowdz
	1147	DO ji=1,kjpindex
	1148	snow_depth(ji) = SUM(snowdz(ji,:))
	1149	ENDDO
[2223]	1150
[2650]	1151	snow_depth_old = snow_depth
	1152
	1153	snowfall_delta = 0.0
[2223]	1154
[2650]	1155	!! 2. incorporate snowfall into snowpack
	1156	DO ji = 1, kjpindex
	1157
[2223]	1158	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
[2650]	1159
	1160	! new snow fall on snowpack
	1161	! NOTE: when the surface temperature is zero, it means that snowfall has no
	1162	! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug)
	1163	psnowhmass(ji) = 0.0
	1164	IF ( (precip_snow(ji) .GT. 0.0) ) THEN
[2223]	1165
[2650]	1166	!calculate
[2223]	1167
[2650]	1168	psnowhmass(ji) = precip_snow(ji)(un-totfrac_nobio(ji)) &
	1169	(xci*(snowtemp(ji,1)-tp_00)-chalfu0)
	1170
	1171	! Snowfall density: Following CROCUS (Pahaut 1976)
	1172	!
	1173	rhosnew(ji) = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00)+ &
	1174	snowfall_c_sn*SQRT(speed))
[2223]	1175
[2650]	1176	! Augment total pack depth:
[3085]	1177	! Note that dsnowfall(ji) can be equal to zero if totfrac_nobio(ji)=1
[2650]	1178	dsnowfall(ji) = (precip_snow(ji)*(un-totfrac_nobio(ji)))/rhosnew(ji) !snowfall thickness (m)
[2223]	1179
[2650]	1180	snow_depth(ji) = snow_depth(ji) + dsnowfall(ji)
[2223]	1181
	1182
[2650]	1183	! Fresh snowfall changes the snowpack density and liquid content in uppermost layer
[2223]	1184
[3085]	1185	IF (dsnowfall(ji) .GT. zero) THEN
[2650]	1186	snowrho(ji,1) = (snowdz(ji,1)snowrho(ji,1) + dsnowfall(ji)rhosnew(ji))/ &
	1187	(snowdz(ji,1)+dsnowfall(ji))
[2223]	1188
[3085]	1189	snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji)
[2223]	1190
[3085]	1191	! Add energy of snowfall to snowpack:
	1192	! Update heat content (J/m2) (therefore the snow temperature
	1193	! and liquid content):
	1194	!
	1195	snowheat(ji,1) = snowheat(ji,1) + psnowhmass(ji)
	1196	!
	1197	! Incorporate snowfall grain size:
	1198	!
	1199	newgrain(ji) = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji)))
	1200
	1201	snowgrain(ji,1) = (snowdz_old(ji,1)snowgrain(ji,1) + dsnowfall(ji)newgrain(ji))/ &
	1202	snowdz(ji,1)
	1203	ENDIF
	1204	ELSE
	1205	dsnowfall(ji) = 0.
[2650]	1206	ENDIF
[2223]	1207
[2650]	1208	! new snow fall on snow free surface.
	1209	! we use the linearization for the new snow fall on snow-free ground
[2223]	1210
[3085]	1211	IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN
[2223]	1212
[2650]	1213	snowfall_delta = 1.0
	1214
	1215	DO jj=1,nsnow
[2223]	1216
[2650]	1217	snowdz(ji,jj) = snowfall_delta*(dsnowfall(ji)/nsnow) + &
	1218	(1.0-snowfall_delta)*snowdz(ji,jj)
[2223]	1219
[2650]	1220	snowheat(ji,jj) = snowfall_delta*(psnowhmass(ji)/nsnow) + &
	1221	(1.0-snowfall_delta)*snowheat(ji,jj)
[2223]	1222
[2650]	1223	snowrho(ji,jj) = snowfall_delta*rhosnew(ji) + &
	1224	(1.0-snowfall_delta)*snowrho(ji,jj)
[2223]	1225
[2650]	1226	snowgrain(ji,jj)= snowfall_delta*newgrain(ji) + &
	1227	(1.0-snowfall_delta)*snowgrain(ji,jj)
[2223]	1228
[2650]	1229	ENDDO
[2223]	1230
	1231
[2650]	1232	ENDIF
[2223]	1233
	1234
[2650]	1235	ENDDO
[2223]	1236
[2650]	1237	END SUBROUTINE explicitsnow_fall
[2223]	1238
[2650]	1239	!!
	1240	!================================================================================================================================
	1241	!! SUBROUTINE : explicitsnow_gone
	1242	!!
	1243	!>\BRIEF Check whether snow is gone
	1244	!!
	1245	!! DESCRIPTION : If so, set thickness (and therefore mass and heat) and liquid
	1246	!! content to zero, and adjust fluxes of water, evaporation and
	1247	!! heat into underlying surface.
	1248	!! RECENT CHANGE(S) : None
	1249	!!
	1250	!! MAIN OUTPUT VARIABLE(S): None
	1251	!!
	1252	!! REFERENCE(S) :
	1253	!!
	1254	!! FLOWCHART : None
	1255	!! \n
	1256	!_
	1257	!================================================================================================================================
[2223]	1258
[2650]	1259	SUBROUTINE explicitsnow_gone(kjpindex,pgflux,&
	1260	snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
[2223]	1261
[2650]	1262	!! 0.1 Input variables
[2223]	1263
[2650]	1264	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
	1265	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pgflux !! Net energy into snow pack(w/m2)
[3410]	1266	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowheat !! Snow heat content (J/m^2)
[2223]	1267
[2650]	1268	!! 0.2 Output variables
[2223]	1269
[2650]	1270	!! 0.3 Modified variables
[2223]	1271
[2650]	1272	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
	1273	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
[3410]	1274	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
[2650]	1275	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
	1276	REAL(r_std),DIMENSION(kjpindex), INTENT(inout) :: grndflux !! Soil/snow interface heat flux (W/m2)
	1277	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: snowmelt !! Snow melt
	1278	REAL(r_std),DIMENSION(kjpindex) :: thrufal !! Water leaving snowpack(kg/m2/s)
[2223]	1279
[2650]	1280	!! 0.4 Local variables
[2223]	1281
[2650]	1282	INTEGER(i_std) :: ji,jj
	1283	REAL(r_std),DIMENSION(kjpindex) :: snowgone_delta
	1284	REAL(r_std),DIMENSION (kjpindex) :: totsnowheat !!snow heat content at each layer
	1285	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
[2223]	1286
[2650]	1287	! first caculate total snowpack snow heat content
	1288	snowgone_delta(:) = un
	1289	thrufal(:)=0.0
	1290	snowmelt(:)=0
	1291	totsnowheat(:) = SUM(snowheat(:,:),2)
	1292
	1293	DO ji = 1, kjpindex
[2223]	1294
[3772]	1295	IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN
[2223]	1296
[2650]	1297	IF( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN
[2223]	1298
[2650]	1299	grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba)
[2223]	1300
[2650]	1301	thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:))
[2223]	1302
[2650]	1303	snowgone_delta(ji) = 0.0
[2223]	1304
[2650]	1305	snowmelt(ji) = snowmelt(ji)+thrufal(ji)
[2223]	1306
[2650]	1307	ENDIF
[2223]	1308
[2650]	1309	! update of snow state (either still present or not)
[2223]	1310
[2650]	1311	DO jj=1,nsnow
	1312	snowdz(ji,jj) = snowdz(ji,jj) *snowgone_delta(ji)
	1313	snowliq(ji,jj) = snowliq(ji,jj) *snowgone_delta(ji)
	1314	snowtemp(ji,jj) = (1.0-snowgone_delta(ji))tp_00 + snowtemp(ji,jj)snowgone_delta(ji)
	1315	ENDDO
[2223]	1316
[3772]	1317	ELSE
	1318	snowdz(ji,:)=0.
	1319	snowliq(ji,:)=0.
	1320	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
[2650]	1321	ENDIF
	1322	ENDDO
[2223]	1323
[2650]	1324	END SUBROUTINE explicitsnow_gone
[2223]	1325
[2650]	1326	!================================================================================================================================
	1327	!! SUBROUTINE : explicitsnow_melt_refrz
	1328	!!
	1329	!>\BRIEF Computes snow melt and refreezing processes within snowpack
	1330	!!
	1331	!! DESCRIPTION :
	1332	!! RECENT CHANGE(S) : None
	1333	!!
	1334	!! MAIN OUTPUT VARIABLE(S): None
	1335	!!
	1336	!! REFERENCE(S) :
	1337	!!
	1338	!! FLOWCHART : None
	1339	!! \n
	1340	!_
	1341	!================================================================================================================================
[2223]	1342
[2650]	1343	SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
	1344	snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
[2223]	1345
[2650]	1346	!! 0.1 Input variables
[2223]	1347
[2650]	1348	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
	1349	REAL(r_std),DIMENSION (kjpindex,nsnow) :: pcapa_snow !! Heat capacity for snow
	1350	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
	1351	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
	1352	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pgflux !! Net energy into snowpack(w/m2)
	1353	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity
[2223]	1354
[2650]	1355	!! 0.2 Output variables
[2223]	1356
[2650]	1357	!! 0.3 Modified variables
[2223]	1358
[2650]	1359	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
	1360	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
[3410]	1361	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow layer density (Kg/m^3)
[2650]	1362	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
	1363	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil
	1364	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: snowmelt !! Snowmelt
[2223]	1365
[2650]	1366	!! 0.4 Local variables
[2223]	1367
[2650]	1368	REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil
	1369	REAL(r_std) :: enerin,melttot,hrain
	1370	REAL(r_std),DIMENSION (nsnow) :: zsnowlwe
	1371	REAL(r_std),DIMENSION (nsnow) :: flowliq
	1372	REAL(r_std),DIMENSION (kjpindex) :: snowmass
	1373	REAL(r_std),DIMENSION (nsnow) :: zphase !! Phase change (from ice to water) (J/m2)
	1374	REAL(r_std),DIMENSION (nsnow) :: zphase2 !! Phase change (from water to ice)
	1375	REAL(r_std),DIMENSION (nsnow) :: zphase3 !! Phase change related with net energy input to snowpack
	1376	REAL(r_std),DIMENSION (nsnow) :: zsnowdz !! Snow layer depth
	1377	REAL(r_std),DIMENSION (nsnow) :: zsnowmelt !! Snow melt (liquid water) (m)
	1378	REAL(r_std),DIMENSION (nsnow) :: zsnowtemp
	1379	REAL(r_std),DIMENSION (nsnow) :: zmeltxs !! Excess melt
	1380	REAL(r_std),DIMENSION (nsnow) :: zwholdmax !! Maximum liquid water holding (m)
	1381	REAL(r_std),DIMENSION (nsnow) :: zcmprsfact !! Compression factor due to densification from melting
	1382	REAL(r_std),DIMENSION (nsnow) :: zscap !! Snow heat capacity (J/m3 K)
	1383	REAL(r_std),DIMENSION (nsnow) :: zsnowliq !! (m)
	1384	REAL(r_std),DIMENSION (nsnow) :: snowtemp_old
	1385	REAL(r_std),DIMENSION (0:nsnow) :: zflowliqt !!(m)
	1386	REAL(r_std) :: zrainfall,zpcpxs
	1387	REAL(r_std) :: ztotwcap
	1388	REAL(r_std),DIMENSION(kjpindex,nsnow) :: snowdz_old,snowliq_old
	1389	INTEGER(i_std) :: jj,ji, iv
	1390	REAL(r_std),DIMENSION(nsnow) :: snowdz_old2
	1391	REAL(r_std),DIMENSION(nsnow) :: zsnowrho
	1392
	1393	!initialize
	1394	snowdz_old = snowdz
	1395	snowliq_old = snowliq
[2223]	1396
[2650]	1397	DO ji = 1, kjpindex
[2223]	1398
	1399
[2650]	1400	snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
[3772]	1401	IF ((snowmass(ji) .GT. min_sechiba)) THEN
[2223]	1402
[2650]	1403	!! 1 snow melting due to positive snowpack snow temperature
	1404
	1405	!! 1.0 total liquid equivalent water content of each snow layer
[2223]	1406
[2650]	1407	zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o
[2223]	1408
[2650]	1409	!! 1.1 phase change (J/m2)
[2223]	1410
[2650]	1411	pcapa_snow(ji,:) = snowrho(ji,:)*xci
[2223]	1412
[2650]	1413	zphase(:) = MIN(pcapa_snow(ji,:)MAX(0.0, snowtemp(ji,:)-tp_00) &
	1414	snowdz(ji,:), &
	1415	MAX(0.0,zsnowlwe(:)-snowliq(ji,:))chalfu0ph2o)
[2223]	1416
[2650]	1417	!! 1.2 update snow liq water and temperature if melting
[2223]	1418
[2650]	1419	zsnowmelt(:) = zphase(:)/(chalfu0*ph2o)
[2223]	1420
[2650]	1421	!! 1.3 cool off the snow layer temperature due to melt
[2223]	1422
[2650]	1423	zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:))
[2223]	1424
[2650]	1425	snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:))
[2223]	1426
[2650]	1427	zmeltxs(:) = (zsnowtemp(:)-snowtemp(ji,:))pcapa_snow(ji,:)snowdz(ji,:)
[2223]	1428
[2650]	1429	!! 1.4 loss of snowpack depth and liquid equivalent water
[2223]	1430
[2650]	1431	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
[2223]	1432
[2650]	1433	zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ &
	1434	(zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:)))
[2223]	1435
[2650]	1436	snowdz(ji,:) = snowdz(ji,:)*zcmprsfact(:)
[2223]	1437
[2650]	1438	snowrho(ji,:) = zsnowlwe(:)*ph2o/snowdz(ji,:)
[2223]	1439
[2650]	1440	snowliq(ji,:) = snowliq(ji,:) + zsnowmelt(:)
[2223]	1441
	1442
[2650]	1443	!! 2 snow refreezing process
	1444	!! 2.1 freeze liquid water in any layer
	1445	zscap(:) = snowrho(ji,:)*xci !J K-1 m-3
	1446	zphase2(:) = MIN(zscap(:)* &
	1447	MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:), &
	1448	snowliq(ji,:)chalfu0ph2o)
[2223]	1449
[2650]	1450	! warm layer and reduce liquid if freezing occurs
	1451	zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
	1452	snowtemp_old(:) = snowtemp(ji,:)
	1453	snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:))
	1454
	1455	! Reduce liquid portion if freezing occurs:
	1456	snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))* &
[2223]	1457	zscap(:)zsnowdz(:)/(chalfu0ph2o) )
[2650]	1458	snowliq(ji,:) = MAX(snowliq(ji,:), 0.0)
	1459	!! 3. thickness change due to snowmelt in excess of holding capacity
	1460	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
	1461	flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:)))
	1462	snowliq(ji,:) = snowliq(ji,:) - flowliq(:)
	1463	snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:)
	1464	! to prevent possible very small negative values (machine
	1465	! prescision as snow vanishes
	1466	snowdz(ji,:) = MAX(0.0, snowdz(ji,:))
[2223]	1467
[2650]	1468	!! 4. liquid water flow
	1469	ztotwcap = SUM(zwholdmax(:))
	1470	! Rain entering snow (m):
	1471	!ORCHIDEE has assumed that all precipitation entering snow has
[5454]	1472	!left the snowpack and finally become runoff in hydrol_soil. Here we just put zrainfall as 0.0
[2650]	1473	!!zrainfall = precip_rain(ji)/ph2o ! rainfall (m)
	1474	zrainfall = 0.0
[2223]	1475
[2650]	1476	zflowliqt(0) = MIN(zrainfall,ztotwcap)
	1477	! Rain assumed to directly pass through the pack to runoff (m):
	1478	zpcpxs = zrainfall - zflowliqt(0)
	1479	!
	1480	DO jj=1,nsnow
	1481	zflowliqt(jj) = flowliq(jj)
	1482	ENDDO
[2223]	1483
[2650]	1484	! translated into a density increase:
	1485	flowliq(:) = 0.0 ! clear this array for work
	1486	zsnowliq(:) = snowliq(ji,:) ! reset liquid water content
	1487	!
[2223]	1488
[2650]	1489	DO jj=1,nsnow
	1490	snowliq(ji,jj) = snowliq(ji,jj) + zflowliqt(jj-1)
	1491	flowliq(jj) = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj)))
	1492	snowliq(ji,jj) = snowliq(ji,jj) - flowliq(jj)
[2223]	1493
[2650]	1494	!Modified by Tao Wang based on previous ISBA-ES scheme
	1495	snowrho(ji,jj) = snowrho(ji,jj) + (snowliq(ji,jj) - zsnowliq(jj))* &
	1496	& ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj))
[2223]	1497
[2650]	1498	zflowliqt(jj) = zflowliqt(jj) + flowliq(jj)
	1499	ENDDO
[2223]	1500
[2650]	1501	snowmelt(ji) = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) * ph2o
[2223]	1502
[2650]	1503	! excess heat from melting, using it to warm underlying ground to conserve energy
	1504	meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2)
[2223]	1505
[2650]	1506	! energy flux into the soil
	1507	grndflux(ji) = grndflux(ji) + meltxs(ji)
[2223]	1508
[3772]	1509	ELSE
	1510	snowdz(ji,:)=0.
	1511	snowliq(ji,:)=0.
	1512	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
[2650]	1513	ENDIF
[2223]	1514
[2650]	1515	ENDDO
[2223]	1516
[2650]	1517	END SUBROUTINE explicitsnow_melt_refrz
[2223]	1518
[2650]	1519	!================================================================================================================================
	1520	!! SUBROUTINE : explicitsnow_levels
	1521	!!
	1522	!>\BRIEF Computes snow discretization based on given total snow depth
	1523	!!
	1524	!! DESCRIPTION :
	1525	!! RECENT CHANGE(S) : None
	1526	!!
	1527	!! MAIN OUTPUT VARIABLE(S): None
	1528	!!
	1529	!! REFERENCE(S) :
	1530	!!
	1531	!! FLOWCHART : None
	1532	!! \n
	1533	!_
	1534	!================================================================================================================================
[2223]	1535	SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz)
[2650]	1536
[2223]	1537	!! 0.1 Input variables
	1538	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
	1539	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_thick !! Total snow depth
	1540
	1541	!! 0.2 Output variables
	1542
	1543	!! 0.3 Modified variables
	1544
[2650]	1545	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
[2223]	1546
	1547	!! 0.4 Local variables
	1548
	1549	INTEGER(i_std) :: il,it,ji
	1550	INTEGER(i_std) :: i,j, ix
	1551	REAL(r_std), DIMENSION(kjpindex) :: xi, xf
	1552	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF1 = (/0.25, 0.50, 0.25/) !! Snow grid parameters
	1553	REAL(r_std), PARAMETER, DIMENSION(2) :: ZSGCOEF2 = (/0.05, 0.34/) !! Snow grid parameters
	1554	REAL(r_std), PARAMETER :: ZSNOWTRANS = 0.20 !! Minimum total snow depth at which surface layer thickness is constant (m)
	1555	REAL(r_std), PARAMETER :: XSNOWCRITD = 0.03 !! (m)
	1556
	1557	DO ji=1,kjpindex
	1558	IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN
	1559
[2650]	1560	snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow)
	1561	snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow)
	1562	snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3)
	1563
	1564	ENDIF
	1565	ENDDO
[2223]	1566
	1567	WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. &
[2650]	1568	snow_thick(:) .GT. (XSNOWCRITD+0.01) )
[2223]	1569	!
[2650]	1570	snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1)
	1571	snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2)
	1572	snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3)
[2223]	1573	!
	1574	END WHERE
	1575
[2650]	1576	DO ji = 1,kjpindex
	1577	IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN
[2223]	1578	snowdz(ji,1) = ZSGCOEF2(1)
	1579	snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1)
[2650]	1580	!When using simple finite differences, limit the thickness
	1581	!factor between the top and 2nd layers to at most 10
[2223]	1582	snowdz(ji,2) = MIN(10*ZSGCOEF2(1), snowdz(ji,2) )
	1583	snowdz(ji,3) = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1)
	1584
[2650]	1585	ENDIF
	1586	ENDDO
	1587
[2223]	1588
	1589	END SUBROUTINE explicitsnow_levels
	1590
[3269]	1591	!!
	1592	!================================================================================================================================
	1593	!! SUBROUTINE : explicitsnow_profile
	1594	!!
	1595	!>\BRIEF
	1596	!!
	1597	!! DESCRIPTION : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile.
	1598	!!
	1599	!! RECENT CHANGE(S) : None
	1600	!!
	1601	!! MAIN OUTPUT VARIABLE(S):
	1602	!!
	1603	!! REFERENCE(S) :
	1604	!!
	1605	!! FLOWCHART : None
	1606	!! \n
	1607	!_
	1608	!================================================================================================================================
	1609	SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
[2223]	1610
[3269]	1611	!! 0. Variables and parameter declaration
	1612
	1613	!! 0.1 Input variables
	1614
	1615	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	1616	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! skin temperature
	1617	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
	1618	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
	1619	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	1620	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness
	1621	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
	1622
	1623	!! 0.3 Modified variable
	1624
	1625	!! 0.2 Output variables
	1626
	1627	!! 0.3 Modified variables
	1628	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout) :: snowtemp
	1629
	1630	!! 0.4 Local variables
	1631
	1632	INTEGER(i_std) :: ji, jg
	1633	!_
	1634	!================================================================================================================================
	1635	!! 1. Computes the snow temperatures ptn.
	1636	DO ji = 1,kjpindex
	1637	IF (SUM(snowdz(ji,:)) .GT. 0) THEN
	1638	snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
	1639	(lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un)
	1640	temp_sol_add(ji) = zero
	1641	DO jg = 1,nsnow-1
	1642	snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg)
	1643	ENDDO
	1644	ENDIF
	1645	ENDDO
	1646	IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done '
	1647
	1648	END SUBROUTINE explicitsnow_profile
	1649
[2223]	1650	END MODULE explicitsnow

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

Download in other formats: