source: LMDZ.3.3/branches/rel-LF/libf/phylmd/stdlevvar.F90 @ 5456

Last change on this file since 5456 was 537, checked in by lmdzadmin, 21 years ago

Petit bug sur le signe de la valeur seuil du gradient tair-tsol
IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.9 KB
RevLine 
[416]1      SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
[509]2 &                         u1, v1, t1, q1, z1, &
3 &                         ts1, qsurf, rugos, psol, pat1, &
4 &                         t_2m, q_2m, u_10m)
[416]5      IMPLICIT NONE
6!-------------------------------------------------------------------------
7!
8! Objet : calcul de la temperature et l'humidite relative a 2m et du
9!         module du vent a 10m a partir des relations de Dyer-Businger et
10!         des equations de Louis.
11!
12! Reference : Hess, Colman et McAvaney (1995)       
13!
14! I. Musat, 01.07.2002
15!-------------------------------------------------------------------------
16!
17! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
18! knon----input-I- nombre de points pour un type de surface
19! nsrf----input-I- indice pour le type de surface; voir indicesol.inc
20! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
21! u1------input-R- vent zonal au 1er niveau du modele
22! v1------input-R- vent meridien au 1er niveau du modele
23! t1------input-R- temperature de l'air au 1er niveau du modele
24! q1------input-R- humidite relative au 1er niveau du modele
25! z1------input-R- geopotentiel au 1er niveau du modele
26! ts1-----input-R- temperature de l'air a la surface
27! qsurf---input-R- humidite relative a la surface
28! rugos---input-R- rugosite
29! psol----input-R- pression au sol
30! pat1----input-R- pression au 1er niveau du modele
31!
32! t_2m---output-R- temperature de l'air a 2m
33! q_2m---output-R- humidite relative a 2m
34! u_10m--output-R- vitesse du vent a 10m
35!
36      INTEGER, intent(in) :: klon, knon, nsrf
37      LOGICAL, intent(in) :: zxli
38      REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
39      REAL, dimension(klon), intent(in) :: qsurf, rugos
40      REAL, dimension(klon), intent(in) :: psol, pat1
41!
42      REAL, dimension(klon), intent(out) :: t_2m, q_2m, u_10m
43!-------------------------------------------------------------------------
44#include "YOMCST.inc"
45!IM PLUS
46#include "YOETHF.inc"
47!
48! Quelques constantes et options:
49!
50! RKAR : constante de von Karman
51      REAL, PARAMETER :: RKAR=0.40
52! niter : nombre iterations calcul "corrector"
[445]53!     INTEGER, parameter :: niter=6, ncon=niter-1
54      INTEGER, parameter :: niter=2, ncon=niter-1
[416]55!
56! Variables locales
57      INTEGER :: i, n
58      REAL :: zref
59      REAL, dimension(klon) :: speed
60! tpot : temperature potentielle
61      REAL, dimension(klon) :: tpot
62      REAL, dimension(klon) :: zri1, cdran
63      REAL, dimension(klon) :: cdram, cdrah
64! ri1 : nb. de Richardson entre la surface --> la 1ere couche
65      REAL, dimension(klon) :: ri1
66      REAL, dimension(klon) :: ustar, testar, qstar
67      REAL, dimension(klon) :: zdte, zdq   
68! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney
69      DOUBLE PRECISION, dimension(klon) :: lmon
70      DOUBLE PRECISION, parameter :: eps=1.0D-20
71      REAL, dimension(klon) :: delu, delte, delq
72      REAL, dimension(klon) :: u_zref, te_zref, q_zref 
73      REAL, dimension(klon) :: temp, pref
74      LOGICAL :: okri
75      REAL, dimension(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
[445]76!convertgence
77      REAL, dimension(klon) :: te_zref_con, q_zref_con
[416]78      REAL, dimension(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
79      REAL, dimension(klon) :: ok_pred, ok_corr
[472]80!     REAL, dimension(klon) :: conv_te, conv_q
[416]81!-------------------------------------------------------------------------
82      DO i=1, knon
83       speed(i)=SQRT(u1(i)**2+v1(i)**2)
84       ri1(i) = 0.0
85      ENDDO
86!
87      okri=.FALSE.
88      CALL coefcdrag(klon, knon, nsrf, zxli, &
[509]89 &                   speed, t1, q1, z1, psol, &
90 &                   ts1, qsurf, rugos, okri, ri1,  &         
91 &                   cdram, cdrah, cdran, zri1, pref)           
[416]92!
93!---------Star variables----------------------------------------------------
94!
95      DO i = 1, knon
96        ri1(i) = zri1(i)
97        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
98        ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
99        zdte(i) = tpot(i) - ts1(i)
[537]100!       PRINT*,'AVANT i,zdte',i,zdte(i)
[509]101!IM cf FH : on prend le max : pour eviter le plantage sur SUN
[537]102!IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
103        zdte(i) = sign(max(abs(zdte(i)),1.e-10),zdte(i))
104!       PRINT*,'APRES i,zdte',i,zdte(i)
105!
[416]106        zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0)
107!
108        testar(i) = (cdrah(i) * zdte(i) * speed(i))/ustar(i)
109        qstar(i) = (cdrah(i) * zdq(i) * speed(i))/ustar(i)
110        lmon(i) = (ustar(i) * ustar(i) * tpot(i))/ &
[509]111 &                (RKAR * RG * testar(i))
[416]112      ENDDO
113!
114!----------First aproximation of variables at zref --------------------------
115      zref = 2.0
116      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
[509]117 &                 ts1, qsurf, rugos, lmon, &
118 &                 ustar, testar, qstar, zref, &
119 &                 delu, delte, delq)
[416]120!
121      DO i = 1, knon
122        u_zref(i) = delu(i)
123        q_zref(i) = max(qsurf(i),0.0) + delq(i)
124        te_zref(i) = ts1(i) + delte(i)
125        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
[445]126        q_zref_p(i) = q_zref(i)
[416]127!       te_zref_p(i) = te_zref(i)
128        temp_p(i) = temp(i)
129      ENDDO
130!
131! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
132!
133      DO n = 1, niter
134!
135        okri=.TRUE.
136        CALL screenc(klon, knon, nsrf, zxli, &
[509]137 &                   u_zref, temp, q_zref, zref, &
138 &                   ts1, qsurf, rugos, psol, &           
139 &                   ustar, testar, qstar, okri, ri1, &
140 &                   pref, delu, delte, delq)
[416]141!
142        DO i = 1, knon
143          u_zref(i) = delu(i)
144          q_zref(i) = delq(i) + max(qsurf(i),0.0)
145          te_zref(i) = delte(i) + ts1(i)
146!
147! return to normal temperature
148!
149          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
150!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
151!                 (1 + RVTMP2 * max(q_zref(i),0.0))
152!
[445]153!IM +++
154!         IF(temp(i).GT.350.) THEN
155!           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
156!         ENDIF
157!IM ---
[416]158!
159        IF(n.EQ.ncon) THEN
[445]160          te_zref_con(i) = te_zref(i)
161          q_zref_con(i) = q_zref(i)
[416]162        ENDIF
163!
164        ENDDO
165!
166      ENDDO
167!
168! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
169!
[472]170!       DO i = 1, knon
171!         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
172!         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
[445]173!IM +++
[433]174!         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
[445]175!           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
176!           q_zref_con(i),q_zref(i),conv_q(i)
[433]177!         ENDIF
[445]178!IM ---
[472]179!       ENDDO
[416]180!
181      DO i = 1, knon
182        q_zref_c(i) = q_zref(i)
183        temp_c(i) = temp(i)
184!
[445]185!       IF(zri1(i).LT.0.) THEN
186!         IF(nsrf.EQ.1) THEN
187!           ok_pred(i)=1.
188!           ok_corr(i)=0.
189!         ELSE
190!           ok_pred(i)=0.
191!           ok_corr(i)=1.
192!         ENDIF
193!       ELSE
194!         ok_pred(i)=0.
195!         ok_corr(i)=1.
196!       ENDIF
[416]197!
[445]198        ok_pred(i)=0.
199        ok_corr(i)=1.
200!
[416]201        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
202        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
[445]203!IM +++
204!       IF(n.EQ.niter) THEN
205!       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
206!         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i)
207!       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
208!         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i)
209!       ENDIF
210!       ENDIF
211!IM ---
[416]212      ENDDO
213!
214!
215!----------First aproximation of variables at zref --------------------------
216!
217      zref = 10.0
218      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
[509]219 &                 ts1, qsurf, rugos, lmon, &
220 &                 ustar, testar, qstar, zref, &
221 &                 delu, delte, delq)
[416]222!
223      DO i = 1, knon
224        u_zref(i) = delu(i)
225        q_zref(i) = max(qsurf(i),0.0) + delq(i)
226        te_zref(i) = ts1(i) + delte(i)
227        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
228!       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
229!                 (1 + RVTMP2 * max(q_zref(i),0.0))
230        u_zref_p(i) = u_zref(i)
231      ENDDO
232!
233! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995
234!
235      DO n = 1, niter
236!
237        okri=.TRUE.
238        CALL screenc(klon, knon, nsrf, zxli, &
[509]239 &                   u_zref, temp, q_zref, zref, &
240 &                   ts1, qsurf, rugos, psol, &
241 &                   ustar, testar, qstar, okri, ri1, &
242 &                   pref, delu, delte, delq)
[416]243!
244        DO i = 1, knon
245          u_zref(i) = delu(i)
246          q_zref(i) = delq(i) + max(qsurf(i),0.0)
247          te_zref(i) = delte(i) + ts1(i)
248          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
249!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
250!                   (1 + RVTMP2 * max(q_zref(i),0.0))
251        ENDDO
252!
253      ENDDO
254!
255      DO i = 1, knon
256        u_zref_c(i) = u_zref(i)
257!
258        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
259      ENDDO
260!
261      RETURN
262      END subroutine stdlevvar
Note: See TracBrowser for help on using the repository browser.