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

Last change on this file since 5320 was 2298, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes -r2237:2291 into testing branch

  • 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.8 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
5                           u1, v1, t1, q1, z1, &
6                           ts1, qsurf, z0m, z0h, 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! z0m, z0h---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, z0m, z0h
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! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag
105
106      CALL cdrag(knon, nsrf, &
107 &                   speed, t1, q1, z1, &
108 &                   psol, ts1, qsurf, z0m, z0h, &
109 &                   cdram, cdrah, zri1, pref)
110
111! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
112     IF (ok_prescr_ust) then
113      DO i = 1, knon
114       print *,'cdram avant=',cdram(i)
115       cdram(i) = ust*ust/speed(i)/speed(i)
116       print *,'cdram ust speed apres=',cdram(i),ust,speed
117      ENDDO
118     ENDIF
119!
120!---------Star variables----------------------------------------------------
121!
122      DO i = 1, knon
123        ri1(i) = zri1(i)
124        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
125        ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
126        zdte(i) = tpot(i) - ts1(i)
127        zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0)
128!
129!
130!IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
131        zdte(i) = sign(max(abs(zdte(i)),1.e-10),zdte(i))
132!
133        testar(i) = (cdrah(i) * zdte(i) * speed(i))/ustar(i)
134        qstar(i) = (cdrah(i) * zdq(i) * speed(i))/ustar(i)
135        lmon(i) = (ustar(i) * ustar(i) * tpot(i))/ &
136 &                (RKAR * RG * testar(i))
137      ENDDO
138!
139!----------First aproximation of variables at zref --------------------------
140      zref = 2.0
141      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
142 &                 ts1, qsurf, z0m, lmon, &
143 &                 ustar, testar, qstar, zref, &
144 &                 delu, delte, delq)
145!
146      DO i = 1, knon
147        u_zref(i) = delu(i)
148        q_zref(i) = max(qsurf(i),0.0) + delq(i)
149        te_zref(i) = ts1(i) + delte(i)
150        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
151        q_zref_p(i) = q_zref(i)
152!       te_zref_p(i) = te_zref(i)
153        temp_p(i) = temp(i)
154      ENDDO
155!
156! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995
157!
158      DO n = 1, niter
159!
160        okri=.TRUE.
161        CALL screenc(klon, knon, nsrf, zxli, &
162 &                   u_zref, temp, q_zref, zref, &
163 &                   ts1, qsurf, z0m, z0h, psol, &           
164 &                   ustar, testar, qstar, okri, ri1, &
165 &                   pref, delu, delte, delq)
166!
167        DO i = 1, knon
168          u_zref(i) = delu(i)
169          q_zref(i) = delq(i) + max(qsurf(i),0.0)
170          te_zref(i) = delte(i) + ts1(i)
171!
172! return to normal temperature
173!
174          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
175!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
176!                 (1 + RVTMP2 * max(q_zref(i),0.0))
177!
178!IM +++
179!         IF(temp(i).GT.350.) THEN
180!           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
181!         ENDIF
182!IM ---
183!
184        IF(n.EQ.ncon) THEN
185          te_zref_con(i) = te_zref(i)
186          q_zref_con(i) = q_zref(i)
187        ENDIF
188!
189        ENDDO
190!
191      ENDDO
192!
193! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
194!
195!       DO i = 1, knon
196!         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
197!         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
198!IM +++
199!         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
200!           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
201!           q_zref_con(i),q_zref(i),conv_q(i)
202!         ENDIF
203!IM ---
204!       ENDDO
205!
206      DO i = 1, knon
207        q_zref_c(i) = q_zref(i)
208        temp_c(i) = temp(i)
209!
210!       IF(zri1(i).LT.0.) THEN
211!         IF(nsrf.EQ.1) THEN
212!           ok_pred(i)=1.
213!           ok_corr(i)=0.
214!         ELSE
215!           ok_pred(i)=0.
216!           ok_corr(i)=1.
217!         ENDIF
218!       ELSE
219!         ok_pred(i)=0.
220!         ok_corr(i)=1.
221!       ENDIF
222!
223        ok_pred(i)=0.
224        ok_corr(i)=1.
225!
226        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
227        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
228!IM +++
229!       IF(n.EQ.niter) THEN
230!       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
231!         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i)
232!       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
233!         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i)
234!       ENDIF
235!       ENDIF
236!IM ---
237      ENDDO
238!
239!
240!----------First aproximation of variables at zref --------------------------
241!
242      zref = 10.0
243      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
244 &                 ts1, qsurf, z0m, lmon, &
245 &                 ustar, testar, qstar, zref, &
246 &                 delu, delte, delq)
247!
248      DO i = 1, knon
249        u_zref(i) = delu(i)
250        q_zref(i) = max(qsurf(i),0.0) + delq(i)
251        te_zref(i) = ts1(i) + delte(i)
252        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
253!       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
254!                 (1 + RVTMP2 * max(q_zref(i),0.0))
255        u_zref_p(i) = u_zref(i)
256      ENDDO
257!
258! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995
259!
260      DO n = 1, niter
261!
262        okri=.TRUE.
263        CALL screenc(klon, knon, nsrf, zxli, &
264 &                   u_zref, temp, q_zref, zref, &
265 &                   ts1, qsurf, z0m, z0h, psol, &
266 &                   ustar, testar, qstar, okri, ri1, &
267 &                   pref, delu, delte, delq)
268!
269        DO i = 1, knon
270          u_zref(i) = delu(i)
271          q_zref(i) = delq(i) + max(qsurf(i),0.0)
272          te_zref(i) = delte(i) + ts1(i)
273          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
274!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
275!                   (1 + RVTMP2 * max(q_zref(i),0.0))
276        ENDDO
277!
278      ENDDO
279!
280      DO i = 1, knon
281        u_zref_c(i) = u_zref(i)
282!
283        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
284!
285!AM
286        q_zref_c(i) = q_zref(i)
287        temp_c(i) = temp(i)
288        t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
289        q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
290!MA
291      ENDDO
292!
293      RETURN
294      END subroutine stdlevvar
Note: See TracBrowser for help on using the repository browser.