source: LMDZ6/branches/Amaury_dev/libf/phylmd/stdlevvar_mod.F90 @ 5157

Last change on this file since 5157 was 5144, checked in by abarral, 7 weeks ago

Put YOMCST.h into modules

  • Property svn:executable set to *
File size: 25.8 KB
Line 
1MODULE stdlevvar_mod
2
3  ! This module contains main procedures for calculation
4  ! of temperature, specific humidity and wind at a reference level
5
6  USE cdrag_mod
7  USE screenp_mod
8  USE screenc_mod
9  IMPLICIT NONE
10
11CONTAINS
12
13  !****************************************************************************************
14
15  !r original routine svn3623
16
17  SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
18          u1, v1, t1, q1, z1, &
19          ts1, qsurf, z0m, z0h, psol, pat1, &
20          t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol)
21    USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
22    USE lmdz_yoethf
23    USE lmdz_yomcst
24
25    IMPLICIT NONE
26    !-------------------------------------------------------------------------
27
28    ! Objet : calcul de la temperature et l'humidite relative a 2m et du
29    !         module du vent a 10m a partir des relations de Dyer-Businger et
30    !         des equations de Louis.
31
32    ! Reference : Hess, Colman et McAvaney (1995)
33
34    ! I. Musat, 01.07.2002
35
36    !AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
37
38    !-------------------------------------------------------------------------
39
40    ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
41    ! knon----input-I- nombre de points pour un type de surface
42    ! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
43    ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
44    ! u1------input-R- vent zonal au 1er niveau du modele
45    ! v1------input-R- vent meridien au 1er niveau du modele
46    ! t1------input-R- temperature de l'air au 1er niveau du modele
47    ! q1------input-R- humidite relative au 1er niveau du modele
48    ! z1------input-R- geopotentiel au 1er niveau du modele
49    ! ts1-----input-R- temperature de l'air a la surface
50    ! qsurf---input-R- humidite relative a la surface
51    ! z0m, z0h---input-R- rugosite
52    ! psol----input-R- pression au sol
53    ! pat1----input-R- pression au 1er niveau du modele
54
55    ! t_2m---output-R- temperature de l'air a 2m
56    ! q_2m---output-R- humidite relative a 2m
57    ! u_10m--output-R- vitesse du vent a 10m
58    !AM
59    ! t_10m--output-R- temperature de l'air a 10m
60    ! q_10m--output-R- humidite specifique a 10m
61    ! ustar--output-R- u*
62
63    INTEGER, INTENT(IN) :: klon, knon, nsrf
64    LOGICAL, INTENT(IN) :: zxli
65    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, t1, q1, z1, ts1
66    REAL, DIMENSION(klon), INTENT(IN) :: qsurf
67    REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h
68    REAL, DIMENSION(klon), INTENT(IN) :: psol, pat1
69
70    REAL, DIMENSION(klon), INTENT(OUT) :: t_2m, q_2m, ustar
71    REAL, DIMENSION(klon), INTENT(OUT) :: u_10m, t_10m, q_10m
72    REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh
73    REAL, DIMENSION(klon), INTENT(IN) :: prain
74    REAL, DIMENSION(klon), INTENT(IN) :: tsol
75    !-------------------------------------------------------------------------
76
77    ! Quelques constantes et options:
78
79    ! RKAR : constante de von Karman
80    REAL, PARAMETER :: RKAR = 0.40
81    ! niter : nombre iterations calcul "corrector"
82    !     INTEGER, parameter :: niter=6, ncon=niter-1
83    INTEGER, parameter :: niter = 2, ncon = niter - 1
84
85    ! Variables locales
86    INTEGER :: i, n
87    REAL :: zref
88    REAL, DIMENSION(klon) :: speed
89    ! tpot : temperature potentielle
90    REAL, DIMENSION(klon) :: tpot
91    REAL, DIMENSION(klon) :: zri1, cdran
92    REAL, DIMENSION(klon) :: cdram, cdrah
93    ! ri1 : nb. de Richardson entre la surface --> la 1ere couche
94    REAL, DIMENSION(klon) :: ri1
95    REAL, DIMENSION(klon) :: testar, qstar
96    REAL, DIMENSION(klon) :: zdte, zdq
97    ! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney
98    DOUBLE PRECISION, DIMENSION(klon) :: lmon
99    DOUBLE PRECISION, parameter :: eps = 1.0D-20
100    REAL, DIMENSION(klon) :: delu, delte, delq
101    REAL, DIMENSION(klon) :: u_zref, te_zref, q_zref
102    REAL, DIMENSION(klon) :: temp, pref
103    LOGICAL :: okri
104    REAL, DIMENSION(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
105    !convertgence
106    REAL, DIMENSION(klon) :: te_zref_con, q_zref_con
107    REAL, DIMENSION(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
108    REAL, DIMENSION(klon) :: ok_pred, ok_corr, zri_zero
109    !     REAL, DIMENSION(klon) :: conv_te, conv_q
110    !-------------------------------------------------------------------------
111    DO i = 1, knon
112      speed(i) = SQRT(u1(i)**2 + v1(i)**2)
113      ri1(i) = 0.0
114    ENDDO
115
116    okri = .FALSE.
117    !      CALL coefcdrag(klon, knon, nsrf, zxli, &
118    ! &                   speed, t1, q1, z1, psol, &
119    ! &                   ts1, qsurf, rugos, okri, ri1,  &
120    ! &                   cdram, cdrah, cdran, zri1, pref)
121    ! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag
122
123    CALL cdrag(knon, nsrf, &
124            speed, t1, q1, z1, &
125            psol, s_pblh, ts1, qsurf, z0m, z0h, &
126            zri_zero, 0, &
127            cdram, cdrah, zri1, pref, prain, tsol, pat1)
128
129    ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
130    IF (ok_prescr_ust) THEN
131      DO i = 1, knon
132        print *, 'cdram avant=', cdram(i)
133        cdram(i) = ust * ust / speed(i) / speed(i)
134        print *, 'cdram ust speed apres=', cdram(i), ust, speed
135      ENDDO
136    ENDIF
137
138    !---------Star variables----------------------------------------------------
139
140    DO i = 1, knon
141      ri1(i) = zri1(i)
142      tpot(i) = t1(i) * (psol(i) / pat1(i))**RKAPPA
143      ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
144      zdte(i) = tpot(i) - ts1(i)
145      zdq(i) = max(q1(i), 0.0) - max(qsurf(i), 0.0)
146
147
148      !IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
149      zdte(i) = sign(max(abs(zdte(i)), 1.e-10), zdte(i))
150
151      testar(i) = (cdrah(i) * zdte(i) * speed(i)) / ustar(i)
152      qstar(i) = (cdrah(i) * zdq(i) * speed(i)) / ustar(i)
153      lmon(i) = (ustar(i) * ustar(i) * tpot(i)) / &
154              (RKAR * RG * testar(i))
155    ENDDO
156
157    !----------First aproximation of variables at zref --------------------------
158    zref = 2.0
159    CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
160            ts1, qsurf, z0m, lmon, &
161            ustar, testar, qstar, zref, &
162            delu, delte, delq)
163
164    DO i = 1, knon
165      u_zref(i) = delu(i)
166      q_zref(i) = max(qsurf(i), 0.0) + delq(i)
167      te_zref(i) = ts1(i) + delte(i)
168      temp(i) = te_zref(i) * (psol(i) / pat1(i))**(-RKAPPA)
169      q_zref_p(i) = q_zref(i)
170      !       te_zref_p(i) = te_zref(i)
171      temp_p(i) = temp(i)
172    ENDDO
173
174    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
175
176    DO n = 1, niter
177
178      okri = .TRUE.
179      CALL screenc(klon, knon, nsrf, zxli, &
180              u_zref, temp, q_zref, zref, &
181              ts1, qsurf, z0m, z0h, psol, &
182              ustar, testar, qstar, okri, ri1, &
183              pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
184
185      DO i = 1, knon
186        u_zref(i) = delu(i)
187        q_zref(i) = delq(i) + max(qsurf(i), 0.0)
188        te_zref(i) = delte(i) + ts1(i)
189
190        ! return to normal temperature
191
192        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
193        !         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
194        !                 (1 + RVTMP2 * max(q_zref(i),0.0))
195
196        !IM +++
197        !         IF(temp(i).GT.350.) THEN
198        !           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
199        !         ENDIF
200        !IM ---
201
202        IF(n==ncon) THEN
203          te_zref_con(i) = te_zref(i)
204          q_zref_con(i) = q_zref(i)
205        ENDIF
206
207      ENDDO
208
209    ENDDO
210
211    ! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
212
213    !       DO i = 1, knon
214    !         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
215    !         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
216    !IM +++
217    !         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
218    !           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
219    !           q_zref_con(i),q_zref(i),conv_q(i)
220    !         ENDIF
221    !IM ---
222    !       ENDDO
223
224    DO i = 1, knon
225      q_zref_c(i) = q_zref(i)
226      temp_c(i) = temp(i)
227
228      !       IF(zri1(i).LT.0.) THEN
229      !         IF(nsrf.EQ.1) THEN
230      !           ok_pred(i)=1.
231      !           ok_corr(i)=0.
232      !         ELSE
233      !           ok_pred(i)=0.
234      !           ok_corr(i)=1.
235      !         ENDIF
236      !       ELSE
237      !         ok_pred(i)=0.
238      !         ok_corr(i)=1.
239      !       ENDIF
240
241      ok_pred(i) = 0.
242      ok_corr(i) = 1.
243
244      t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
245      q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
246      !IM +++
247      !       IF(n.EQ.niter) THEN
248      !       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
249      !         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i)
250      !       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
251      !         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i)
252      !       ENDIF
253      !       ENDIF
254      !IM ---
255    ENDDO
256
257
258    !----------First aproximation of variables at zref --------------------------
259
260    zref = 10.0
261    CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
262            ts1, qsurf, z0m, lmon, &
263            ustar, testar, qstar, zref, &
264            delu, delte, delq)
265
266    DO i = 1, knon
267      u_zref(i) = delu(i)
268      q_zref(i) = max(qsurf(i), 0.0) + delq(i)
269      te_zref(i) = ts1(i) + delte(i)
270      temp(i) = te_zref(i) * (psol(i) / pat1(i))**(-RKAPPA)
271      !       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
272      !                 (1 + RVTMP2 * max(q_zref(i),0.0))
273      u_zref_p(i) = u_zref(i)
274    ENDDO
275
276    ! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995
277
278    DO n = 1, niter
279
280      okri = .TRUE.
281      CALL screenc(klon, knon, nsrf, zxli, &
282              u_zref, temp, q_zref, zref, &
283              ts1, qsurf, z0m, z0h, psol, &
284              ustar, testar, qstar, okri, ri1, &
285              pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
286
287      DO i = 1, knon
288        u_zref(i) = delu(i)
289        q_zref(i) = delq(i) + max(qsurf(i), 0.0)
290        te_zref(i) = delte(i) + ts1(i)
291        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
292        !         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
293        !                   (1 + RVTMP2 * max(q_zref(i),0.0))
294      ENDDO
295
296    ENDDO
297
298    DO i = 1, knon
299      u_zref_c(i) = u_zref(i)
300
301      u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
302
303      !AM
304      q_zref_c(i) = q_zref(i)
305      temp_c(i) = temp(i)
306      t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
307      q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
308      !MA
309    ENDDO
310
311  END SUBROUTINE stdlevvar
312
313  SUBROUTINE stdlevvarn(klon, knon, nsrf, zxli, &
314          u1, v1, t1, q1, z1, &
315          ts1, qsurf, z0m, z0h, psol, pat1, &
316          t_2m, q_2m, t_10m, q_10m, u_10m, ustar, &
317          n2mout)
318
319    USE lmdz_ioipsl_getin_p, ONLY: getin_p
320    USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf
321    USE lmdz_yoethf
322    USE lmdz_yomcst
323
324    IMPLICIT NONE
325    !-------------------------------------------------------------------------
326
327    ! Objet : calcul de la temperature et l'humidite relative a 2m et du
328    !         module du vent a 10m a partir des relations de Dyer-Businger et
329    !         des equations de Louis.
330
331    ! Reference : Hess, Colman et McAvaney (1995)
332
333    ! I. Musat, 01.07.2002
334
335    !AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
336
337    !-------------------------------------------------------------------------
338
339    ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
340    ! knon----input-I- nombre de points pour un type de surface
341    ! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
342    ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
343    ! u1------input-R- vent zonal au 1er niveau du modele
344    ! v1------input-R- vent meridien au 1er niveau du modele
345    ! t1------input-R- temperature de l'air au 1er niveau du modele
346    ! q1------input-R- humidite relative au 1er niveau du modele
347    ! z1------input-R- geopotentiel au 1er niveau du modele
348    ! ts1-----input-R- temperature de l'air a la surface
349    ! qsurf---input-R- humidite relative a la surface
350    ! z0m, z0h---input-R- rugosite
351    ! psol----input-R- pression au sol
352    ! pat1----input-R- pression au 1er niveau du modele
353
354    ! t_2m---output-R- temperature de l'air a 2m
355    ! q_2m---output-R- humidite relative a 2m
356    ! u_2m--output-R- vitesse du vent a 2m
357    ! u_10m--output-R- vitesse du vent a 10m
358    ! ustar--output-R- u*
359    !AM
360    ! t_10m--output-R- temperature de l'air a 10m
361    ! q_10m--output-R- humidite specifique a 10m
362
363    INTEGER, INTENT(IN) :: klon, knon, nsrf
364    LOGICAL, INTENT(IN) :: zxli
365    REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, t1, q1, ts1, z1
366    REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h
367    REAL, DIMENSION(klon), INTENT(IN) :: qsurf
368    REAL, DIMENSION(klon), INTENT(IN) :: psol, pat1
369
370    REAL, DIMENSION(klon), INTENT(OUT) :: t_2m, q_2m, ustar
371    REAL, DIMENSION(klon), INTENT(OUT) :: u_10m, t_10m, q_10m
372    INTEGER, DIMENSION(klon, 6), INTENT(OUT) :: n2mout
373
374    REAL, DIMENSION(klon) :: u_2m
375    REAL, DIMENSION(klon) :: cdrm2m, cdrh2m, ri2m
376    REAL, DIMENSION(klon) :: cdram, cdrah, zri1
377    REAL, DIMENSION(klon) :: cdmn1, cdhn1, fm1, fh1
378    REAL, DIMENSION(klon) :: cdmn2m, cdhn2m, fm2m, fh2m
379    REAL, DIMENSION(klon) :: ri2m_new
380    REAL, DIMENSION(klon) :: s_pblh
381    REAL, DIMENSION(klon) :: prain
382    REAL, DIMENSION(klon) :: tsol
383    !-------------------------------------------------------------------------
384
385    ! Quelques constantes et options:
386
387    ! RKAR : constante de von Karman
388    REAL, PARAMETER :: RKAR = 0.40
389    ! niter : nombre iterations calcul "corrector"
390    !     INTEGER, parameter :: niter=6, ncon=niter-1
391    !IM 071020     INTEGER, parameter :: niter=2, ncon=niter-1
392    INTEGER, parameter :: niter = 2, ncon = niter
393    !     INTEGER, parameter :: niter=6, ncon=niter
394
395    ! Variables locales
396    INTEGER :: i, n
397    REAL :: zref
398    REAL, DIMENSION(klon) :: speed
399    ! tpot : temperature potentielle
400    REAL, DIMENSION(klon) :: tpot
401    REAL, DIMENSION(klon) :: cdran
402    ! ri1 : nb. de Richardson entre la surface --> la 1ere couche
403    REAL, DIMENSION(klon) :: ri1
404    DOUBLE PRECISION, parameter :: eps = 1.0D-20
405    REAL, DIMENSION(klon) :: delu, delte, delq
406    REAL, DIMENSION(klon) :: delh, delm
407    REAL, DIMENSION(klon) :: delh_new, delm_new
408    REAL, DIMENSION(klon) :: u_zref, te_zref, q_zref
409    REAL, DIMENSION(klon) :: u_zref_pnew, te_zref_pnew, q_zref_pnew
410    REAL, DIMENSION(klon) :: temp, pref
411    REAL, DIMENSION(klon) :: temp_new, pref_new
412    LOGICAL :: okri
413    REAL, DIMENSION(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
414    REAL, DIMENSION(klon) :: u_zref_p_new, te_zref_p_new, temp_p_new, q_zref_p_new
415    !convergence
416    REAL, DIMENSION(klon) :: te_zref_con, q_zref_con
417    REAL, DIMENSION(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
418    REAL, DIMENSION(klon) :: ok_pred, ok_corr
419
420    REAL, DIMENSION(klon) :: cdrm10m, cdrh10m, ri10m
421    REAL, DIMENSION(klon) :: cdmn10m, cdhn10m, fm10m, fh10m
422    REAL, DIMENSION(klon) :: cdn2m, cdn1, zri_zero
423    REAL :: CEPDUE, zdu2
424    INTEGER :: nzref, nz1
425    LOGICAL, DIMENSION(klon) :: ok_t2m_toosmall, ok_t2m_toobig
426    LOGICAL, DIMENSION(klon) :: ok_q2m_toosmall, ok_q2m_toobig
427    LOGICAL, DIMENSION(klon) :: ok_u2m_toobig
428    LOGICAL, DIMENSION(klon) :: ok_t10m_toosmall, ok_t10m_toobig
429    LOGICAL, DIMENSION(klon) :: ok_q10m_toosmall, ok_q10m_toobig
430    LOGICAL, DIMENSION(klon) :: ok_u10m_toobig
431    INTEGER, DIMENSION(klon, 6) :: n10mout
432
433    !-------------------------------------------------------------------------
434    CEPDUE = 0.1
435
436    ! n2mout : compteur des pas de temps ou t2m,q2m ou u2m sont en dehors des intervalles
437    !          [tsurf, temp], [qsurf, q1] ou [0, speed]
438    ! n10mout : compteur des pas de temps ou t10m,q10m ou u10m sont en dehors des intervalles
439    !          [tsurf, temp], [qsurf, q1] ou [0, speed]
440
441    n2mout(:, :) = 0
442    n10mout(:, :) = 0
443
444    DO i = 1, knon
445      speed(i) = MAX(SQRT(u1(i)**2 + v1(i)**2), CEPDUE)
446      ri1(i) = 0.0
447    ENDDO
448
449    okri = .FALSE.
450    CALL cdrag(knon, nsrf, &
451            speed, t1, q1, z1, &
452            psol, s_pblh, ts1, qsurf, z0m, z0h, &
453            zri_zero, 0, &
454            cdram, cdrah, zri1, pref, prain, tsol, pat1)
455
456    DO i = 1, knon
457      ri1(i) = zri1(i)
458      tpot(i) = t1(i) * (psol(i) / pat1(i))**RKAPPA
459      zdu2 = MAX(CEPDUE * CEPDUE, speed(i)**2)
460      ustar(i) = sqrt(cdram(i) * zdu2)
461
462    ENDDO
463
464    !----------First aproximation of variables at zref --------------------------
465    zref = 2.0
466
467    ! calcul first-guess en utilisant dans les calculs à 2m
468    ! le Richardson de la premiere couche atmospherique
469
470    CALL screencn(klon, knon, nsrf, zxli, &
471            speed, tpot, q1, zref, &
472            ts1, qsurf, z0m, z0h, psol, &
473            cdram, cdrah, okri, &
474            ri1, 1, &
475            pref_new, delm_new, delh_new, ri2m, &
476            s_pblh, prain, tsol, pat1)
477
478    DO i = 1, knon
479      u_zref(i) = delm_new(i) * speed(i)
480      u_zref_p(i) = u_zref(i)
481      q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
482              max(qsurf(i), 0.0) * (1 - delh_new(i))
483      q_zref_p(i) = q_zref(i)
484      te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
485      te_zref_p(i) = te_zref(i)
486
487      ! return to normal temperature
488      temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
489      temp_p(i) = temp(i)
490
491      ! compteurs ici
492
493      ok_t2m_toosmall(i) = te_zref(i)<tpot(i).AND. &
494              te_zref(i)<ts1(i)
495      ok_t2m_toobig(i) = te_zref(i)>tpot(i).AND. &
496              te_zref(i)>ts1(i)
497      ok_q2m_toosmall(i) = q_zref(i)<q1(i).AND. &
498              q_zref(i)<qsurf(i)
499      ok_q2m_toobig(i) = q_zref(i)>q1(i).AND. &
500              q_zref(i)>qsurf(i)
501      ok_u2m_toobig(i) = u_zref(i)>speed(i)
502
503      IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
504        n2mout(i, 1) = n2mout(i, 1) + 1
505      ENDIF
506      IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
507        n2mout(i, 3) = n2mout(i, 3) + 1
508      ENDIF
509      IF(ok_u2m_toobig(i)) THEN
510        n2mout(i, 5) = n2mout(i, 5) + 1
511      ENDIF
512
513      IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
514              ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
515              ok_u2m_toobig(i)) THEN
516        delm_new(i) = min(max(delm_new(i), 0.), 1.)
517        delh_new(i) = min(max(delh_new(i), 0.), 1.)
518        u_zref(i) = delm_new(i) * speed(i)
519        u_zref_p(i) = u_zref(i)
520        q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
521                max(qsurf(i), 0.0) * (1 - delh_new(i))
522        q_zref_p(i) = q_zref(i)
523        te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
524        te_zref_p(i) = te_zref(i)
525
526        ! return to normal temperature
527        temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
528        temp_p(i) = temp(i)
529      ENDIF
530
531    ENDDO
532
533    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
534
535    DO n = 1, niter
536
537      okri = .TRUE.
538      CALL screencn(klon, knon, nsrf, zxli, &
539              u_zref, temp, q_zref, zref, &
540              ts1, qsurf, z0m, z0h, psol, &
541              cdram, cdrah, okri, &
542              ri1, 0, &
543              pref, delm, delh, ri2m, &
544              s_pblh, prain, tsol, pat1)
545
546      DO i = 1, knon
547        u_zref(i) = delm(i) * speed(i)
548        q_zref(i) = delh(i) * max(q1(i), 0.0) + &
549                max(qsurf(i), 0.0) * (1 - delh(i))
550        te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
551
552        ! return to normal temperature
553        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
554
555        ! compteurs ici
556
557        ok_t2m_toosmall(i) = te_zref(i)<tpot(i).AND. &
558                te_zref(i)<ts1(i)
559        ok_t2m_toobig(i) = te_zref(i)>tpot(i).AND. &
560                te_zref(i)>ts1(i)
561        ok_q2m_toosmall(i) = q_zref(i)<q1(i).AND. &
562                q_zref(i)<qsurf(i)
563        ok_q2m_toobig(i) = q_zref(i)>q1(i).AND. &
564                q_zref(i)>qsurf(i)
565        ok_u2m_toobig(i) = u_zref(i)>speed(i)
566
567        IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
568          n2mout(i, 2) = n2mout(i, 2) + 1
569        ENDIF
570        IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
571          n2mout(i, 4) = n2mout(i, 4) + 1
572        ENDIF
573        IF(ok_u2m_toobig(i)) THEN
574          n2mout(i, 6) = n2mout(i, 6) + 1
575        ENDIF
576
577        IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
578                ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
579                ok_u2m_toobig(i)) THEN
580          delm(i) = min(max(delm(i), 0.), 1.)
581          delh(i) = min(max(delh(i), 0.), 1.)
582          u_zref(i) = delm(i) * speed(i)
583          q_zref(i) = delh(i) * max(q1(i), 0.0) + &
584                  max(qsurf(i), 0.0) * (1 - delh(i))
585          te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
586          temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
587        ENDIF
588
589        IF(n==ncon) THEN
590          te_zref_con(i) = te_zref(i)
591          q_zref_con(i) = q_zref(i)
592        ENDIF
593
594      ENDDO
595
596    ENDDO
597
598    DO i = 1, knon
599      q_zref_c(i) = q_zref(i)
600      temp_c(i) = temp(i)
601
602      ok_pred(i) = 0.
603      ok_corr(i) = 1.
604
605      t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
606      q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
607
608      u_zref_c(i) = u_zref(i)
609      u_2m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
610    ENDDO
611
612
613    !----------First aproximation of variables at zref --------------------------
614
615    zref = 10.0
616
617    CALL screencn(klon, knon, nsrf, zxli, &
618            speed, tpot, q1, zref, &
619            ts1, qsurf, z0m, z0h, psol, &
620            cdram, cdrah, okri, &
621            ri1, 1, &
622            pref_new, delm_new, delh_new, ri10m, &
623            s_pblh, prain, tsol, pat1)
624
625    DO i = 1, knon
626      u_zref(i) = delm_new(i) * speed(i)
627      q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
628              max(qsurf(i), 0.0) * (1 - delh_new(i))
629      te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
630      temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
631      u_zref_p(i) = u_zref(i)
632
633      ! compteurs ici
634
635      ok_t10m_toosmall(i) = te_zref(i)<tpot(i).AND. &
636              te_zref(i)<ts1(i)
637      ok_t10m_toobig(i) = te_zref(i)>tpot(i).AND. &
638              te_zref(i)>ts1(i)
639      ok_q10m_toosmall(i) = q_zref(i)<q1(i).AND. &
640              q_zref(i)<qsurf(i)
641      ok_q10m_toobig(i) = q_zref(i)>q1(i).AND. &
642              q_zref(i)>qsurf(i)
643      ok_u10m_toobig(i) = u_zref(i)>speed(i)
644
645      IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
646        n10mout(i, 1) = n10mout(i, 1) + 1
647      ENDIF
648      IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
649        n10mout(i, 3) = n10mout(i, 3) + 1
650      ENDIF
651      IF(ok_u10m_toobig(i)) THEN
652        n10mout(i, 5) = n10mout(i, 5) + 1
653      ENDIF
654
655      IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
656              ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
657              ok_u10m_toobig(i)) THEN
658        delm_new(i) = min(max(delm_new(i), 0.), 1.)
659        delh_new(i) = min(max(delh_new(i), 0.), 1.)
660        u_zref(i) = delm_new(i) * speed(i)
661        u_zref_p(i) = u_zref(i)
662        q_zref(i) = delh_new(i) * max(q1(i), 0.0) + &
663                max(qsurf(i), 0.0) * (1 - delh_new(i))
664        te_zref(i) = delh_new(i) * tpot(i) + ts1(i) * (1 - delh_new(i))
665        temp(i) = te_zref(i) * (psol(i) / pref_new(i))**(-RKAPPA)
666      ENDIF
667
668    ENDDO
669
670    ! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
671
672    DO n = 1, niter
673
674      okri = .TRUE.
675      CALL screencn(klon, knon, nsrf, zxli, &
676              u_zref, temp, q_zref, zref, &
677              ts1, qsurf, z0m, z0h, psol, &
678              cdram, cdrah, okri, &
679              ri1, 0, &
680              pref, delm, delh, ri10m, &
681              s_pblh, prain, tsol, pat1)
682
683      DO i = 1, knon
684        u_zref(i) = delm(i) * speed(i)
685        q_zref(i) = delh(i) * max(q1(i), 0.0) + &
686                max(qsurf(i), 0.0) * (1 - delh(i))
687        te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
688
689        ! return to normal temperature
690        temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
691
692        ! compteurs ici
693
694        ok_t10m_toosmall(i) = te_zref(i)<tpot(i).AND. &
695                te_zref(i)<ts1(i)
696        ok_t10m_toobig(i) = te_zref(i)>tpot(i).AND. &
697                te_zref(i)>ts1(i)
698        ok_q10m_toosmall(i) = q_zref(i)<q1(i).AND. &
699                q_zref(i)<qsurf(i)
700        ok_q10m_toobig(i) = q_zref(i)>q1(i).AND. &
701                q_zref(i)>qsurf(i)
702        ok_u10m_toobig(i) = u_zref(i)>speed(i)
703
704        IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
705          n10mout(i, 2) = n10mout(i, 2) + 1
706        ENDIF
707        IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
708          n10mout(i, 4) = n10mout(i, 4) + 1
709        ENDIF
710        IF(ok_u10m_toobig(i)) THEN
711          n10mout(i, 6) = n10mout(i, 6) + 1
712        ENDIF
713
714        IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
715                ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
716                ok_u10m_toobig(i)) THEN
717          delm(i) = min(max(delm(i), 0.), 1.)
718          delh(i) = min(max(delh(i), 0.), 1.)
719          u_zref(i) = delm(i) * speed(i)
720          q_zref(i) = delh(i) * max(q1(i), 0.0) + &
721                  max(qsurf(i), 0.0) * (1 - delh(i))
722          te_zref(i) = delh(i) * tpot(i) + ts1(i) * (1 - delh(i))
723          temp(i) = te_zref(i) * (psol(i) / pref(i))**(-RKAPPA)
724        ENDIF
725
726        IF(n==ncon) THEN
727          te_zref_con(i) = te_zref(i)
728          q_zref_con(i) = q_zref(i)
729        ENDIF
730
731      ENDDO
732
733    ENDDO
734
735    DO i = 1, knon
736      q_zref_c(i) = q_zref(i)
737      temp_c(i) = temp(i)
738
739      ok_pred(i) = 0.
740      ok_corr(i) = 1.
741
742      t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
743      q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
744
745      u_zref_c(i) = u_zref(i)
746      u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
747    ENDDO
748
749  END SUBROUTINE stdlevvarn
750
751END MODULE stdlevvar_mod
Note: See TracBrowser for help on using the repository browser.