source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/qsat_moisture.f90

Last change on this file was 3924, checked in by albert.jornet, 8 years ago

Merge: from trunk revisions [3643:3883/trunk/ORCHIDEE]

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