source: LMDZ4/trunk/libf/phylmd/stdlevvar.F90 @ 728

Last change on this file since 728 was 644, checked in by Laurent Fairhead, 20 years ago

Synchronisation avec tous les diagnostiques de Ionela IM
Inclusion du slab ocean IM
LF

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