source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/qsat_moisture.f90 @ 277

Last change on this file since 277 was 64, checked in by didier.solyga, 14 years ago

Import first version of ORCHIDEE_EXT
File size: 10.7 KB
Line 
1! 29/11/2010 : New module used for computing qsat
2! previously in constantes_veg
3! We decide to put these routines in sechiba
4
5MODULE qsat_moisture
6
7  USE defprec
8  USE constantes
9  USE IOIPSL
10
11!-
12  IMPLICIT NONE
13
14  !-
15  LOGICAL,SAVE :: l_qsat_first=.TRUE.
16
17  ! Size of local array to keep saturated humidity
18  ! at each temperature level
19  INTEGER(i_std),PARAMETER :: max_temp=370
20
21  ! Minimum temperature for saturated humidity
22  INTEGER(i_std),PARAMETER :: min_temp=100
23
24  ! Local array to keep saturated humidity at each temperature level
25  REAL(r_std),DIMENSION(max_temp),SAVE :: qsfrict
26
27  CONTAINS
28!===
29    SUBROUTINE qsatcalc (kjpindex,temp_in,pres_in,qsat_out)
30!---------------------------------------------------------------------
31      ! input value
32      ! Domain size
33      INTEGER(i_std),INTENT(in) :: kjpindex
34      ! Temperature in degre Kelvin
35      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: temp_in
36      ! Pressure
37      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: pres_in
38    ! output value
39      ! Result
40      REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: qsat_out
41      !-
42    ! local variables
43      INTEGER(i_std), DIMENSION(kjpindex)       :: jt
44      INTEGER(i_std)                            :: ji
45      REAL(r_std),DIMENSION(kjpindex)     :: zz_a, zz_b, zz_f
46      INTEGER(i_std)                     :: nbad
47      INTEGER(i_std),DIMENSION(1)        :: lo
48!---------------------------------------------------------------------
49      IF (l_qsat_first) THEN
50         CALL qsfrict_init
51         l_qsat_first = .FALSE.
52      ENDIF
53      !-
54      ! 1. computes qsat interpolation into two successive temperature
55      !-
56      jt = INT(temp_in(:))
57      !-
58      nbad = COUNT(jt(:) >= max_temp-1)
59      IF (nbad > 0) THEN
60         WRITE(numout,*) ' qsatcalc: temperature too high at ', &
61              &    nbad, ' points.'
62         IF (.NOT.diag_qsat) THEN
63            CALL ipslerr(2,'qsatcalc','diffuco', '', &
64 &                     'temperature incorect.')
65         ELSE
66            lo = MAXLOC(temp_in(:))
67            WRITE(numout,*) &
68                 &     'Maximum temperature ( ',MAXVAL(temp_in),') found at ',lo(1)
69            WHERE (jt(:) >= max_temp-1)   jt(:) = max_temp-1
70         ENDIF
71      ENDIF
72      !-
73      nbad = COUNT(jt(:) <= min_temp)
74      IF (nbad > 0) THEN
75         WRITE(numout,*) ' qsatcalc: temperature too low at ', &
76              &    nbad, ' points.'
77         IF (.NOT.diag_qsat) THEN
78            CALL ipslerr(2,'qsatcalc','diffuco', '', &
79                 &                   'temperature incorect.')
80         ELSE
81            lo = MINLOC(temp_in(:))
82            WRITE(numout,*) &
83                 &     'Minimum temperature ( ',MINVAL(temp_in),') found at ',lo(1)
84            WHERE (jt(:) <= min_temp)   jt(:) = min_temp
85         ENDIF
86      ENDIF
87      !-
88      DO ji = 1, kjpindex
89         zz_f(ji) = temp_in(ji)-FLOAT(jt(ji))
90         zz_a(ji) = qsfrict(jt(ji))
91         zz_b(ji) = qsfrict(jt(ji)+1)
92      ENDDO
93      !-
94      ! 2. interpolates between this two values
95      !-
96      DO ji = 1, kjpindex
97         qsat_out(ji) = ((zz_b(ji)-zz_a(ji))*zz_f(ji)+zz_a(ji))/pres_in(ji)
98      ENDDO
99      !----------------------
100    END SUBROUTINE qsatcalc
101    !===
102    FUNCTION qsat (temp_in,pres_in) RESULT (qsat_result)
103      !!--------------------------------------------------------------------
104      !! FUNCTION qsat (temp_in, pres_in) RESULT (qsat_result)
105      !!--------------------------------------------------------------------
106      REAL(r_std),INTENT(in) :: temp_in  ! Temperature in degre Kelvin
107      REAL(r_std),INTENT(in) :: pres_in  ! Pressure
108      REAL(r_std) :: qsat_result
109      !-
110      INTEGER(i_std)        :: jt
111      REAL(r_std)     :: zz_a,zz_b,zz_f
112      !---------------------------------------------------------------------
113      IF (l_qsat_first) THEN
114         CALL qsfrict_init
115         l_qsat_first = .FALSE.
116      ENDIF
117      !-
118      ! 1. computes qsat interpolation into two successive temperature
119      !-
120      jt = INT(temp_in)
121      !-
122      IF (jt >= max_temp-1) THEN
123         WRITE(numout,*) &
124              &   ' We stop. temperature too BIG : ',temp_in, &
125              &   ' approximation for : ',jt
126         IF (.NOT.diag_qsat) THEN
127            CALL ipslerr(2,'qsat','', '',&
128                 &                   'temperature incorect.')
129         ELSE
130            qsat_result = 999999.
131            RETURN
132         ENDIF
133      ENDIF
134      !-
135      IF (jt <= min_temp ) THEN
136         WRITE(numout,*) &
137              &   ' We stop. temperature too SMALL : ',temp_in, &
138              &   ' approximation for : ',jt
139         IF (.NOT.diag_qsat) THEN
140            CALL ipslerr(2,'qsat','', '',&
141                 &                   'temperature incorect.')
142         ELSE
143            qsat_result = -999999.
144            RETURN
145         ENDIF
146      ENDIF
147      !-
148      zz_f = temp_in-FLOAT(jt)
149      zz_a = qsfrict(jt)
150      zz_b = qsfrict(jt+1)
151      !-
152      ! 2. interpolates between this two values
153      !-
154      qsat_result = ((zz_b-zz_a)*zz_f+zz_a)/pres_in
155      !----------------
156    END FUNCTION qsat
157    !===
158    SUBROUTINE dev_qsatcalc (kjpindex,temp_in,pres_in,dev_qsat_out)
159      !---------------------------------------------------------------------
160      ! Domain size
161      INTEGER(i_std),INTENT(in)                  :: kjpindex
162      ! Temperature in degre Kelvin
163      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: temp_in
164      ! Pressure
165      REAL(r_std),DIMENSION(kjpindex),INTENT(in)  :: pres_in
166      ! Result
167      REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: dev_qsat_out
168      !-
169      INTEGER(i_std),DIMENSION(kjpindex) :: jt
170      INTEGER(i_std)                     :: ji
171      REAL(r_std),DIMENSION(kjpindex)     :: zz_a, zz_b, zz_c, zz_f
172      INTEGER(i_std)                     :: nbad
173      !---------------------------------------------------------------------
174      IF (l_qsat_first) THEN
175         CALL qsfrict_init
176         l_qsat_first = .FALSE.
177      ENDIF
178      !-
179      ! 1. computes qsat interpolation into two successive temperature
180      !-
181      jt = INT(temp_in(:)+undemi)
182      !-
183      nbad = COUNT( jt(:) >= max_temp-1 )
184      IF (nbad > 0) THEN
185         WRITE(numout,*) &
186              &   ' dev_qsatcalc: temperature too high at ',nbad,' points.'
187         IF (.NOT.diag_qsat) THEN
188            CALL ipslerr(3,'dev_qsatcalc','', '', &
189                 &                   'temperature incorect.')
190         ELSE
191            WHERE (jt(:) >= max_temp-1)   jt(:) = max_temp-1
192         ENDIF
193      ENDIF
194      !-
195      nbad = COUNT( jt(:) <= min_temp )
196      IF (nbad > 0) THEN
197         WRITE(numout,*) &
198              &   ' dev_qsatcalc: temperature too low at ',nbad,' points.'
199         IF (.NOT.diag_qsat) THEN
200            CALL ipslerr(3,'dev_qsatcalc', '', '',&
201                 &                   'temperature incorect.')
202         ELSE
203            WHERE (jt(:) <= min_temp)   jt(:) = min_temp
204         ENDIF
205      ENDIF
206      !-
207      DO ji=1,kjpindex
208         zz_f(ji) = temp_in(ji)+undemi-FLOAT(jt(ji))
209         zz_a(ji) = qsfrict(jt(ji)-1)
210         zz_b(ji) = qsfrict(jt(ji))
211         zz_c(ji) = qsfrict(jt(ji)+1)
212      ENDDO
213      !-
214      ! 2. interpolates between this two values
215      !-
216      DO ji = 1, kjpindex
217         dev_qsat_out(ji) = &
218              &   ((zz_c(ji)-deux*zz_b(ji)+zz_a(ji))*(zz_f(ji)-un) + &
219              &                         zz_c(ji)-zz_b(ji))/pres_in(ji)
220      ENDDO
221      !--------------------------
222    END SUBROUTINE dev_qsatcalc
223    !===
224    FUNCTION dev_qsat (temp_in,pres_in) RESULT (dev_qsat_result)
225      !!--------------------------------------------------------------------
226      !! FUNCTION dev_qsat (temp_in, pres_in) RESULT (dev_qsat_result)
227      !! computes deviation of qsat
228      !!--------------------------------------------------------------------
229      REAL(r_std),INTENT(in)  :: pres_in    ! Pressure
230      REAL(r_std),INTENT(in)  :: temp_in    ! Temperture in degre Kelvin
231      REAL(r_std) :: dev_qsat_result
232      !-
233      INTEGER(i_std)        :: jt
234      REAL(r_std)     :: zz_a, zz_b, zz_c, zz_f
235      !---------------------------------------------------------------------
236      IF (l_qsat_first) THEN
237         CALL qsfrict_init
238         l_qsat_first = .FALSE.
239      ENDIF
240      !-
241      ! 1. computes qsat deviation interpolation
242      !    into two successive temperature
243      !-
244      jt = INT(temp_in+undemi)
245      !-
246      IF (jt >= max_temp-1) THEN
247         WRITE(numout,*) &
248              &   ' We stop. temperature too HIGH : ',temp_in, &
249              &   ' approximation for : ',jt
250         IF (.NOT.diag_qsat) THEN
251            CALL ipslerr(3,'dev_qsat','', '',&
252                 &                   'temperature incorect.')
253         ELSE
254            dev_qsat_result = 999999.
255            RETURN
256         ENDIF
257      ENDIF
258      !-
259      IF (jt <= min_temp ) THEN
260         WRITE(numout,*) &
261              &   ' We stop. temperature too LOW : ',temp_in, &
262              &   ' approximation for : ',jt
263         IF (.NOT.diag_qsat) THEN
264            CALL ipslerr(3,'dev_qsat','', '',&
265                 &                    'temperature incorect.')
266         ELSE
267            dev_qsat_result = -999999.
268            RETURN
269         ENDIF
270      ENDIF
271      !-
272      zz_f = temp_in+undemi-FLOAT(jt)
273      zz_a = qsfrict(jt-1)
274      zz_b = qsfrict(jt)
275      zz_c = qsfrict(jt+1)
276      !-
277      ! 2. interpolates
278      !-
279      dev_qsat_result=((zz_c-deux*zz_b+zz_a)*(zz_f-un)+zz_c-zz_b)/pres_in
280      !--------------------
281    END FUNCTION dev_qsat
282    !===
283    SUBROUTINE qsfrict_init
284      !!--------------------------------------------------------------------
285      !! The qsfrict_init routine initialises qsfrict array
286      !! to store precalculated value for qsat
287      !!--------------------------------------------------------------------
288      INTEGER(i_std)        :: ji
289      REAL(r_std)     :: zrapp,zcorr,ztemperature,zqsat
290      !---------------------------------------------------------------------
291      ! initialisation
292      zrapp = msmlr_h2o/msmlr_air
293      zcorr = 0.00320991_r_std
294      ! computes saturated humidity one time and store in qsfrict local array
295      DO ji=100,max_temp
296         ztemperature = FLOAT(ji)
297         IF (ztemperature < 273._r_std) THEN
298            zqsat = zrapp*10.0_r_std**(2.07023_r_std-zcorr*ztemperature &
299                 &             -2484.896/ztemperature+3.56654*LOG10(ztemperature))
300         ELSE
301            zqsat = zrapp*10.0**(23.8319-2948.964/ztemperature &
302                 &              -5.028*LOG10(ztemperature) &
303                 &              -29810.16*EXP(-0.0699382*ztemperature) &
304                 &              +25.21935*EXP(-2999.924/ztemperature))
305         ENDIF
306         qsfrict (ji) = zqsat
307      ENDDO
308      !-
309      qsfrict(1:100) = zero
310      !-
311      IF (long_print) WRITE (numout,*) ' qsfrict_init done'
312      !--------------------------
313    END SUBROUTINE qsfrict_init
314
315
316
317END MODULE qsat_moisture
Note: See TracBrowser for help on using the repository browser.