/[lmdze]/trunk/phylmd/stdlevvar.f
ViewVC logotype

Annotation of /trunk/phylmd/stdlevvar.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 3 months ago) by guez
Original Path: trunk/libf/phylmd/stdlevvar.f90
File size: 9533 byte(s)
Initial import
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/stdlevvar.F90,v 1.3 2005/05/25 13:10:09 fairhead Exp $
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     use YOMCST
9     use yoethf
10     IMPLICIT NONE
11     !-------------------------------------------------------------------------
12     !
13     ! Objet : calcul de la temperature et l'humidite relative a 2m et du
14     ! module du vent a 10m a partir des relations de Dyer-Businger et
15     ! des equations de Louis.
16     !
17     ! Reference : Hess, Colman et McAvaney (1995)
18     !
19     ! I. Musat, 01.07.2002
20     !
21     !AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
22     !
23     !-------------------------------------------------------------------------
24     !
25     ! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
26     ! knon----input-I- nombre de points pour un type de surface
27     ! nsrf----input-I- indice pour le type de surface; voir indicesol.inc
28     ! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
29     ! u1------input-R- vent zonal au 1er niveau du modele
30     ! v1------input-R- vent meridien au 1er niveau du modele
31     ! t1------input-R- temperature de l'air au 1er niveau du modele
32     ! q1------input-R- humidite relative au 1er niveau du modele
33     ! z1------input-R- geopotentiel au 1er niveau du modele
34     ! ts1-----input-R- temperature de l'air a la surface
35     ! qsurf---input-R- humidite relative a la surface
36     ! rugos---input-R- rugosite
37     ! psol----input-R- pression au sol
38     ! pat1----input-R- pression au 1er niveau du modele
39     !
40     ! t_2m---output-R- temperature de l'air a 2m
41     ! q_2m---output-R- humidite relative a 2m
42     ! u_10m--output-R- vitesse du vent a 10m
43     !AM
44     ! t_10m--output-R- temperature de l'air a 10m
45     ! q_10m--output-R- humidite specifique a 10m
46     ! ustar--output-R- u*
47     !
48     INTEGER, intent(in) :: klon, knon, nsrf
49     LOGICAL, intent(in) :: zxli
50     REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
51     REAL, dimension(klon), intent(in) :: qsurf, rugos
52     REAL, dimension(klon), intent(in) :: psol, pat1
53     !
54     REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar
55     REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
56     !-------------------------------------------------------------------------
57     !IM PLUS
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
77     REAL, dimension(klon) :: testar, qstar
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)
111     zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0)
112     !
113     !
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     !
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)
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
275     ENDDO
276     !
277     RETURN
278     END subroutine stdlevvar

  ViewVC Help
Powered by ViewVC 1.1.21