source: branches/publications/ORCHILEAK-Gommet/src_sechiba/qsat_moisture.f90 @ 7746

Last change on this file since 7746 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
24MODULE 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)-deux*zz_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-deux*zz_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 = zrapp*10.0_r_std**(2.07023_r_std-zcorr*ztemperature &
563                 &             -2484.896/ztemperature+3.56654*LOG10(ztemperature)) ! Equilibrium water vapor - solid
564         ELSE
565            zqsat = zrapp*10.0**(23.8319-2948.964/ztemperature & 
566                 &              -5.028*LOG10(ztemperature) &
567                 &              -29810.16*EXP(-0.0699382*ztemperature) &
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 = ZHOLDMAXR*PSNOWDZ*ZSNOWRHO/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         - chalfu0*PSNOWRHO(:,:) ) + chalfu0*ph2o*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         -chalfu0*PSNOWRHO(:) ) + chalfu0*ph2o*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         chalfu0*PSNOWRHO(:))*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 
1198END MODULE qsat_moisture
Note: See TracBrowser for help on using the repository browser.