Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_sechiba/qsat_moisture.f90 @ 7346

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

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

Download in other formats: