source: LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90 @ 4257

Last change on this file since 4257 was 3839, checked in by musat, 4 years ago

Ajout ustar pour diagnostics pbl

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