source: LMDZ6/branches/blowing_snow/libf/phylmd/stdlevvar_mod.F90 @ 5453

Last change on this file since 5453 was 4478, checked in by evignon, 21 months ago

commission sur le melange turbulent et les cdrag suite au dernier atelier tke

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