source: LMDZ6/branches/contrails/libf/phylmd/stdlevvar_mod.f90 @ 5445

Last change on this file since 5445 was 5301, checked in by abarral, 8 weeks ago

Turn tsoilnudge.h fcg_gcssold.h flux_arp.h into module

  • Property svn:executable set to *
File size: 26.1 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 yomcst_mod_h
8      USE cdrag_mod
9  USE screenp_mod
10  USE screenc_mod
11  IMPLICIT NONE
12
13CONTAINS
14!
15!****************************************************************************************
16!
17!r original routine svn3623
18!
19      SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
20                           u1, v1, t1, q1, z1, &
21                           ts1, qsurf, z0m, z0h, psol, pat1, &
22                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol)
23        USE yoethf_mod_h
24        USE flux_arp_mod_h
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.EQ.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      RETURN
312      END subroutine stdlevvar
313!
314      SUBROUTINE stdlevvarn(klon, knon, nsrf, zxli, &
315                           u1, v1, t1, q1, z1, &
316                           ts1, qsurf, z0m, z0h, psol, pat1, &
317                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar, &
318                           n2mout)
319!
320      USE yomcst_mod_h
321      USE ioipsl_getin_p_mod, ONLY : getin_p
322      USE yoethf_mod_h
323      USE flux_arp_mod_h
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!
457      DO i = 1, knon
458        ri1(i) = zri1(i)
459        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
460        zdu2 = MAX(CEPDUE*CEPDUE, speed(i)**2)
461        ustar(i) = sqrt(cdram(i) * zdu2)
462!
463      ENDDO
464!
465!----------First aproximation of variables at zref --------------------------
466      zref = 2.0
467!
468! calcul first-guess en utilisant dans les calculs à 2m
469! le Richardson de la premiere couche atmospherique
470!
471       CALL screencn(klon, knon, nsrf, zxli, &
472 &                   speed, tpot, q1, zref, &
473 &                   ts1, qsurf, z0m, z0h, psol, &           
474 &                   cdram, cdrah,  okri, &
475 &                   ri1, 1, &
476 &                   pref_new, delm_new, delh_new, ri2m, &
477 &                   s_pblh, prain, tsol, pat1      )
478!
479       DO i = 1, knon
480         u_zref(i) = delm_new(i)*speed(i)
481         u_zref_p(i) = u_zref(i)
482         q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
483         &           max(qsurf(i),0.0)*(1-delh_new(i))
484         q_zref_p(i) = q_zref(i)
485         te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
486         te_zref_p(i) = te_zref(i)
487!
488! return to normal temperature
489         temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
490         temp_p(i) = temp(i)
491!
492! compteurs ici
493!
494         ok_t2m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
495         & te_zref(i).LT.ts1(i)
496         ok_t2m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
497         & te_zref(i).GT.ts1(i)
498         ok_q2m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
499         & q_zref(i).LT.qsurf(i)
500         ok_q2m_toobig(i)=q_zref(i).GT.q1(i).AND. &
501         & q_zref(i).GT.qsurf(i)
502         ok_u2m_toobig(i)=u_zref(i).GT.speed(i)
503!
504         IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
505             n2mout(i,1)=n2mout(i,1)+1
506         ENDIF
507         IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
508             n2mout(i,3)=n2mout(i,3)+1
509         ENDIF
510         IF(ok_u2m_toobig(i)) THEN
511             n2mout(i,5)=n2mout(i,5)+1
512         ENDIF
513!
514         IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
515          & ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
516          & ok_u2m_toobig(i)) THEN
517             delm_new(i)=min(max(delm_new(i),0.),1.)
518             delh_new(i)=min(max(delh_new(i),0.),1.)
519             u_zref(i) = delm_new(i)*speed(i)
520             u_zref_p(i) = u_zref(i)
521             q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
522         &               max(qsurf(i),0.0)*(1-delh_new(i))
523             q_zref_p(i) = q_zref(i)
524             te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
525             te_zref_p(i) = te_zref(i)
526!
527! return to normal temperature
528             temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
529             temp_p(i) = temp(i)
530         ENDIF
531!
532       ENDDO
533!
534! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
535!
536      DO n = 1, niter
537!
538        okri=.TRUE.
539        CALL screencn(klon, knon, nsrf, zxli, &
540 &                   u_zref, temp, q_zref, zref, &
541 &                   ts1, qsurf, z0m, z0h, psol, &
542 &                   cdram, cdrah,  okri, &
543 &                   ri1, 0, &
544 &                   pref, delm, delh, ri2m, &
545 &                   s_pblh, prain, tsol, pat1      )
546!
547        DO i = 1, knon
548          u_zref(i) = delm(i)*speed(i)
549          q_zref(i) = delh(i)*max(q1(i),0.0) + &
550          &           max(qsurf(i),0.0)*(1-delh(i))
551          te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
552!
553! return to normal temperature
554          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
555!
556! compteurs ici
557!
558          ok_t2m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
559          & te_zref(i).LT.ts1(i)
560          ok_t2m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
561          & te_zref(i).GT.ts1(i)
562          ok_q2m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
563          & q_zref(i).LT.qsurf(i)
564          ok_q2m_toobig(i)=q_zref(i).GT.q1(i).AND. &
565          & q_zref(i).GT.qsurf(i)
566          ok_u2m_toobig(i)=u_zref(i).GT.speed(i)
567!
568          IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
569              n2mout(i,2)=n2mout(i,2)+1
570          ENDIF
571          IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
572              n2mout(i,4)=n2mout(i,4)+1
573          ENDIF
574          IF(ok_u2m_toobig(i)) THEN
575              n2mout(i,6)=n2mout(i,6)+1
576          ENDIF
577!
578          IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
579           & ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
580           & ok_u2m_toobig(i)) THEN
581              delm(i)=min(max(delm(i),0.),1.)
582              delh(i)=min(max(delh(i),0.),1.)
583              u_zref(i) = delm(i)*speed(i)
584              q_zref(i) = delh(i)*max(q1(i),0.0) + &
585          &           max(qsurf(i),0.0)*(1-delh(i))
586              te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
587              temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
588          ENDIF
589!
590!
591          IF(n.EQ.ncon) THEN
592            te_zref_con(i) = te_zref(i)
593            q_zref_con(i) = q_zref(i)
594          ENDIF
595!
596        ENDDO
597!
598      ENDDO
599!
600      DO i = 1, knon
601        q_zref_c(i) = q_zref(i)
602        temp_c(i) = temp(i)
603!
604        ok_pred(i)=0.
605        ok_corr(i)=1.
606!
607        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
608        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
609!
610        u_zref_c(i) = u_zref(i)
611        u_2m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
612      ENDDO
613!
614!
615!----------First aproximation of variables at zref --------------------------
616!
617      zref = 10.0
618!
619       CALL screencn(klon, knon, nsrf, zxli, &
620 &                   speed, tpot, q1, zref, &
621 &                   ts1, qsurf, z0m, z0h, psol, &           
622 &                   cdram, cdrah,  okri, &
623 &                   ri1, 1, &
624 &                   pref_new, delm_new, delh_new, ri10m, &
625 &                   s_pblh, prain, tsol, pat1      )
626!
627       DO i = 1, knon
628         u_zref(i) = delm_new(i)*speed(i)
629         q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
630         &           max(qsurf(i),0.0)*(1-delh_new(i))
631         te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
632         temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
633         u_zref_p(i) = u_zref(i)
634!
635! compteurs ici
636!
637         ok_t10m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
638         & te_zref(i).LT.ts1(i)
639         ok_t10m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
640         & te_zref(i).GT.ts1(i)
641         ok_q10m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
642         & q_zref(i).LT.qsurf(i)
643         ok_q10m_toobig(i)=q_zref(i).GT.q1(i).AND. &
644         & q_zref(i).GT.qsurf(i)
645         ok_u10m_toobig(i)=u_zref(i).GT.speed(i)
646!
647         IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
648             n10mout(i,1)=n10mout(i,1)+1
649         ENDIF
650         IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
651             n10mout(i,3)=n10mout(i,3)+1
652         ENDIF
653         IF(ok_u10m_toobig(i)) THEN
654             n10mout(i,5)=n10mout(i,5)+1
655         ENDIF
656!
657         IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
658          & ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
659          & ok_u10m_toobig(i)) THEN
660             delm_new(i)=min(max(delm_new(i),0.),1.)
661             delh_new(i)=min(max(delh_new(i),0.),1.)
662             u_zref(i) = delm_new(i)*speed(i)
663             u_zref_p(i) = u_zref(i)
664             q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
665         &               max(qsurf(i),0.0)*(1-delh_new(i))
666             te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
667             temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
668         ENDIF
669!
670       ENDDO
671!
672! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
673!
674      DO n = 1, niter
675!
676        okri=.TRUE.
677        CALL screencn(klon, knon, nsrf, zxli, &
678 &                   u_zref, temp, q_zref, zref, &
679 &                   ts1, qsurf, z0m, z0h, psol, &
680 &                   cdram, cdrah,  okri, &
681 &                   ri1, 0, &
682 &                   pref, delm, delh, ri10m, &
683 &                   s_pblh, prain, tsol, pat1      )
684!
685        DO i = 1, knon
686          u_zref(i) = delm(i)*speed(i)
687          q_zref(i) = delh(i)*max(q1(i),0.0) + &
688          &           max(qsurf(i),0.0)*(1-delh(i))
689          te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
690!
691! return to normal temperature
692          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
693!
694! compteurs ici
695!
696          ok_t10m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
697          & te_zref(i).LT.ts1(i)
698          ok_t10m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
699          & te_zref(i).GT.ts1(i)
700          ok_q10m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
701          & q_zref(i).LT.qsurf(i)
702          ok_q10m_toobig(i)=q_zref(i).GT.q1(i).AND. &
703          & q_zref(i).GT.qsurf(i)
704          ok_u10m_toobig(i)=u_zref(i).GT.speed(i)
705!
706          IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
707              n10mout(i,2)=n10mout(i,2)+1
708          ENDIF
709          IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
710              n10mout(i,4)=n10mout(i,4)+1
711          ENDIF
712          IF(ok_u10m_toobig(i)) THEN
713              n10mout(i,6)=n10mout(i,6)+1
714          ENDIF
715!
716          IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
717           & ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
718           & ok_u10m_toobig(i)) THEN
719              delm(i)=min(max(delm(i),0.),1.)
720              delh(i)=min(max(delh(i),0.),1.)
721              u_zref(i) = delm(i)*speed(i)
722              q_zref(i) = delh(i)*max(q1(i),0.0) + &
723          &           max(qsurf(i),0.0)*(1-delh(i))
724              te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
725              temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
726          ENDIF
727!
728!
729          IF(n.EQ.ncon) THEN
730            te_zref_con(i) = te_zref(i)
731            q_zref_con(i) = q_zref(i)
732          ENDIF
733!
734        ENDDO
735!
736      ENDDO
737!
738      DO i = 1, knon
739        q_zref_c(i) = q_zref(i)
740        temp_c(i) = temp(i)
741!
742        ok_pred(i)=0.
743        ok_corr(i)=1.
744!
745        t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
746        q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
747!
748        u_zref_c(i) = u_zref(i)
749        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
750      ENDDO
751!
752      RETURN
753      END subroutine stdlevvarn
754
755END MODULE stdlevvar_mod
Note: See TracBrowser for help on using the repository browser.