source: LMDZ5/branches/testing/libf/phylmd/stdlevvar.F90 @ 2176

Last change on this file since 2176 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
5                           u1, v1, t1, q1, z1, &
6                           ts1, qsurf, rugos, psol, pat1, &
7                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar)
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
18!
19!AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
20!
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 indice_sol_mod.F90
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
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*
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!
52      REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar
53      REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
54!-------------------------------------------------------------------------
55      include "flux_arp.h"
56      include "YOMCST.h"
57!IM PLUS
58      include "YOETHF.h"
59!
60! Quelques constantes et options:
61!
62! RKAR : constante de von Karman
63      REAL, PARAMETER :: RKAR=0.40
64! niter : nombre iterations calcul "corrector"
65!     INTEGER, parameter :: niter=6, ncon=niter-1
66      INTEGER, parameter :: niter=2, ncon=niter-1
67!
68! Variables locales
69      INTEGER :: i, n
70      REAL :: zref
71      REAL, dimension(klon) :: speed
72! tpot : temperature potentielle
73      REAL, dimension(klon) :: tpot
74      REAL, dimension(klon) :: zri1, cdran
75      REAL, dimension(klon) :: cdram, cdrah
76! ri1 : nb. de Richardson entre la surface --> la 1ere couche
77      REAL, dimension(klon) :: ri1
78      REAL, dimension(klon) :: testar, qstar
79      REAL, dimension(klon) :: zdte, zdq   
80! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney
81      DOUBLE PRECISION, dimension(klon) :: lmon
82      DOUBLE PRECISION, parameter :: eps=1.0D-20
83      REAL, dimension(klon) :: delu, delte, delq
84      REAL, dimension(klon) :: u_zref, te_zref, q_zref 
85      REAL, dimension(klon) :: temp, pref
86      LOGICAL :: okri
87      REAL, dimension(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
88!convertgence
89      REAL, dimension(klon) :: te_zref_con, q_zref_con
90      REAL, dimension(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
91      REAL, dimension(klon) :: ok_pred, ok_corr
92!     REAL, dimension(klon) :: conv_te, conv_q
93!-------------------------------------------------------------------------
94      DO i=1, knon
95       speed(i)=SQRT(u1(i)**2+v1(i)**2)
96       ri1(i) = 0.0
97      ENDDO
98!
99      okri=.FALSE.
100      CALL coefcdrag(klon, knon, nsrf, zxli, &
101 &                   speed, t1, q1, z1, psol, &
102 &                   ts1, qsurf, rugos, okri, ri1,  &         
103 &                   cdram, cdrah, cdran, zri1, pref)           
104! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
105     IF (ok_prescr_ust) then
106      DO i = 1, knon
107       print *,'cdram avant=',cdram(i)
108       cdram(i) = ust*ust/speed(i)/speed(i)
109       print *,'cdram ust speed apres=',cdram(i),ust,speed
110      ENDDO
111     ENDIF
112!
113!---------Star variables----------------------------------------------------
114!
115      DO i = 1, knon
116        ri1(i) = zri1(i)
117        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
118        ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
119        zdte(i) = tpot(i) - ts1(i)
120        zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0)
121!
122!
123!IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
124        zdte(i) = sign(max(abs(zdte(i)),1.e-10),zdte(i))
125!
126        testar(i) = (cdrah(i) * zdte(i) * speed(i))/ustar(i)
127        qstar(i) = (cdrah(i) * zdq(i) * speed(i))/ustar(i)
128        lmon(i) = (ustar(i) * ustar(i) * tpot(i))/ &
129 &                (RKAR * RG * testar(i))
130      ENDDO
131!
132!----------First aproximation of variables at zref --------------------------
133      zref = 2.0
134      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
135 &                 ts1, qsurf, rugos, lmon, &
136 &                 ustar, testar, qstar, zref, &
137 &                 delu, delte, delq)
138!
139      DO i = 1, knon
140        u_zref(i) = delu(i)
141        q_zref(i) = max(qsurf(i),0.0) + delq(i)
142        te_zref(i) = ts1(i) + delte(i)
143        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
144        q_zref_p(i) = q_zref(i)
145!       te_zref_p(i) = te_zref(i)
146        temp_p(i) = temp(i)
147      ENDDO
148!
149! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
150!
151      DO n = 1, niter
152!
153        okri=.TRUE.
154        CALL screenc(klon, knon, nsrf, zxli, &
155 &                   u_zref, temp, q_zref, zref, &
156 &                   ts1, qsurf, rugos, psol, &           
157 &                   ustar, testar, qstar, okri, ri1, &
158 &                   pref, delu, delte, delq)
159!
160        DO i = 1, knon
161          u_zref(i) = delu(i)
162          q_zref(i) = delq(i) + max(qsurf(i),0.0)
163          te_zref(i) = delte(i) + ts1(i)
164!
165! return to normal temperature
166!
167          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
168!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
169!                 (1 + RVTMP2 * max(q_zref(i),0.0))
170!
171!IM +++
172!         IF(temp(i).GT.350.) THEN
173!           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
174!         ENDIF
175!IM ---
176!
177        IF(n.EQ.ncon) THEN
178          te_zref_con(i) = te_zref(i)
179          q_zref_con(i) = q_zref(i)
180        ENDIF
181!
182        ENDDO
183!
184      ENDDO
185!
186! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
187!
188!       DO i = 1, knon
189!         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
190!         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
191!IM +++
192!         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
193!           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
194!           q_zref_con(i),q_zref(i),conv_q(i)
195!         ENDIF
196!IM ---
197!       ENDDO
198!
199      DO i = 1, knon
200        q_zref_c(i) = q_zref(i)
201        temp_c(i) = temp(i)
202!
203!       IF(zri1(i).LT.0.) THEN
204!         IF(nsrf.EQ.1) THEN
205!           ok_pred(i)=1.
206!           ok_corr(i)=0.
207!         ELSE
208!           ok_pred(i)=0.
209!           ok_corr(i)=1.
210!         ENDIF
211!       ELSE
212!         ok_pred(i)=0.
213!         ok_corr(i)=1.
214!       ENDIF
215!
216        ok_pred(i)=0.
217        ok_corr(i)=1.
218!
219        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
220        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
221!IM +++
222!       IF(n.EQ.niter) THEN
223!       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
224!         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i)
225!       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
226!         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i)
227!       ENDIF
228!       ENDIF
229!IM ---
230      ENDDO
231!
232!
233!----------First aproximation of variables at zref --------------------------
234!
235      zref = 10.0
236      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
237 &                 ts1, qsurf, rugos, lmon, &
238 &                 ustar, testar, qstar, zref, &
239 &                 delu, delte, delq)
240!
241      DO i = 1, knon
242        u_zref(i) = delu(i)
243        q_zref(i) = max(qsurf(i),0.0) + delq(i)
244        te_zref(i) = ts1(i) + delte(i)
245        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
246!       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
247!                 (1 + RVTMP2 * max(q_zref(i),0.0))
248        u_zref_p(i) = u_zref(i)
249      ENDDO
250!
251! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995
252!
253      DO n = 1, niter
254!
255        okri=.TRUE.
256        CALL screenc(klon, knon, nsrf, zxli, &
257 &                   u_zref, temp, q_zref, zref, &
258 &                   ts1, qsurf, rugos, psol, &
259 &                   ustar, testar, qstar, okri, ri1, &
260 &                   pref, delu, delte, delq)
261!
262        DO i = 1, knon
263          u_zref(i) = delu(i)
264          q_zref(i) = delq(i) + max(qsurf(i),0.0)
265          te_zref(i) = delte(i) + ts1(i)
266          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
267!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
268!                   (1 + RVTMP2 * max(q_zref(i),0.0))
269        ENDDO
270!
271      ENDDO
272!
273      DO i = 1, knon
274        u_zref_c(i) = u_zref(i)
275!
276        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
277!
278!AM
279        q_zref_c(i) = q_zref(i)
280        temp_c(i) = temp(i)
281        t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
282        q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
283!MA
284      ENDDO
285!
286      RETURN
287      END subroutine stdlevvar
Note: See TracBrowser for help on using the repository browser.