Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_sechiba/qsat_moisture.f90 @ 7346

Last change on this file since 7346 was 5530, checked in by anne-sofie.lanso, 6 years ago
UPDATED ipslerr in qsatcalc to crash instead of giving warnings, when there is problems for the temperature. In particular, when you run with the multi-layer energy budget these problems will arise, and NaN can propergate to other variables in the code.
Property svn:keywords set to `Date Revision`
File size: 41.3 KB

Line
1	! =================================================================================================================================
2	! MODULE : qsat_moisture
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2011)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF "qsat_moisture" module contains public tools functions like qsat, dev_qsat.
10	!!
11	!!\n DESCRIPTION: This module is the result of the splitting of constantes_veg.\n
12	!! As the subroutines qsatcalc, dev_qsatcalc are used only by enerbil and diffuco, they are part of SECHIBA
13	!! component.
14	!!
15	!! REFERENCE(S) :
16	!!
17	!! SVN :
18	!! $HeadURL: $
19	!! $Date$
20	!! $Revision$
21	!! \n
22	!_ ================================================================================================================================
23
24	MODULE qsat_moisture
25
26	USE defprec
27	USE constantes
28	USE IOIPSL
29	USE constantes_soil
30
31	IMPLICIT NONE
32	PRIVATE
33	PUBLIC qsatcalc, dev_qsatcalc
34	PUBLIC snow3lhold_2d, snow3lhold_1d, snow3lhold_0d
35	PUBLIC snow3lheat_2d, snow3lheat_1d
36	PUBLIC snow3lscap_2d, snow3lscap_1d
37	PUBLIC snow3ltemp_2d, snow3ltemp_1d
38	PUBLIC snow3lgrain_2d, snow3lgrain_1d, snow3lgrain_0d
39	PUBLIC snow3lliq_2d, snow3lliq_1d
40
41
42	LOGICAL,SAVE :: l_qsat_first=.TRUE. !! First call to qsat subroutines and functions (true/false)
43	!$OMP THREADPRIVATE(l_qsat_first)
44
45
46	INTEGER(i_std),PARAMETER :: max_temp=370 !! Maximum temperature for saturated humidity (K) and also used as
47	!! the size of local array to keep saturated humidity (unitless)
48
49	INTEGER(i_std),PARAMETER :: min_temp=100 !! Minimum temperature for saturated humidity (K)
50
51	REAL(r_std),DIMENSION(max_temp),SAVE :: qsfrict !! Array to keep water vapor pressure at saturation for each temperature level
52	!! (hPa)
53	!$OMP THREADPRIVATE(qsfrict)
54
55	CONTAINS
56
57	!! ================================================================================================================================
58	!! SUBROUTINE : qsatcalc
59	!!
60	!>\BRIEF This routine calculates the saturated humidity using the pressure
61	!! and the temperature for all pixels.
62	!!
63	!! DESCRIPTION : This routine interpolates qsat between temperatures by the following formula :
64	!! \latexonly
65	!! \input{qsatcalc.tex}
66	!! \endlatexonly
67	!! \n
68	!!
69	!! RECENT CHANGE(S): None
70	!!
71	!! MAIN OUTPUT VARIABLE(S) : qsat_out
72	!!
73	!! REFERENCE(S) : None
74	!!
75	!! FLOWCHART : None
76	!! \n
77	!_ ================================================================================================================================
78
79	SUBROUTINE qsatcalc (kjpindex,temp_in,pres_in,qsat_out)
80
81	IMPLICIT NONE
82
83	!! 0. Variables and parameters declaration
84
85	!! 0.1 Input variables
86
87	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (unitless)
88	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_in !! Temperature in degre Kelvin (K)
89	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: pres_in !! Pressure (hPa)
90
91	!! 0.2 Output variables
92
93	REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: qsat_out !! Saturated humidity at the surface (kg of water/kg of air)
94
95	!! 0.4 Local variables
96
97	INTEGER(i_std), DIMENSION(kjpindex) :: jt !! Temporary array stocking the truncated temperatures in Kelvin
98	!!(converted into integers)
99	INTEGER(i_std) :: ji !! indices (unitless)
100	REAL(r_std),DIMENSION(kjpindex) :: zz_a, zz_b, zz_f !! Temporary variables
101	INTEGER(i_std) :: nbad !! Number of points where the temperature is too high or too low
102	INTEGER(i_std),DIMENSION(1) :: lo !! Temporary vector to mark the position of the highest temperature
103	!! or the lowest temperature over all the pixels in jt (unitless)
104	!_ ================================================================================================================================
105
106	!-
107	!! 1.Initialize qsfrict array if needed
108	!-
109	IF (l_qsat_first) THEN
110	!-
111	CALL qsfrict_init
112	l_qsat_first = .FALSE.
113	!-
114	ENDIF !(l_qsat_first)
115
116	!-
117	!! 2. Computes qsat interpolation into two successive temperature
118	!-
119	jt = INT(temp_in(:))
120
121	!! 2.1 Diagnostic pixels where the temperature is too high
122	nbad = COUNT(jt(:) >= max_temp-1)
123
124	!! When running with the 10 layer multi-layer energy budget, problems
125	! aries here in qsatcalc. This can result in NaN that propergates through
126	! to other variable. Therefore, instead of just giving a warning, when the
127	! temperature is incorrect in qsatcal, we will now let the model crash. The
128	! problems we have with mleb and qsatcalc is most likely caused by the
129	! fact that we are not running with tall trees (it should run fine with
130	! tall trees). Thus, once we activate the mleb 10 layer agian, we need to
131	! look into this matter (ASL).
132	IF (nbad > 0) THEN
133	WRITE(numout,*) ' qsatcalc: temperature too high at ', &
134	& nbad, ' points.'
135	!-
136	IF (.NOT.diag_qsat) THEN
137	CALL ipslerr_p(3,'qsatcalc','diffuco', '', &
138	& 'temperature incorect.') ! Warning message
139	ELSE
140	lo = MAXLOC(temp_in(:))
141	WRITE(numout,*) &
142	& 'Maximum temperature ( ',MAXVAL(temp_in),') found at ',lo(1)
143	WHERE (jt(:) >= max_temp-1) jt(:) = max_temp-1
144	ENDIF !(.NOT.diag_qsat)
145	!-
146	ENDIF ! (nbad > 0)
147
148
149	!! 2.2 Diagnostic pixels where the temperature is too low
150	nbad = COUNT(jt(:) <= min_temp)
151
152	IF (nbad > 0) THEN
153	WRITE(numout,*) ' qsatcalc: temperature too low at ', &
154	& nbad, ' points.'
155	!-
156	IF (.NOT.diag_qsat) THEN
157	CALL ipslerr_p(3,'qsatcalc','diffuco', '', &
158	& 'temperature incorect.') ! Warning message
159	ELSE
160	lo = MINLOC(temp_in(:))
161	WRITE(numout,*) &
162	& 'Minimum temperature ( ',MINVAL(temp_in),') found at ',lo(1)
163	WHERE (jt(:) <= min_temp) jt(:) = min_temp
164	ENDIF !(.NOT.diag_qsat)
165	!-
166	ENDIF! (nbad > 0)
167
168	!! 2.3 Temporary variables needed for interpolation
169	DO ji = 1, kjpindex ! Loop over # pixels
170
171	zz_f(ji) = temp_in(ji)-FLOAT(jt(ji))
172	zz_a(ji) = qsfrict(jt(ji))
173	zz_b(ji) = qsfrict(jt(ji)+1)
174
175	ENDDO ! Loop over # pixels
176
177	!-
178	!! 3. Interpolation between these two values
179	!-
180	DO ji = 1, kjpindex ! Loop over # pixels
181
182	qsat_out(ji) = ((zz_b(ji)-zz_a(ji))*zz_f(ji)+zz_a(ji))/pres_in(ji)
183
184	ENDDO ! Loop over # pixels
185
186
187	END SUBROUTINE qsatcalc
188
189	!! ================================================================================================================================
190	!! FUNCTION : [DISPENSABLE] qsat_old
191	!!
192	!>\BRIEF This function computes deviation the saturated humidity with the pressure
193	!! and the temperature for a scalar.
194	!!
195	!! DESCRIPTION : This routine is obsolete : replaced by the subroutine qsatcalc. \n
196	!! qsat is interpolated by : \n
197	!! \latexonly
198	!! \input{qsat.tex}
199	!! \endlatexonly
200	!!
201	!! RECENT CHANGE(S): None\n
202	!!
203	!! RETURN VALUE : qsat_result
204	!!
205	!! REFERENCE(S) : None
206	!!
207	!! FLOWCHART : None
208	!! \n
209	!_ ================================================================================================================================
210
211	FUNCTION qsat_old (temp_in,pres_in) RESULT (qsat_result)
212
213	IMPLICIT NONE
214
215	!! 0. Variables and parameters declaration
216
217	!! 0.1 Input variables
218
219	REAL(r_std),INTENT(in) :: temp_in !! Temperature (K)
220	REAL(r_std),INTENT(in) :: pres_in !! Pressure (hPa)
221
222	!! 0.2 Result
223
224	REAL(r_std) :: qsat_result !! Saturated humidity calculated at the surface (kg/kg)
225
226	!! 0.4 Local variables
227
228	INTEGER(i_std) :: jt !! Temporary scalar stocking the truncated temperature in Kelvin
229	!! (converted into integer)
230	REAL(r_std) :: zz_a,zz_b,zz_f !! Temporary scalar variables
231
232	!_ ================================================================================================================================
233
234	!-
235	!! 1.Initialize qsfrict array if needed
236	!-
237	IF (l_qsat_first) THEN
238	!-
239	CALL qsfrict_init
240	l_qsat_first = .FALSE.
241	!-
242	ENDIF
243
244	!-
245	!! 2. Computes qsat interpolation into two successive temperatures
246	!-
247	jt = INT(temp_in)
248
249	!! 2.1 Is the temperature too high ?
250	IF (jt >= max_temp-1) THEN
251	WRITE(numout,*) &
252	& ' We stop. temperature too BIG : ',temp_in, &
253	& ' approximation for : ',jt
254	!-
255	IF (.NOT.diag_qsat) THEN
256	CALL ipslerr_p(2,'qsat','', '',&
257	& 'temperature incorect.') ! Warning message
258	ELSE
259	qsat_result = 999999.
260	RETURN
261	ENDIF !(.NOT.diag_qsat)
262	!-
263	ENDIF !(jt >= max_temp-1)
264
265	!! 2.2 Is the temperature too low ?
266	IF (jt <= min_temp ) THEN
267	WRITE(numout,*) &
268	& ' We stop. temperature too SMALL : ',temp_in, &
269	& ' approximation for : ',jt
270	!-
271	IF (.NOT.diag_qsat) THEN
272	CALL ipslerr_p(2,'qsat','', '',&
273	& 'temperature incorect.')
274	ELSE
275	qsat_result = -999999.
276	RETURN
277	ENDIF!(.NOT.diag_qsat)
278	!-
279	ENDIF !(jt <= min_temp )
280
281	!! 2.3 Temporary variables needed for interpolation
282	zz_f = temp_in-FLOAT(jt)
283	zz_a = qsfrict(jt)
284	zz_b = qsfrict(jt+1)
285
286	!! 3. Interpolates between these two values
287
288	qsat_result = ((zz_b-zz_a)*zz_f+zz_a)/pres_in
289
290
291	END FUNCTION qsat_old
292
293
294	!! ================================================================================================================================
295	!! SUBROUTINE : dev_qsatcalc
296	!!
297	!>\BRIEF This routine calculates the deviation of the saturated humidity qsat.
298	!!
299	!! DESCRIPTION : The deviation of qsat is calculated by :
300	!! \latexonly
301	!! \input{dev_qsatcalc.tex}
302	!! \endlatexonly
303	!!
304	!! RECENT CHANGE(S): None
305	!!
306	!! MAIN OUTPUT VARIABLE(S) : dev_qsat_out
307	!!
308	!! REFERENCE(S) : None
309	!!
310	!! FLOWCHART : None
311	!!
312	!! FLOWCHART :
313	!! \latexonly
314	!! \includegraphics[scale = 1]{pheno_moigdd.png}
315	!! \endlatexonly
316	!! \n
317	!_ ================================================================================================================================
318
319	SUBROUTINE dev_qsatcalc (kjpindex,temp_in,pres_in,dev_qsat_out)
320
321	IMPLICIT NONE
322
323	!! 0. Variables and parameters declaration
324
325	!! 0.1 Input variables
326
327	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (unitless)
328	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_in !! Temperature (K)
329	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: pres_in !! Pressure (hPa)
330
331
332	!! 0.2 Output variables
333
334	REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: dev_qsat_out !! Result (??units??)
335
336
337	!! 0.4 Local variables
338
339	INTEGER(i_std),DIMENSION(kjpindex) :: jt !! Temporary array stocking the truncated temperatures
340	!! in Kelvin (converted into integers)
341	INTEGER(i_std) :: ji !! Indice (unitless)
342	REAL(r_std),DIMENSION(kjpindex) :: zz_a, zz_b, zz_c, zz_f !! Temporary vector variables
343	INTEGER(i_std) :: nbad !! Number of points where the temperature is too high or too low
344
345	!_ ================================================================================================================================
346
347	!-
348	!! 1.Initialize qsfrict array if needed
349	!-
350	IF (l_qsat_first) THEN
351	!-
352	CALL qsfrict_init
353	l_qsat_first = .FALSE.
354	!-
355	ENDIF
356
357	!-
358	!! 2. Compute qsat interpolation into two successive temperature
359	!-
360	jt = INT(temp_in(:)+undemi)
361
362	!! 2.1 Pixels where the temperature is too high
363	nbad = COUNT( jt(:) >= max_temp-1 )
364
365	IF (nbad > 0) THEN
366	WRITE(numout,*) &
367	& ' dev_qsatcalc: temperature too high at ',nbad,' points.'
368	!-
369	IF (.NOT.diag_qsat) THEN
370	CALL ipslerr_p(3,'dev_qsatcalc','', '', &
371	& 'temperature incorect.') ! Fatal error
372	ELSE
373	WHERE (jt(:) >= max_temp-1) jt(:) = max_temp-1
374	ENDIF !(.NOT.diag_qsat)
375	!-
376	ENDIF !(nbad > 0)
377
378	!! 2.2 Pixels where the temperature is too low
379	nbad = COUNT( jt(:) <= min_temp )
380
381	IF (nbad > 0) THEN
382	WRITE(numout,*) &
383	& ' dev_qsatcalc: temperature too low at ',nbad,' points.'
384	!-
385	IF (.NOT.diag_qsat) THEN
386	CALL ipslerr_p(3,'dev_qsatcalc', '', '',&
387	& 'temperature incorect.') ! Fatal error
388	ELSE
389	WHERE (jt(:) <= min_temp) jt(:) = min_temp
390	ENDIF !(.NOT.diag_qsat)
391	!-
392	ENDIF !(nbad > 0)
393
394	!! 2.3 Temporary variables needed for interpolation
395	DO ji=1,kjpindex ! Loop over # pixels
396
397	zz_f(ji) = temp_in(ji)+undemi-FLOAT(jt(ji))
398	zz_a(ji) = qsfrict(jt(ji)-1)
399	zz_b(ji) = qsfrict(jt(ji))
400	zz_c(ji) = qsfrict(jt(ji)+1)
401
402	ENDDO ! Loop over # pixels
403
404	!-
405	!! 3. Interpolates between these two values
406	!-
407	DO ji = 1, kjpindex ! Loop over # pixels
408
409	dev_qsat_out(ji) = &
410	& ((zz_c(ji)-deuxzz_b(ji)+zz_a(ji))(zz_f(ji)-un) + &
411	& zz_c(ji)-zz_b(ji))/pres_in(ji)
412
413	ENDDO ! Loop over # pixels
414
415
416	END SUBROUTINE dev_qsatcalc
417
418	!! ================================================================================================================================
419	!! FUNCTION : [DISPENSABLE] dev_qsat_old
420	!!
421	!>\BRIEF This function computes deviation of qsat.
422	!!
423	!! DESCRIPTION : The deviation of qsat is calculated by :
424	!! \latexonly
425	!! \input{dev_qsat.tex}
426	!! \endlatexonly
427	!!
428	!! RECENT CHANGE(S): None
429	!!
430	!! RETURN VALUE : dev_qsat_result
431	!!
432	!! REFERENCE(S) : None
433	!!
434	!! FLOWCHART : None
435	!! \n
436	!_ ================================================================================================================================
437
438	FUNCTION dev_qsat_old (temp_in,pres_in) RESULT (dev_qsat_result)
439
440	IMPLICIT NONE
441
442	!! 0. Variables and parameters declaration
443
444	!! 0.1 Input variables
445
446	REAL(r_std),INTENT(in) :: pres_in !! Pressure (hPa)
447	REAL(r_std),INTENT(in) :: temp_in !! Temperture (K)
448
449	!! 0.2 Result
450
451	REAL(r_std) :: dev_qsat_result !! (??units??) !!
452
453	!! 0.4 Local variables
454
455	INTEGER(i_std) :: jt !! Index (unitless)
456	REAL(r_std) :: zz_a, zz_b, zz_c, zz_f !! Temporary scalars
457
458	!_ ================================================================================================================================
459
460	!-
461	!! 1.Initialize qsfrict array if needed
462	!-
463	IF (l_qsat_first) THEN
464	!-
465	CALL qsfrict_init
466	l_qsat_first = .FALSE.
467	!-
468	ENDIF
469
470	!-
471	!! 2. computes qsat deviation interpolation
472	!! into two successive temperature
473	!-
474	jt = INT(temp_in+undemi)
475
476	!! 2.1 Is the temperature too high ?
477	IF (jt >= max_temp-1) THEN
478	!-
479	WRITE(numout,*) &
480	& ' We stop. temperature too HIGH : ',temp_in, &
481	& ' approximation for : ',jt
482	IF (.NOT.diag_qsat) THEN
483	CALL ipslerr_p(3,'dev_qsat','', '',&
484	& 'temperature incorect.') ! Fatal error
485	ELSE
486	dev_qsat_result = 999999.
487	RETURN
488	ENDIF !(.NOT.diag_qsat)
489	!-
490	ENDIF !(jt >= max_temp-1)
491	!-
492	!! 2.2 Is the temperature too low ?
493	IF (jt <= min_temp ) THEN
494	WRITE(numout,*) &
495	& ' We stop. temperature too LOW : ',temp_in, &
496	& ' approximation for : ',jt
497	!-
498	IF (.NOT.diag_qsat) THEN
499	CALL ipslerr_p(3,'dev_qsat','', '',&
500	& 'temperature incorect.')
501	ELSE
502	dev_qsat_result = -999999.
503	RETURN
504	ENDIF !(.NOT.diag_qsat)
505	!-
506	ENDIF !(jt <= min_temp )
507
508	!! 2.3 Temporary variables for interpolation
509	zz_f = temp_in+undemi-FLOAT(jt)
510	zz_a = qsfrict(jt-1)
511	zz_b = qsfrict(jt)
512	zz_c = qsfrict(jt+1)
513
514	!-
515	!! 3. Interpolate
516	!-
517	dev_qsat_result=((zz_c-deuxzz_b+zz_a)(zz_f-un)+zz_c-zz_b)/pres_in
518
519	END FUNCTION dev_qsat_old
520
521
522	!! ================================================================================================================================
523	!! SUBROUTINE : qsfrict_init
524	!!
525	!>\BRIEF The qsfrict_init routine initialises qsfrict array to store
526	!! precalculated values for qsat by using Goff-Gratch equations.
527	!!
528	!! DESCRIPTION : This routine calculates the specific humidity qsat as a function of temperature in
529	!! Kelvin by using the modified Goff-Gratch equations(1946): \n
530	!! \latexonly
531	!! \input{goff_gratch.tex}
532	!! \endlatexonly
533	!! qsfrict is initialized by the following formulas : \n
534	!! \latexonly
535	!! \input{qsfrict_init.tex}
536	!! \endlatexonly
537	!! These values are used by the subroutines qsatcalc, dev_qsat. \n
538	!!
539	!! RECENT CHANGE(S): None
540	!!
541	!! MAIN OUTPUT VARIABLE(S): ::qsfrict
542	!!
543	!! REFERENCE(S) :
544	!! - Algorithme d'un ensemble de paramÃ©trisation physique (1998),
545	!! Note de Laurent Li dÃ©crivant les paramÃ©trisations physiques incluses dans le modÃšle (pdf),
546	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques
547	!! - Goff, J. A., and S. Gratch (1946) Low-pressure properties of water from â160 to 212 Â°F, in Transactions of the
548	!! American Society of Heating and Ventilating Engineers, pp 95â122, presented at the 52nd annual meeting of the
549	!! American Society of Heating and Ventilating Engineers, New York, 1946.
550	!!
551	!! FLOWCHART : None
552	!! \n
553	!_ ================================================================================================================================
554
555	SUBROUTINE qsfrict_init
556
557	IMPLICIT NONE
558
559	!! 0. Variables and parameters declaration
560
561	!! 0.4 Local variables
562
563	INTEGER(i_std) :: ji !! Indice(unitless)
564	REAL(r_std) :: zrapp,zcorr,ztemperature,zqsat !! Temporary vector variables
565
566	!_ ================================================================================================================================
567
568	!! 1. Initialisation
569	zrapp = msmlr_h2o/msmlr_air
570	zcorr = 0.00320991_r_std
571
572	!! 2. Computes saturated humidity one time and store in qsfrict local array
573	DO ji=100,max_temp ! Loop over size(qsfrict) : each position of qsfrict matches a temperature
574
575	ztemperature = FLOAT(ji)
576	!-
577	IF (ztemperature < 273._r_std) THEN
578	zqsat = zrapp10.0_r_std(2.07023_r_std-zcorrztemperature &
579	& -2484.896/ztemperature+3.56654*LOG10(ztemperature)) ! Equilibrium water vapor - solid
580	ELSE
581	zqsat = zrapp10.0*(23.8319-2948.964/ztemperature &
582	& -5.028*LOG10(ztemperature) &
583	& -29810.16EXP(-0.0699382ztemperature) &
584	& +25.21935*EXP(-2999.924/ztemperature)) ! Equilibrium water vapor - liquid
585	ENDIF !(ztemperature < 273._r_std)
586	!-
587	qsfrict (ji) = zqsat
588
589	ENDDO ! Loop over size(qsfrict)
590
591	!! 3. Set to zero the non-computed values
592	qsfrict(1:100) = zero
593	!-
594	IF (printlev>=3) WRITE (numout,*) ' qsfrict_init done'
595
596
597	END SUBROUTINE qsfrict_init
598
599	!!
600	!================================================================================================================================
601	!! FUNCTION : snow3lhold_2d
602	!!
603	!>\BRIEF Calculate the maximum liquid water holding capacity of
604	!! snow layer(s)
605	!! DESCRIPTION :
606	!!
607	!! RECENT CHANGE(S): None
608	!!
609	!! MAIN OUTPUT VARIABLE(S): :: PWHOLDMAX
610	!!
611	!! REFERENCE(S) :
612	!!
613	!! FLOWCHART : None
614	!! \n
615	!_
616	!================================================================================================================================
617
618	FUNCTION snow3lhold_2d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
619
620	!! 0.1 Input variables
621	REAL(r_std), DIMENSION(:,:), INTENT(IN) :: PSNOWDZ !! Snow depth
622	REAL(r_std), DIMENSION(:,:), INTENT(IN) :: PSNOWRHO !! Snow density
623
624	!! 0.2 Output variables
625	REAL(r_std), DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PWHOLDMAX !! Maximum Water holding capacity
626
627	!! 0.3 Modified variables
628
629	!! 0.4 Local variables
630
631	REAL(r_std), DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZHOLDMAXR, ZSNOWRHO
632
633
634	! Evaluate capacity using upper density limit:
635	ZSNOWRHO(:,:) = MIN(xrhosmax, PSNOWRHO(:,:))
636
637	! Maximum ratio of liquid to SWE:
638	ZHOLDMAXR(:,:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)* &
639	MAX(0.,xsnowrhohold-ZSNOWRHO(:,:))/xsnowrhohold
640
641	! Maximum liquid water holding capacity of the snow (m):
642	PWHOLDMAX(:,:) = ZHOLDMAXR(:,:)PSNOWDZ(:,:)ZSNOWRHO(:,:)/ph2o
643	WHERE(ZSNOWRHO(:,:) .GE. xrhosmax) PWHOLDMAX(:,:) = 0.0
644
645	END FUNCTION snow3lhold_2d
646
647
648	!!
649	!================================================================================================================================
650	!! FUNCTION : snow3lhold_1d
651	!!
652	!>\BRIEF Calculate the maximum liquid water holding capacity of
653	!! snow layer(s)
654	!! DESCRIPTION :
655	!!
656	!! RECENT CHANGE(S): None
657	!!
658	!! MAIN OUTPUT VARIABLE(S): ::
659	!!
660	!! REFERENCE(S) :
661	!!
662	!! FLOWCHART : None
663	!! \n
664	!_
665	!================================================================================================================================
666
667	FUNCTION snow3lhold_1d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
668
669	!! 0.1 Input variables
670	REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ !! Snow depth
671	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! Snow density
672
673	!! 0.2 Output variables
674	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PWHOLDMAX !! Maximum Water holding capacity
675
676	!! 0.3 Modified variables
677
678	!! 0.4 Local variables
679	REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZHOLDMAXR, ZSNOWRHO
680
681
682	! Evaluate capacity using upper density limit:
683	ZSNOWRHO(:) = MIN(xrhosmax, PSNOWRHO(:))
684
685	! Maximum ratio of liquid to SWE:
686	ZHOLDMAXR(:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)* &
687	MAX(0.,xsnowrhohold-ZSNOWRHO(:))/xsnowrhohold
688
689	! Maximum liquid water holding capacity of the snow (m):
690	PWHOLDMAX(:) = ZHOLDMAXR(:)PSNOWDZ(:)ZSNOWRHO(:)/ph2o
691
692	WHERE(ZSNOWRHO(:) .GE. xrhosmax) PWHOLDMAX(:)=0.0
693
694	END FUNCTION snow3lhold_1d
695
696	!!
697	!================================================================================================================================
698	!! FUNCTION : snow3lhold_0d
699	!!
700	!>\BRIEF Calculate the maximum liquid water holding capacity of
701	!! snow layer(s)
702	!! DESCRIPTION :
703	!!
704	!! RECENT CHANGE(S): None
705	!!
706	!! MAIN OUTPUT VARIABLE(S): ::
707	!!
708	!! REFERENCE(S) :
709	!!
710	!! FLOWCHART : None
711	!! \n
712	!_
713	!================================================================================================================================
714
715	FUNCTION snow3lhold_0d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
716
717	!! 0.1 Input variables
718	REAL(r_std), INTENT(IN) :: PSNOWRHO !!
719	!! Snow density
720	REAL(r_std), INTENT(IN) :: PSNOWDZ !!
721	!! Snow depth
722
723	!! 0.2 Output variables
724	REAL(r_std) :: PWHOLDMAX !!
725	!! Maximum water holding capacity
726
727	!! 0.3 Modified variables
728
729	!! 0.4 Local variables
730	REAL(r_std) :: ZHOLDMAXR, ZSNOWRHO
731
732
733	! Evaluate capacity using upper density limit:
734	ZSNOWRHO = MIN(xrhosmax, PSNOWRHO)
735
736	! Maximum ratio of liquid to SWE:
737	ZHOLDMAXR = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*&
738	& MAX(0.,xsnowrhohold-ZSNOWRHO)/xsnowrhohold
739
740	! Maximum liquid water holding capacity of the snow (m):
741	PWHOLDMAX = ZHOLDMAXRPSNOWDZZSNOWRHO/ph2o
742
743	IF (ZSNOWRHO .GE. xrhosmax) PWHOLDMAX = 0.0
744
745	END FUNCTION snow3lhold_0d
746
747	!!
748	!================================================================================================================================
749	!! FUNCTION : snow3lheat_2d
750	!!
751	!>\BRIEF Compute snow heat content (J m-2) from snow mass and liquid
752	!! water content and temperature.
753	!! snow layer(s)
754	!! DESCRIPTION :
755	!!
756	!! RECENT CHANGE(S): None
757	!!
758	!! MAIN OUTPUT VARIABLE(S): ::
759	!!
760	!! REFERENCE(S) :
761	!!
762	!! FLOWCHART : None
763	!! \n
764	!_
765	!================================================================================================================================
766
767	FUNCTION snow3lheat_2d(PSNOWLIQ,PSNOWRHO,PSNOWDZ,PSNOWTEMP) RESULT(PSNOWHEAT)
768
769	!! 0.1 Input variables
770	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO !! layer density (kg m-3)
771	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ !! layer thickness (m)
772	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWLIQ !! liquid water content (m)
773	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP !! layer temperature (K)
774
775	!! 0.2 Output variables
776	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PSNOWHEAT !! heat content (enthalpy) (J m-2)
777
778	!! 0.3 Modified variables
779
780	!! 0.4 Local variables
781	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSCAP !! snow heat capacity (J K-1 m-3)
782
783	ZSCAP(:,:) = snow3lscap_2d(PSNOWRHO)
784
785	! snow heat content (heat required to melt the snowpack) or enthalpy (J m-2)
786	PSNOWHEAT(:,:) = PSNOWDZ(:,:)( ZSCAP(:,:)(PSNOWTEMP(:,:)-tp_00) &
787	- chalfu0PSNOWRHO(:,:) ) + chalfu0ph2o*PSNOWLIQ(:,:)
788
789	END FUNCTION snow3lheat_2d
790
791	!!
792	!================================================================================================================================
793	!! FUNCTION : snow3lheat_1d
794	!!
795	!>\BRIEF Compute snow heat content (J m-2) from snow mass and liquid
796	!! water content and temperature.
797	!! snow layer(s)
798	!! DESCRIPTION :
799	!!
800	!! RECENT CHANGE(S): None
801	!!
802	!! MAIN OUTPUT VARIABLE(S): ::
803	!!
804	!! REFERENCE(S) :
805	!!
806	!! FLOWCHART : None
807	!! \n
808	!_
809	!================================================================================================================================
810
811	FUNCTION snow3lheat_1d(PSNOWLIQ,PSNOWRHO,PSNOWDZ,PSNOWTEMP) RESULT(PSNOWHEAT)
812
813	!! 0.1 Input variables
814	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! layer density (kg m-3)
815	REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ !! layer thickness (m)
816	REAL, DIMENSION(:), INTENT(IN) :: PSNOWLIQ !! liquid water content (m)
817	REAL, DIMENSION(:), INTENT(IN) :: PSNOWTEMP !! layer temperature (K)
818
819	!! 0.2 Output variables
820	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PSNOWHEAT !! heat content (enthalpy) (J m-2)
821
822	!! 0.3 Modified variables
823
824	!! 0.4 Local variables
825	REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZSCAP !! snow heat capacity (J K-1 m-3)
826
827
828	ZSCAP(:) = snow3lscap_1d(PSNOWRHO)
829
830	! snow heat content (heat required to melt the snowpack) or enthalpy (J m-2)
831	PSNOWHEAT(:) = PSNOWDZ(:)( ZSCAP(:)(PSNOWTEMP(:)-tp_00) &
832	-chalfu0PSNOWRHO(:) ) + chalfu0ph2o*PSNOWLIQ(:)
833
834	END FUNCTION snow3lheat_1d
835
836	!!
837	!================================================================================================================================
838	!! FUNCTION : snow3lscap_2d
839	!!
840	!>\BRIEF Calculate the heat capacity of a snow layer.
841	!!
842	!! DESCRIPTION :
843	!!
844	!! RECENT CHANGE(S): None
845	!!
846	!! MAIN OUTPUT VARIABLE(S): ::
847	!!
848	!! REFERENCE(S) : The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
849	!!
850	!! FLOWCHART : None
851	!! \n
852	!_
853	!================================================================================================================================
854	FUNCTION snow3lscap_2d(PSNOWRHO) RESULT(PSCAP)
855
856	!! 0.1 Input variables
857	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO !! Snow density
858
859	!! 0.2 Output variables
860	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PSCAP !! Heat capacity (J K-1 m-3)
861
862	PSCAP(:,:) = PSNOWRHO(:,:)*xci
863
864	END FUNCTION snow3lscap_2d
865
866	!!
867	!================================================================================================================================
868	!! FUNCTION : snow3lscap_1d
869	!!
870	!>\BRIEF Calculate the heat capacity of a snow layer.
871	!!
872	!! DESCRIPTION :
873	!!
874	!! RECENT CHANGE(S): None
875	!!
876	!! MAIN OUTPUT VARIABLE(S): ::
877	!!
878	!! REFERENCE(S) : The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
879	!!
880	!! FLOWCHART : None
881	!! \n
882	!_
883	!================================================================================================================================
884	FUNCTION snow3lscap_1d(PSNOWRHO) RESULT(PSCAP)
885
886	!! 0.1 Input variables
887	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! Snow density
888
889	!! 0.2 Output variables
890	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PSCAP !! Heat capacity (J K-1 m-3)
891
892	PSCAP(:) = PSNOWRHO(:)*xci
893
894	END FUNCTION snow3lscap_1d
895
896
897	!!
898	!================================================================================================================================
899	!! FUNCTION : snow3ltemp_2d
900	!!
901	!>\BRIEF Diagnose snow temperature (K) from heat content (J m-2)
902	!!
903	!! DESCRIPTION :
904	!!
905	!! RECENT CHANGE(S): None
906	!!
907	!! MAIN OUTPUT VARIABLE(S): ::
908	!!
909	!! REFERENCE(S) :
910	!!
911	!! FLOWCHART : None
912	!! \n
913	!_
914	!================================================================================================================================
915	FUNCTION snow3ltemp_2d(PSNOWHEAT,PSNOWRHO,PSNOWDZ) RESULT(PSNOWTEMP)
916
917	!! 0.1 Input variables
918	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO !! layer density (kg m-3)
919	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ !! layer thickness (m)
920	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWHEAT !! heat content (J m-2)
921
922	!! 0.2 Output variables
923	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PSNOWTEMP !! layer temperature (K)
924
925	!! 0.3 Modified variables
926
927	!! 0.4 Local variables
928	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSCAP !! snow heat capacity (J K-1 m-3)
929
930	ZSCAP(:,:) = snow3lscap_2d(PSNOWRHO)
931
932	PSNOWTEMP(:,:) = tp_00 + ( ((PSNOWHEAT(:,:)/PSNOWDZ(:,:)) &
933	+ chalfu0*PSNOWRHO(:,:))/ZSCAP(:,:) )
934
935	PSNOWTEMP(:,:) = MIN(tp_00, PSNOWTEMP(:,:))
936
937	END FUNCTION snow3ltemp_2d
938
939	!!
940	!================================================================================================================================
941	!! FUNCTION : snow3ltemp_1d
942	!!
943	!>\BRIEF Diagnose snow temperature (K) from heat content (J m-2)
944	!!
945	!! DESCRIPTION :
946	!!
947	!! RECENT CHANGE(S): None
948	!!
949	!! MAIN OUTPUT VARIABLE(S): ::
950	!!
951	!! REFERENCE(S) :
952	!!
953	!! FLOWCHART : None
954	!! \n
955	!_
956	!================================================================================================================================
957	FUNCTION snow3ltemp_1d(PSNOWHEAT,PSNOWRHO,PSNOWDZ) RESULT(PSNOWTEMP)
958
959	!! 0.1 Input variables
960	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! layer density (kg m-3)
961	REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ !! layer thickness (m)
962	REAL, DIMENSION(:), INTENT(IN) :: PSNOWHEAT !! heat content (J m-2)
963
964	!! 0.2 Output variables
965	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PSNOWTEMP !! layer temperature (K)
966
967	!! 0.3 Modified variables
968
969	!! 0.4 Local variables
970	REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZSCAP !! snow heat capacity (J K-1 m-3)
971
972	ZSCAP(:) = snow3lscap_1d(PSNOWRHO)
973
974	PSNOWTEMP(:) = tp_00 + ( ((PSNOWHEAT(:)/PSNOWDZ(:)) &
975	+ chalfu0*PSNOWRHO(:))/ZSCAP(:) )
976
977	PSNOWTEMP(:) = MIN(tp_00, PSNOWTEMP(:))
978	WHERE(PSNOWTEMP(:) .LE. 100) PSNOWTEMP(:) = tp_00
979
980	END FUNCTION snow3ltemp_1d
981
982	!!
983	!================================================================================================================================
984	!! FUNCTION : snow3lgrain_2d
985	!!
986	!>\BRIEF Calculate the grain size (m) for initialization
987	!!
988	!! DESCRIPTION :
989	!!
990	!! RECENT CHANGE(S): None
991	!!
992	!! MAIN OUTPUT VARIABLE(S): ::
993	!!
994	!! REFERENCE(S) : Loth and Graf 1993
995	!!
996	!! FLOWCHART : None
997	!! \n
998	!_
999	!================================================================================================================================
1000	FUNCTION snow3lgrain_2d(PSNOWRHO) RESULT(PDSGRAIN)
1001
1002	!! 0.1 Input variables
1003	REAL(r_std), DIMENSION(:,:), INTENT(IN) :: PSNOWRHO !! Snow density
1004
1005	!! 0.2 Output variables
1006	REAL(r_std), DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PDSGRAIN !! Snow grain size
1007
1008	!! 0.3 Modified variables
1009
1010	!! 0.4 Local variables
1011	REAL(r_std), PARAMETER :: ZSNOWRAD_AGRAIN = 1.6e-4 !! (m)
1012	REAL(r_std), PARAMETER :: ZSNOWRAD_BGRAIN = 1.1e-13 !! (m13/kg4)
1013	REAL(r_std), PARAMETER :: ZDSGRAIN_MAX = 2.796e-3 !! (m)
1014
1015	! grain size in m:
1016
1017	PDSGRAIN(:,:) = ZSNOWRAD_AGRAIN + ZSNOWRAD_BGRAIN(PSNOWRHO(:,:)*4)
1018	PDSGRAIN(:,:) = MIN(ZDSGRAIN_MAX, PDSGRAIN(:,:))
1019
1020	END FUNCTION snow3lgrain_2d
1021
1022	!!
1023	!================================================================================================================================
1024	!! FUNCTION : snow3lgrain_1d
1025	!!
1026	!>\BRIEF Calculate the grain size (m) for initialization
1027	!!
1028	!! DESCRIPTION :
1029	!!
1030	!! RECENT CHANGE(S): None
1031	!!
1032	!! MAIN OUTPUT VARIABLE(S): ::
1033	!!
1034	!! REFERENCE(S) : Loth and Graf 1993
1035	!!
1036	!! FLOWCHART : None
1037	!! \n
1038	!_
1039	!================================================================================================================================
1040	FUNCTION snow3lgrain_1d(PSNOWRHO) RESULT(PDSGRAIN)
1041
1042	!! 0.1 Input variables
1043	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! Snow density
1044
1045	!! 0.2 Output variables
1046	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PDSGRAIN !! Snow grain size
1047
1048	!! 0.3 Modified variables
1049
1050	!! 0.4 Local variables
1051	REAL, PARAMETER :: ZSNOWRAD_AGRAIN = 1.6e-4 !! (m)
1052	REAL, PARAMETER :: ZSNOWRAD_BGRAIN = 1.1e-13 !! (m13/kg4)
1053	REAL, PARAMETER :: ZDSGRAIN_MAX = 2.796e-3!! (m)
1054
1055	! grain size in m:
1056
1057	PDSGRAIN(:) = ZSNOWRAD_AGRAIN + ZSNOWRAD_BGRAIN(PSNOWRHO(:)*4)
1058	PDSGRAIN(:) = MIN(ZDSGRAIN_MAX, PDSGRAIN(:))
1059
1060	END FUNCTION snow3lgrain_1d
1061
1062	!================================================================================================================================
1063	!! FUNCTION : snow3lgrain_0d
1064	!!
1065	!>\BRIEF Calculate the grain size (m) for initialization
1066	!!
1067	!! DESCRIPTION :
1068	!!
1069	!! RECENT CHANGE(S): None
1070	!!
1071	!! MAIN OUTPUT VARIABLE(S): ::
1072	!!
1073	!! REFERENCE(S) : Loth and Graf 1993
1074	!!
1075	!! FLOWCHART : None
1076	!! \n
1077	!_
1078	!================================================================================================================================
1079	FUNCTION snow3lgrain_0d(PSNOWRHO) RESULT(PDSGRAIN)
1080
1081	!! 0.1 Input variables
1082	REAL(r_std), INTENT(IN) :: PSNOWRHO !! Snow density
1083
1084	!! 0.2 Output variables
1085	REAL(r_std) :: PDSGRAIN !! Snow grain size
1086
1087	!! 0.3 Modified variables
1088
1089	!! 0.4 Local variables
1090	REAL, PARAMETER :: ZSNOWRAD_AGRAIN = 1.6e-4 !! (m)
1091	REAL, PARAMETER :: ZSNOWRAD_BGRAIN = 1.1e-13 !! (m13/kg4)
1092	REAL, PARAMETER :: ZDSGRAIN_MAX = 2.796e-3!! (m)
1093
1094	! grain size in m:
1095
1096	PDSGRAIN = ZSNOWRAD_AGRAIN + ZSNOWRAD_BGRAIN(PSNOWRHO*4)
1097	PDSGRAIN = MIN(ZDSGRAIN_MAX, PDSGRAIN)
1098
1099	END FUNCTION snow3lgrain_0d
1100
1101	!================================================================================================================================
1102	!! FUNCTION : snow3lliq_2d
1103	!!
1104	!>\BRIEF Diagnose snow liquid water content from temperature (K) and
1105	!! heat content (J m-2)
1106	!!
1107	!! DESCRIPTION : Diagnose snow liquid water content from temperature (K)
1108	!! and heat content (J m-2). Note, need to evaluate SNOWTEMP from
1109	!! SNOW3LTEMP before calling this function (i.e. using same
1110	!! heat content, mass and diagnosed temperature).
1111	!!
1112	!! RECENT CHANGE(S): None
1113	!!
1114	!! MAIN OUTPUT VARIABLE(S): ::
1115	!!
1116	!! REFERENCE(S) :
1117	!!
1118	!! FLOWCHART : None
1119	!! \n
1120	!_
1121	!================================================================================================================================
1122	FUNCTION snow3lliq_2d(PSNOWHEAT,PSNOWRHO,PSNOWDZ,PSNOWTEMP)&
1123	& RESULT(PSNOWLIQ)
1124
1125	!! 0.1 Input variables
1126	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO &
1127	& !! layer density (kg m-3)
1128	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ &
1129	& !! layer thickness (m)
1130	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWHEAT&
1131	& !! heat content (J m-2)
1132	REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWTEMP&
1133	& !! layer temperature (K)
1134
1135	!! 0.2 Output variables
1136	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: PSNOWLIQ &
1137	& ! liquid water content (m)
1138
1139	!! 0.3 Modified variables
1140
1141	!! 0.4 Local variables
1142	REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSCAP &
1143	& !! snow heat capacity (J K-1 m-3)
1144
1145	ZSCAP(:,:) = snow3lscap_2d(PSNOWRHO)
1146
1147	! The result of the full heat balance equation: if the sum
1148	! equals zero,
1149	! then no liquid. If an imbalance occurs, this represents
1150	! liquid water content.
1151
1152	PSNOWLIQ(:,:) = ( ((tp_00-PSNOWTEMP(:,:))*ZSCAP(:,:) + chalfu0&
1153	&PSNOWRHO(:,:))PSNOWDZ(:,:) + PSNOWHEAT(:,:) ) /(chalfu0&
1154	&*ph2o)
1155
1156	! just a numerical check:
1157
1158	PSNOWLIQ(:,:) = MAX(0.0, PSNOWLIQ(:,:))
1159
1160	END FUNCTION snow3lliq_2d
1161
1162	!================================================================================================================================
1163	!! FUNCTION : snow3lliq_1d
1164	!!
1165	!>\BRIEF Diagnose snow liquid water content from temperature (K) and
1166	!! heat content (J m-2)
1167	!!
1168	!! DESCRIPTION : Diagnose snow liquid water content from temperature (K)
1169	!! and heat content (J m-2). Note, need to evaluate SNOWTEMP from
1170	!! SNOW3LTEMP before calling this function (i.e. using same
1171	!! heat content, mass and diagnosed temperature).
1172	!!
1173	!! RECENT CHANGE(S): None
1174	!!
1175	!! MAIN OUTPUT VARIABLE(S): ::
1176	!!
1177	!! REFERENCE(S) :
1178	!!
1179	!! FLOWCHART : None
1180	!! \n
1181	!_
1182	!================================================================================================================================
1183	FUNCTION snow3lliq_1d(PSNOWHEAT,PSNOWRHO,PSNOWDZ,PSNOWTEMP) RESULT(PSNOWLIQ)
1184
1185	!! 0.1 Input variables
1186	REAL, DIMENSION(:), INTENT(IN) :: PSNOWRHO !! layer density (kg m-3)
1187	REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ !! layer thickness (m)
1188	REAL, DIMENSION(:), INTENT(IN) :: PSNOWHEAT !! heat content (J m-2)
1189	REAL, DIMENSION(:), INTENT(IN) :: PSNOWTEMP !! layer temperature (K)
1190
1191	!! 0.2 Output variables
1192	REAL, DIMENSION(SIZE(PSNOWRHO)) :: PSNOWLIQ !! liquid water content (m)
1193
1194	!! 0.3 Modified variables
1195
1196	!! 0.4 Local variables
1197	REAL, DIMENSION(SIZE(PSNOWRHO)) :: ZSCAP !! snow heat capacity (J K-1 m-3)
1198
1199	ZSCAP(:) = snow3lscap_1d(PSNOWRHO)
1200
1201	! The result of the full heat balance equation: if the sum equals zero,
1202	! then no liquid. If an imbalance occurs, this represents liquid water content.
1203	!
1204	PSNOWLIQ(:) = ( ((tp_00-PSNOWTEMP(:))*ZSCAP(:) + &
1205	chalfu0PSNOWRHO(:))PSNOWDZ(:) + PSNOWHEAT(:) ) &
1206	/(chalfu0*ph2o)
1207
1208	! just a numerical check:
1209
1210	PSNOWLIQ(:) = MAX(0.0, PSNOWLIQ(:))
1211
1212	END FUNCTION snow3lliq_1d
1213
1214	END MODULE qsat_moisture

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

Download in other formats: