source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/qsat_moisture.f90 @ 6268

Last change on this file since 6268 was 737, checked in by didier.solyga, 12 years ago

Correct units for bqsb and gqsb. Improve documentation
  • Property svn:keywords set to Date Revision
File size: 19.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
30!-
31  IMPLICIT NONE
32!-
33
34  LOGICAL,SAVE :: l_qsat_first=.TRUE.                !! First call to qsat subroutines and functions (true/false)
35
36
37  INTEGER(i_std),PARAMETER :: max_temp=370           !! Maximum temperature for saturated humidity (K) and also used as
38                                                     !! the size of local array to keep saturated humidity (unitless)
39
40  INTEGER(i_std),PARAMETER :: min_temp=100           !! Minimum temperature for saturated humidity (K)
41
42  REAL(r_std),DIMENSION(max_temp),SAVE :: qsfrict    !! Array to keep saturated humidity at each temperature level
43                                                     !! (kg of water/kg of air)
44
45  CONTAINS
46
47!! ================================================================================================================================
48!! SUBROUTINE   : qsatcalc
49!!
50!>\BRIEF          This routine calculates the saturated humidity using the pressure
51!!                and the temperature for all pixels.
52!!
53!! DESCRIPTION : This routine interpolates qsat between temperatures by the following formula :
54!!               \latexonly
55!!               \input{qsatcalc.tex}
56!!               \endlatexonly
57!!               \n
58!!
59!! RECENT CHANGE(S): None
60!!
61!! MAIN OUTPUT VARIABLE(S) : qsat_out
62!!
63!! REFERENCE(S) : None
64!!
65!! FLOWCHART    : None
66!! \n
67!_ ================================================================================================================================
68
69    SUBROUTINE qsatcalc (kjpindex,temp_in,pres_in,qsat_out)
70
71      IMPLICIT NONE
72
73      !! 0. Variables and parameters declaration
74     
75      !! 0.1 Input variables
76     
77      INTEGER(i_std),INTENT(in) :: kjpindex                   !! Domain size (unitless)
78      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: temp_in  !! Temperature in degre Kelvin (K)
79      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: pres_in  !! Pressure (Pa)
80     
81      !! 0.2 Output variables
82     
83      REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: qsat_out !! Saturated humidity at the surface (kg of water/kg of air)
84     
85      !! 0.4 Local variables
86     
87      INTEGER(i_std), DIMENSION(kjpindex) :: jt               !! Temporary array stocking the truncated temperatures in Kelvin
88                                                              !!(converted into integers)
89      INTEGER(i_std)                      :: ji               !! indices (unitless)
90      REAL(r_std),DIMENSION(kjpindex)     :: zz_a, zz_b, zz_f !! Temporary variables
91      INTEGER(i_std)                      :: nbad             !! Number of points where the temperature is too high or too low
92      INTEGER(i_std),DIMENSION(1)         :: lo               !! Temporary vector to mark the position of the highest temperature
93                                                              !! or the lowest temperature over all the pixels in jt (unitless)
94!_ ================================================================================================================================
95
96      !-
97      !! 1.Initialize qsfrict array if needed
98      !-
99      IF (l_qsat_first) THEN
100         !-
101         CALL qsfrict_init
102         l_qsat_first = .FALSE.
103         !-
104      ENDIF !(l_qsat_first)
105     
106      !-
107      !! 2. Computes qsat interpolation into two successive temperature
108      !-
109      jt = INT(temp_in(:))
110       
111      !! 2.1 Diagnostic pixels where the temperature is too high
112      nbad = COUNT(jt(:) >= max_temp-1)
113     
114      IF (nbad > 0) THEN
115         WRITE(numout,*) ' qsatcalc: temperature too high at ', &
116              &    nbad, ' points.'
117         !-
118         IF (.NOT.diag_qsat) THEN
119            CALL ipslerr(2,'qsatcalc','diffuco', '', &
120                 &                     'temperature incorect.') ! Warning message
121         ELSE
122            lo = MAXLOC(temp_in(:))
123            WRITE(numout,*) &
124                 &     'Maximum temperature ( ',MAXVAL(temp_in),') found at ',lo(1)
125            WHERE (jt(:) >= max_temp-1)   jt(:) = max_temp-1
126         ENDIF !(.NOT.diag_qsat)
127         !-
128      ENDIF ! (nbad > 0)
129 
130     
131      !! 2.2 Diagnostic pixels where the temperature is too low
132      nbad = COUNT(jt(:) <= min_temp)
133     
134      IF (nbad > 0) THEN
135         WRITE(numout,*) ' qsatcalc: temperature too low at ', &
136              &    nbad, ' points.'
137         !-
138         IF (.NOT.diag_qsat) THEN
139            CALL ipslerr(2,'qsatcalc','diffuco', '', &
140                 &                   'temperature incorect.') ! Warning message
141         ELSE
142            lo = MINLOC(temp_in(:))
143            WRITE(numout,*) &
144                 &     'Minimum temperature ( ',MINVAL(temp_in),') found at ',lo(1)
145            WHERE (jt(:) <= min_temp)   jt(:) = min_temp
146         ENDIF !(.NOT.diag_qsat)
147         !-
148      ENDIF! (nbad > 0)
149 
150      !! 2.3 Temporary variables needed for interpolation
151      DO ji = 1, kjpindex ! Loop over # pixels
152         
153         zz_f(ji) = temp_in(ji)-FLOAT(jt(ji))
154         zz_a(ji) = qsfrict(jt(ji))
155         zz_b(ji) = qsfrict(jt(ji)+1)
156         
157      ENDDO ! Loop over # pixels
158     
159      !-
160      !! 3. Interpolation between these two values
161      !-
162      DO ji = 1, kjpindex ! Loop over # pixels
163         
164         qsat_out(ji) = ((zz_b(ji)-zz_a(ji))*zz_f(ji)+zz_a(ji))/pres_in(ji)
165         
166      ENDDO  ! Loop over # pixels
167
168
169    END SUBROUTINE qsatcalc
170 
171!! ================================================================================================================================
172!! FUNCTION     : [DISPENSABLE] qsat
173!!
174!>\BRIEF         This function computes deviation the saturated humidity with the pressure
175!! and the temperature for a scalar.
176!!
177!! DESCRIPTION : This routine is obsolete : replaced by the subroutine qsatcalc. \n
178!!               qsat is interpolated by : \n
179!!               \latexonly
180!!               \input{qsat.tex}
181!!               \endlatexonly
182!!
183!! RECENT CHANGE(S): None\n
184!!
185!! RETURN VALUE : qsat_result
186!!
187!! REFERENCE(S) : None
188!!
189!! FLOWCHART    : None
190!! \n
191!_ ================================================================================================================================   
192
193    FUNCTION qsat (temp_in,pres_in) RESULT (qsat_result)
194
195      IMPLICIT NONE
196
197      !! 0. Variables and parameters declaration
198
199      !! 0.1 Input variables
200
201      REAL(r_std),INTENT(in) :: temp_in   !! Temperature (K)
202      REAL(r_std),INTENT(in) :: pres_in   !! Pressure    (Pa)
203
204      !! 0.2 Result
205     
206      REAL(r_std) :: qsat_result          !! Saturated humidity calculated at the surface (kg/kg)
207     
208      !! 0.4 Local variables
209     
210      INTEGER(i_std)   :: jt              !! Temporary scalar stocking the truncated temperature in Kelvin
211                                          !! (converted into integer)
212      REAL(r_std)      :: zz_a,zz_b,zz_f  !! Temporary scalar variables     
213 
214!_ ================================================================================================================================
215
216      !-
217      !! 1.Initialize qsfrict array if needed
218      !-
219      IF (l_qsat_first) THEN
220         !-
221         CALL qsfrict_init
222         l_qsat_first = .FALSE.
223         !-
224      ENDIF
225
226      !-
227      !! 2. Computes qsat interpolation into two successive temperatures
228      !-
229      jt = INT(temp_in)
230
231      !! 2.1 Is the temperature too high ?
232      IF (jt >= max_temp-1) THEN
233         WRITE(numout,*) &
234              &   ' We stop. temperature too BIG : ',temp_in, &
235              &   ' approximation for : ',jt
236         !-
237         IF (.NOT.diag_qsat) THEN
238            CALL ipslerr(2,'qsat','', '',&
239                 &                   'temperature incorect.') ! Warning message
240         ELSE
241            qsat_result = 999999.
242            RETURN
243         ENDIF !(.NOT.diag_qsat)
244         !-
245      ENDIF !(jt >= max_temp-1)
246
247      !! 2.2 Is the temperature too low ?
248      IF (jt <= min_temp ) THEN
249         WRITE(numout,*) &
250              &   ' We stop. temperature too SMALL : ',temp_in, &
251              &   ' approximation for : ',jt
252         !-
253         IF (.NOT.diag_qsat) THEN
254            CALL ipslerr(2,'qsat','', '',&
255                 &                   'temperature incorect.')
256         ELSE
257            qsat_result = -999999.
258            RETURN
259         ENDIF!(.NOT.diag_qsat)
260         !-
261      ENDIF !(jt <= min_temp )
262
263      !! 2.3 Temporary variables needed for interpolation
264      zz_f = temp_in-FLOAT(jt)
265      zz_a = qsfrict(jt)
266      zz_b = qsfrict(jt+1)
267
268      !! 3. Interpolates between these two values
269
270      qsat_result = ((zz_b-zz_a)*zz_f+zz_a)/pres_in
271
272
273    END FUNCTION qsat
274
275
276!! ================================================================================================================================
277!! SUBROUTINE   : dev_qsatcalc
278!!
279!>\BRIEF         This routine calculates the deviation of the saturated humidity qsat.
280!!
281!! DESCRIPTION : The deviation of qsat is calculated by :
282!!               \latexonly
283!!               \input{dev_qsatcalc.tex}
284!!               \endlatexonly
285!!
286!! RECENT CHANGE(S): None
287!!
288!! MAIN OUTPUT VARIABLE(S) : dev_qsat_out
289!!
290!! REFERENCE(S) : None
291!!
292!! FLOWCHART    : None
293!!
294!! FLOWCHART    :
295!! \latexonly
296!! \includegraphics[scale = 1]{pheno_moigdd.png}
297!! \endlatexonly
298!! \n
299!_ ================================================================================================================================
300
301    SUBROUTINE dev_qsatcalc (kjpindex,temp_in,pres_in,dev_qsat_out)
302
303      IMPLICIT NONE
304
305      !! 0. Variables and parameters declaration
306     
307      !! 0.1 Input variables
308
309      INTEGER(i_std),INTENT(in)                   :: kjpindex       !! Domain size (unitless)
310      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: temp_in        !! Temperature (K)
311      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: pres_in        !! Pressure
312
313
314      !! 0.2 Output variables
315     
316      REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: dev_qsat_out   !! Result (??units??)
317
318     
319      !! 0.4 Local variables
320     
321      INTEGER(i_std),DIMENSION(kjpindex) :: jt                      !! Temporary array stocking the truncated temperatures
322                                                                    !! in Kelvin (converted into integers)
323      INTEGER(i_std)                     :: ji                      !! Indice (unitless)
324      REAL(r_std),DIMENSION(kjpindex)    :: zz_a, zz_b, zz_c, zz_f  !! Temporary vector variables
325      INTEGER(i_std)                     :: nbad                    !! Number of points where the temperature is too high or too low
326
327!_ ================================================================================================================================
328
329      !-
330      !! 1.Initialize qsfrict array if needed
331      !-
332      IF (l_qsat_first) THEN
333         !-
334         CALL qsfrict_init
335         l_qsat_first = .FALSE.
336         !-
337      ENDIF
338
339      !-
340      !! 2. Compute qsat interpolation into two successive temperature
341      !-
342      jt = INT(temp_in(:)+undemi)
343     
344      !! 2.1 Pixels where the temperature is too high
345      nbad = COUNT( jt(:) >= max_temp-1 )
346     
347      IF (nbad > 0) THEN
348         WRITE(numout,*) &
349              &   ' dev_qsatcalc: temperature too high at ',nbad,' points.'
350         !-
351         IF (.NOT.diag_qsat) THEN
352            CALL ipslerr(3,'dev_qsatcalc','', '', &
353                 &                   'temperature incorect.') ! Fatal error
354         ELSE
355            WHERE (jt(:) >= max_temp-1)   jt(:) = max_temp-1
356         ENDIF !(.NOT.diag_qsat)
357         !-
358      ENDIF !(nbad > 0)
359     
360      !! 2.2 Pixels where the temperature is too low
361      nbad = COUNT( jt(:) <= min_temp )
362     
363      IF (nbad > 0) THEN
364         WRITE(numout,*) &
365              &   ' dev_qsatcalc: temperature too low at ',nbad,' points.'
366         !-
367         IF (.NOT.diag_qsat) THEN
368            CALL ipslerr(3,'dev_qsatcalc', '', '',&
369                 &                   'temperature incorect.') ! Fatal error
370         ELSE
371            WHERE (jt(:) <= min_temp)   jt(:) = min_temp
372         ENDIF !(.NOT.diag_qsat)
373         !-
374      ENDIF !(nbad > 0)
375
376      !! 2.3 Temporary variables needed for interpolation
377      DO ji=1,kjpindex  ! Loop over # pixels
378         
379         zz_f(ji) = temp_in(ji)+undemi-FLOAT(jt(ji))
380         zz_a(ji) = qsfrict(jt(ji)-1)
381         zz_b(ji) = qsfrict(jt(ji))
382         zz_c(ji) = qsfrict(jt(ji)+1)
383         
384      ENDDO ! Loop over # pixels
385     
386      !-
387      !! 3. Interpolates between these two values
388      !-
389      DO ji = 1, kjpindex ! Loop over # pixels
390         
391         dev_qsat_out(ji) = &
392              &   ((zz_c(ji)-deux*zz_b(ji)+zz_a(ji))*(zz_f(ji)-un) + &
393              &                         zz_c(ji)-zz_b(ji))/pres_in(ji)
394         
395      ENDDO ! Loop over # pixels
396
397
398    END SUBROUTINE dev_qsatcalc
399
400!! ================================================================================================================================
401!! FUNCTION     : [DISPENSABLE] dev_qsat
402!!
403!>\BRIEF         This function computes deviation of qsat.
404!!
405!! DESCRIPTION : The deviation of qsat is calculated by :
406!!               \latexonly
407!!               \input{dev_qsat.tex}
408!!               \endlatexonly
409!!
410!! RECENT CHANGE(S): None
411!!
412!! RETURN VALUE : dev_qsat_result
413!!
414!! REFERENCE(S) : None
415!!
416!! FLOWCHART    : None
417!! \n
418!_ ================================================================================================================================
419
420    FUNCTION dev_qsat (temp_in,pres_in) RESULT (dev_qsat_result)
421
422      IMPLICIT NONE
423
424      !! 0. Variables and parameters declaration
425   
426      !! 0.1 Input variables
427
428      REAL(r_std),INTENT(in)  :: pres_in          !! Pressure (Pa)
429      REAL(r_std),INTENT(in)  :: temp_in          !! Temperture (K)
430     
431      !! 0.2 Result
432
433      REAL(r_std) :: dev_qsat_result              !! (??units??) !!
434     
435      !! 0.4 Local variables
436
437      INTEGER(i_std)  :: jt                       !! Index (unitless)
438      REAL(r_std)     :: zz_a, zz_b, zz_c, zz_f   !! Temporary scalars   
439
440!_ ================================================================================================================================
441 
442      !-
443      !! 1.Initialize qsfrict array if needed
444      !-
445      IF (l_qsat_first) THEN
446         !-
447         CALL qsfrict_init
448         l_qsat_first = .FALSE.
449         !-
450      ENDIF
451
452      !-
453      !! 2. computes qsat deviation interpolation
454      !!    into two successive temperature
455      !-
456      jt = INT(temp_in+undemi)
457
458      !! 2.1 Is the temperature too high ?
459      IF (jt >= max_temp-1) THEN
460         !-
461         WRITE(numout,*) &
462              &   ' We stop. temperature too HIGH : ',temp_in, &
463              &   ' approximation for : ',jt
464         IF (.NOT.diag_qsat) THEN
465            CALL ipslerr(3,'dev_qsat','', '',&
466                 &                   'temperature incorect.') ! Fatal error
467         ELSE
468            dev_qsat_result = 999999.
469            RETURN
470         ENDIF !(.NOT.diag_qsat)
471         !-
472      ENDIF !(jt >= max_temp-1)
473      !-
474      !! 2.2 Is the temperature too low ?
475      IF (jt <= min_temp ) THEN
476         WRITE(numout,*) &
477              &   ' We stop. temperature too LOW : ',temp_in, &
478              &   ' approximation for : ',jt
479         !-
480         IF (.NOT.diag_qsat) THEN
481            CALL ipslerr(3,'dev_qsat','', '',&
482                 &                    'temperature incorect.')
483         ELSE
484            dev_qsat_result = -999999.
485            RETURN
486         ENDIF !(.NOT.diag_qsat)
487         !-
488      ENDIF !(jt <= min_temp )
489
490      !! 2.3  Temporary variables for interpolation
491      zz_f = temp_in+undemi-FLOAT(jt)
492      zz_a = qsfrict(jt-1)
493      zz_b = qsfrict(jt)
494      zz_c = qsfrict(jt+1)
495
496      !-
497      !! 3. Interpolate
498      !-
499      dev_qsat_result=((zz_c-deux*zz_b+zz_a)*(zz_f-un)+zz_c-zz_b)/pres_in
500
501    END FUNCTION dev_qsat
502 
503
504!! ================================================================================================================================
505!! SUBROUTINE   : qsfrict_init
506!!
507!>\BRIEF         The qsfrict_init routine initialises qsfrict array to store
508!!               precalculated values for qsat by using Goff-Gratch equations. 
509!!
510!! DESCRIPTION : This routine calculates the specific humidity qsat as a function of temperature in
511!!               Kelvin by using the modified Goff-Gratch equations(1946): \n
512!!               \latexonly
513!!               \input{goff_gratch.tex}
514!!               \endlatexonly
515!!               qsfrict is initialized by the following formulas : \n
516!!               \latexonly
517!!               \input{qsfrict_init.tex}
518!!               \endlatexonly
519!!               These values are used by the subroutines qsatcalc, dev_qsat. \n
520!!
521!! RECENT CHANGE(S): None
522!!
523!! MAIN OUTPUT VARIABLE(S): ::qsfrict
524!!
525!! REFERENCE(S) :
526!! - Algorithme d'un ensemble de paramétrisation physique (1998),
527!! Note de Laurent Li décrivant les paramétrisations physiques incluses dans le modÚle (pdf),
528!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques
529!! - Goff, J. A., and S. Gratch (1946) Low-pressure properties of water from −160 to 212 °F, in Transactions of the
530!! American Society of Heating and Ventilating Engineers, pp 95–122, presented at the 52nd annual meeting of the
531!! American Society of Heating and Ventilating Engineers, New York, 1946.
532!!
533!! FLOWCHART    : None
534!! \n
535!_ ================================================================================================================================
536
537    SUBROUTINE qsfrict_init
538     
539      IMPLICIT NONE
540
541      !! 0. Variables and parameters declaration
542     
543      !! 0.4 Local variables
544
545      INTEGER(i_std) :: ji                             !! Indice(unitless)
546      REAL(r_std)    :: zrapp,zcorr,ztemperature,zqsat !! Temporary vector variables     
547
548!_ ================================================================================================================================
549
550      !! 1. Initialisation
551      zrapp = msmlr_h2o/msmlr_air
552      zcorr = 0.00320991_r_std
553     
554      !! 2. Computes saturated humidity one time and store in qsfrict local array
555      DO ji=100,max_temp ! Loop over size(qsfrict) : each position of qsfrict matches a temperature
556         
557         ztemperature = FLOAT(ji)
558         !-
559         IF (ztemperature < 273._r_std) THEN
560            zqsat = zrapp*10.0_r_std**(2.07023_r_std-zcorr*ztemperature &
561                 &             -2484.896/ztemperature+3.56654*LOG10(ztemperature)) ! Equilibrium water vapor - solid
562         ELSE
563            zqsat = zrapp*10.0**(23.8319-2948.964/ztemperature & 
564                 &              -5.028*LOG10(ztemperature) &
565                 &              -29810.16*EXP(-0.0699382*ztemperature) &
566                 &              +25.21935*EXP(-2999.924/ztemperature)) ! Equilibrium water vapor - liquid
567         ENDIF !(ztemperature < 273._r_std)
568         !-
569         qsfrict (ji) = zqsat
570         
571      ENDDO ! Loop over size(qsfrict)
572
573      !! 3. Set to zero the non-computed values
574      qsfrict(1:100) = zero 
575      !-
576      IF (long_print) WRITE (numout,*) ' qsfrict_init done'
577
578
579    END SUBROUTINE qsfrict_init
580
581
582END MODULE qsat_moisture
Note: See TracBrowser for help on using the repository browser.