/[lmdze]/trunk/Sources/dyn3d/Guide/guide.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/Guide/guide.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (show annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/guide.f90
File size: 10762 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

1 MODULE guide_m
2
3 ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4 ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5
6 IMPLICIT NONE
7
8 REAL tau_min_u, tau_max_u
9 REAL tau_min_v, tau_max_v
10 REAL tau_min_t, tau_max_t
11 REAL tau_min_q, tau_max_q
12 REAL tau_min_p, tau_max_p
13 REAL aire_min, aire_max
14
15
16 LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p
17 LOGICAL ncep, ini_anal
18 INTEGER online
19
20 CONTAINS
21
22 SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
23
24 ! Author: F.Hourdin
25
26 USE comconst, ONLY : cpp, daysec, dtvr, kappa
27 USE comvert, ONLY : ap, bp, preff, presnivs
28 USE conf_gcm_m, ONLY : day_step, iperiod
29 USE comgeom, ONLY : aire, rlatu, rlonv
30 USE dimens_m, ONLY : jjm, llm
31 USE exner_hyb_m, ONLY : exner_hyb
32 USE inigrads_m, ONLY : inigrads
33 use netcdf, only: nf90_nowrite, nf90_open, nf90_close
34 use nr_util, only: pi
35 USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
36 USE q_sat_m, ONLY : q_sat
37 USE serre, ONLY : clat, clon
38 use tau2alpha_m, only: tau2alpha
39
40 INCLUDE 'netcdf.inc'
41
42 ! variables dynamiques
43 REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
44 REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
45 REAL q(ip1jmp1, llm) ! temperature potentielle
46 REAL ps(ip1jmp1) ! pression au sol
47 REAL masse(ip1jmp1, llm) ! masse d'air
48
49 ! common passe pour des sorties
50 REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
51 COMMON /comdxdy/dxdys, dxdyu, dxdyv
52
53 ! variables dynamiques pour les reanalyses.
54 REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
55 REAL tetarea1(ip1jmp1, llm) ! temp pot reales
56 REAL qrea1(ip1jmp1, llm) ! temp pot reales
57 REAL psrea1(ip1jmp1) ! ps
58 REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
59 REAL tetarea2(ip1jmp1, llm) ! temp pot reales
60 REAL qrea2(ip1jmp1, llm) ! temp pot reales
61 REAL masserea2(ip1jmp1, llm) ! masse
62 REAL psrea2(ip1jmp1) ! ps
63
64 REAL alpha_q(ip1jmp1)
65 REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)
66 REAL alpha_u(ip1jmp1), alpha_v(ip1jm)
67 REAL dday_step, toto, reste, itau_test
68 INTEGER step_rea, count_no_rea
69
70 INTEGER ilon, ilat
71 REAL factt, ztau(ip1jmp1)
72
73 INTEGER, INTENT (IN) :: itau
74 INTEGER ij, l
75 INTEGER ncidpl, varidpl, nlev, status
76 INTEGER rcod, rid
77 REAL ditau, tau, a
78 SAVE nlev
79
80 ! TEST SUR QSAT
81 REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
82 REAL pkf(ip1jmp1, llm)
83 REAL pres(ip1jmp1, llm)
84
85 REAL qsat(ip1jmp1, llm)
86 REAL unskap
87 REAL tnat(ip1jmp1, llm)
88
89 LOGICAL:: first = .TRUE.
90
91 SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
92 SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
93
94 SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
95 SAVE step_rea, count_no_rea
96
97 CHARACTER (10) file
98 INTEGER igrads
99 REAL dtgrads
100 SAVE igrads, dtgrads
101 DATA igrads, dtgrads/2, 100./
102
103 !-----------------------------------------------------------------------
104
105 PRINT *, 'Call sequence information: guide'
106
107 ! calcul de l'humidite saturante
108
109 forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
110 CALL massdair(p, masse)
111 PRINT *, 'OK1'
112 CALL exner_hyb(ps, p, pks, pk, pkf)
113 PRINT *, 'OK2'
114 tnat(:, :) = pk(:, :)*teta(:, :)/cpp
115 PRINT *, 'OK3'
116 unskap = 1./kappa
117 pres(:, :) = preff*(pk(:, :)/cpp)**unskap
118 PRINT *, 'OK4'
119 qsat = q_sat(tnat, pres)
120
121 ! initialisations pour la lecture des reanalyses.
122 ! alpha determine la part des injections de donnees a chaque etape
123 ! alpha=1 signifie pas d'injection
124 ! alpha=0 signifie injection totale
125
126 PRINT *, 'ONLINE=', online
127 IF (online==-1) THEN
128 RETURN
129 END IF
130
131 IF (first) THEN
132
133 PRINT *, 'initialisation du guide '
134 CALL conf_guide
135 PRINT *, 'apres conf_guide'
136
137 file = 'guide'
138 CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
139 180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
140
141 PRINT *, &
142 '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
143
144 IF (online==-1) RETURN
145 IF (online==1) THEN
146
147 ! Constantes de temps de rappel en jour
148 ! 0.1 c'est en gros 2h30.
149 ! 1e10 est une constante infinie donc en gros pas de guidage
150
151 ! coordonnees du centre du zoom
152 CALL coordij(clon, clat, ilon, ilat)
153 ! aire de la maille au centre du zoom
154 aire_min = aire(ilon+(ilat-1)*iip1)
155 ! aire maximale de la maille
156 aire_max = 0.
157 DO ij = 1, ip1jmp1
158 aire_max = max(aire_max, aire(ij))
159 END DO
160 ! factt = pas de temps en fraction de jour
161 factt = dtvr*iperiod/daysec
162
163 CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
164 CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
165 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
166 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
167 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
168
169 CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
170 CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
171 CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
172
173 ! Cas ou on force exactement par les variables analysees
174 ELSE
175 alpha_t = 0.
176 alpha_u = 0.
177 alpha_v = 0.
178 alpha_p = 0.
179 ! physic=.false.
180 END IF
181
182 itau_test = 1001
183 step_rea = 1
184 count_no_rea = 0
185 ncidpl = -99
186
187 ! itau_test montre si l'importation a deja ete faite au rang itau
188 ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
189 if (guide_u) then
190 if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
191 endif
192
193 if (guide_v) then
194 if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
195 endif
196
197 if (guide_T) then
198 if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
199 endif
200
201 if (guide_Q) then
202 if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
203 endif
204
205 IF (ncep) THEN
206 status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
207 ELSE
208 status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
209 END IF
210 status = nf_inq_dimlen(ncidpl, rid, nlev)
211 PRINT *, 'nlev', nlev
212 rcod = nf90_close(ncidpl)
213 ! Lecture du premier etat des reanalyses.
214 CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
215 masserea2, psrea2, 1, nlev)
216 qrea2(:, :) = max(qrea2(:, :), 0.1)
217
218
219 ! Debut de l'integration temporelle:
220 END IF ! first
221
222 ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
223
224 ditau = real(itau)
225 dday_step = real(day_step)
226 WRITE (*, *) 'ditau, dday_step'
227 WRITE (*, *) ditau, dday_step
228 toto = 4*ditau/dday_step
229 reste = toto - aint(toto)
230
231 IF (reste==0.) THEN
232 IF (itau_test==itau) THEN
233 WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
234 STOP
235 ELSE
236 vcovrea1(:, :) = vcovrea2(:, :)
237 ucovrea1(:, :) = ucovrea2(:, :)
238 tetarea1(:, :) = tetarea2(:, :)
239 qrea1(:, :) = qrea2(:, :)
240
241 PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
242 count_no_rea, ' non lectures'
243 step_rea = step_rea + 1
244 itau_test = itau
245 CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
246 qrea2, masserea2, psrea2, 1, nlev)
247 qrea2(:, :) = max(qrea2(:, :), 0.1)
248 factt = dtvr*iperiod/daysec
249 ztau(:) = factt/max(alpha_t(:), 1.E-10)
250 CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
251 CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
252 CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
253 CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
254 CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
255 CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
256 CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
257 CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
258 CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
259 CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
260 CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
261
262 CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
263
264 END IF
265 ELSE
266 count_no_rea = count_no_rea + 1
267 END IF
268
269 ! Guidage
270 ! x_gcm = a * x_gcm + (1-a) * x_reanalyses
271
272 IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
273
274 ditau = real(itau)
275 dday_step = real(day_step)
276
277
278 tau = 4*ditau/dday_step
279 tau = tau - aint(tau)
280
281 ! ucov
282 IF (guide_u) THEN
283 DO l = 1, llm
284 DO ij = 1, ip1jmp1
285 a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
286 ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
287 IF (first .AND. ini_anal) ucov(ij, l) = a
288 END DO
289 END DO
290 END IF
291
292 IF (guide_t) THEN
293 DO l = 1, llm
294 DO ij = 1, ip1jmp1
295 a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
296 teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
297 IF (first .AND. ini_anal) teta(ij, l) = a
298 END DO
299 END DO
300 END IF
301
302 ! P
303 IF (guide_p) THEN
304 DO ij = 1, ip1jmp1
305 a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
306 ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
307 IF (first .AND. ini_anal) ps(ij) = a
308 END DO
309 forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
310 CALL massdair(p, masse)
311 END IF
312
313
314 ! q
315 IF (guide_q) THEN
316 DO l = 1, llm
317 DO ij = 1, ip1jmp1
318 a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
319 ! hum relative en % -> hum specif
320 a = qsat(ij, l)*a*0.01
321 q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
322 IF (first .AND. ini_anal) q(ij, l) = a
323 END DO
324 END DO
325 END IF
326
327 ! vcov
328 IF (guide_v) THEN
329 DO l = 1, llm
330 DO ij = 1, ip1jm
331 a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
332 vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
333 IF (first .AND. ini_anal) vcov(ij, l) = a
334 END DO
335 IF (first .AND. ini_anal) vcov(ij, l) = a
336 END DO
337 END IF
338
339 first = .FALSE.
340
341 END SUBROUTINE guide
342
343 END MODULE guide_m

  ViewVC Help
Powered by ViewVC 1.1.21