/[lmdze]/trunk/libf/dyn3d/guide.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/guide.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 29 - (show annotations)
Tue Mar 30 10:44:42 2010 UTC (14 years, 1 month ago) by guez
File size: 14872 byte(s)
In "leapfrog", transformed some arrays with a single dimension for horizontal
position into arrays with two horizontal dimensions. Simplified the
computation of potential temperature and surface pressure at the
poles.

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 REAL :: tau_min_u, tau_max_u
7 REAL :: tau_min_v, tau_max_v
8 REAL :: tau_min_t, tau_max_t
9 REAL :: tau_min_q, tau_max_q
10 REAL :: tau_min_p, tau_max_p
11 REAL :: aire_min, aire_max
12
13
14 LOGICAL :: guide_u, guide_v, guide_t, guide_q, guide_p
15 REAL :: lat_min_guide, lat_max_guide
16
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 dimens_m, ONLY : jjm, llm
27 USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
28 USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
29 USE comvert, ONLY : ap, bp, preff, presnivs
30 USE conf_gcm_m, ONLY : day_step, iperiod
31 USE comgeom, ONLY : aire, rlatu, rlonv
32 USE serre, ONLY : clat, clon
33 USE q_sat_m, ONLY : q_sat
34 USE exner_hyb_m, ONLY : exner_hyb
35 USE pression_m, ONLY : pression
36 USE inigrads_m, ONLY : inigrads
37 use netcdf, only: nf90_nowrite, nf90_open, nf90_close
38
39 IMPLICIT NONE
40
41 INCLUDE 'netcdf.inc'
42
43 ! variables dynamiques
44 REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
45 REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
46 REAL :: q(ip1jmp1, llm) ! temperature potentielle
47 REAL :: ps(ip1jmp1) ! pression au sol
48 REAL :: masse(ip1jmp1, llm) ! masse d'air
49
50 ! common passe pour des sorties
51 REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
52 COMMON /comdxdy/dxdys, dxdyu, dxdyv
53
54 ! variables dynamiques pour les reanalyses.
55 REAL :: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
56 REAL :: tetarea1(ip1jmp1, llm) ! temp pot reales
57 REAL :: qrea1(ip1jmp1, llm) ! temp pot reales
58 REAL :: psrea1(ip1jmp1) ! ps
59 REAL :: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
60 REAL :: tetarea2(ip1jmp1, llm) ! temp pot reales
61 REAL :: qrea2(ip1jmp1, llm) ! temp pot reales
62 REAL :: masserea2(ip1jmp1, llm) ! masse
63 REAL :: psrea2(ip1jmp1) ! ps
64
65 REAL :: alpha_q(ip1jmp1)
66 REAL :: alpha_t(ip1jmp1), alpha_p(ip1jmp1)
67 REAL :: alpha_u(ip1jmp1), alpha_v(ip1jm)
68 REAL :: dday_step, toto, reste, itau_test
69 INTEGER :: step_rea, count_no_rea
70
71 INTEGER :: ilon, ilat
72 REAL :: factt, ztau(ip1jmp1)
73
74 INTEGER, INTENT (IN) :: itau
75 INTEGER :: ij, l
76 INTEGER :: ncidpl, varidpl, nlev, status
77 INTEGER :: rcod, rid
78 REAL :: ditau, tau, a
79 SAVE nlev
80
81 ! TEST SUR QSAT
82 REAL :: p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
83 REAL :: pkf(ip1jmp1, llm)
84 REAL :: pres(ip1jmp1, llm)
85
86 REAL :: qsat(ip1jmp1, llm)
87 REAL :: unskap
88 REAL :: tnat(ip1jmp1, llm)
89
90
91 LOGICAL :: first
92 SAVE first
93 DATA first/ .TRUE./
94
95 SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
96 SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
97
98 SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
99 SAVE step_rea, count_no_rea
100
101 CHARACTER (10) :: file
102 INTEGER :: igrads
103 REAL :: dtgrads
104 SAVE igrads, dtgrads
105 DATA igrads, dtgrads/2, 100./
106
107 !-----------------------------------------------------------------------
108
109 PRINT *, 'Call sequence information: guide'
110
111 ! calcul de l'humidite saturante
112
113 CALL pression(ip1jmp1, ap, bp, ps, p)
114 CALL massdair(p, masse)
115 PRINT *, 'OK1'
116 CALL exner_hyb(ps, p, pks, pk, pkf)
117 PRINT *, 'OK2'
118 tnat(:, :) = pk(:, :)*teta(:, :)/cpp
119 PRINT *, 'OK3'
120 unskap = 1./kappa
121 pres(:, :) = preff*(pk(:, :)/cpp)**unskap
122 PRINT *, 'OK4'
123 qsat = q_sat(tnat, pres)
124
125 ! initialisations pour la lecture des reanalyses.
126 ! alpha determine la part des injections de donnees a chaque etape
127 ! alpha=1 signifie pas d'injection
128 ! alpha=0 signifie injection totale
129
130 PRINT *, 'ONLINE=', online
131 IF (online==-1) THEN
132 RETURN
133 END IF
134
135 IF (first) THEN
136
137 PRINT *, 'initialisation du guide '
138 CALL conf_guide
139 PRINT *, 'apres conf_guide'
140
141 file = 'guide'
142 CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
143 180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
144
145 PRINT *, &
146 '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
147
148 IF (online==-1) RETURN
149 IF (online==1) THEN
150
151 ! Constantes de temps de rappel en jour
152 ! 0.1 c'est en gros 2h30.
153 ! 1e10 est une constante infinie donc en gros pas de guidage
154
155 ! coordonnees du centre du zoom
156 CALL coordij(clon, clat, ilon, ilat)
157 ! aire de la maille au centre du zoom
158 aire_min = aire(ilon+(ilat-1)*iip1)
159 ! aire maximale de la maille
160 aire_max = 0.
161 DO ij = 1, ip1jmp1
162 aire_max = max(aire_max, aire(ij))
163 END DO
164 ! factt = pas de temps en fraction de jour
165 factt = dtvr*iperiod/daysec
166
167 CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
168 CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
169 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
170 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
171 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
172
173 CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
174 CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
175 CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
176
177 ! Cas ou on force exactement par les variables analysees
178 ELSE
179 alpha_t = 0.
180 alpha_u = 0.
181 alpha_v = 0.
182 alpha_p = 0.
183 ! physic=.false.
184 END IF
185
186 itau_test = 1001
187 step_rea = 1
188 count_no_rea = 0
189 ncidpl = -99
190
191 ! itau_test montre si l'importation a deja ete faite au rang itau
192 ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
193 if (guide_u) then
194 if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
195 endif
196
197 if (guide_v) then
198 if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
199 endif
200
201 if (guide_T) then
202 if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
203 endif
204
205 if (guide_Q) then
206 if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
207 endif
208
209 IF (ncep) THEN
210 status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
211 ELSE
212 status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
213 END IF
214 status = nf_inq_dimlen(ncidpl, rid, nlev)
215 PRINT *, 'nlev', nlev
216 rcod = nf90_close(ncidpl)
217 ! Lecture du premier etat des reanalyses.
218 CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
219 masserea2, psrea2, 1, nlev)
220 qrea2(:, :) = max(qrea2(:, :), 0.1)
221
222
223 ! Debut de l'integration temporelle:
224 END IF ! first
225
226 ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
227
228 ditau = real(itau)
229 dday_step = real(day_step)
230 WRITE (*, *) 'ditau, dday_step'
231 WRITE (*, *) ditau, dday_step
232 toto = 4*ditau/dday_step
233 reste = toto - aint(toto)
234
235 IF (reste==0.) THEN
236 IF (itau_test==itau) THEN
237 WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
238 STOP
239 ELSE
240 vcovrea1(:, :) = vcovrea2(:, :)
241 ucovrea1(:, :) = ucovrea2(:, :)
242 tetarea1(:, :) = tetarea2(:, :)
243 qrea1(:, :) = qrea2(:, :)
244
245 PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
246 count_no_rea, ' non lectures'
247 step_rea = step_rea + 1
248 itau_test = itau
249 CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
250 qrea2, masserea2, psrea2, 1, nlev)
251 qrea2(:, :) = max(qrea2(:, :), 0.1)
252 factt = dtvr*iperiod/daysec
253 ztau(:) = factt/max(alpha_t(:), 1.E-10)
254 CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
255 CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
256 CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
257 CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
258 CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
259 CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
260 CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
261 CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
262 CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
263 CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
264 CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
265
266 CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
267
268 END IF
269 ELSE
270 count_no_rea = count_no_rea + 1
271 END IF
272
273 ! Guidage
274 ! x_gcm = a * x_gcm + (1-a) * x_reanalyses
275
276 IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
277
278 ditau = real(itau)
279 dday_step = real(day_step)
280
281
282 tau = 4*ditau/dday_step
283 tau = tau - aint(tau)
284
285 ! ucov
286 IF (guide_u) THEN
287 DO l = 1, llm
288 DO ij = 1, ip1jmp1
289 a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
290 ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
291 IF (first .AND. ini_anal) ucov(ij, l) = a
292 END DO
293 END DO
294 END IF
295
296 IF (guide_t) THEN
297 DO l = 1, llm
298 DO ij = 1, ip1jmp1
299 a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
300 teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
301 IF (first .AND. ini_anal) teta(ij, l) = a
302 END DO
303 END DO
304 END IF
305
306 ! P
307 IF (guide_p) THEN
308 DO ij = 1, ip1jmp1
309 a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
310 ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
311 IF (first .AND. ini_anal) ps(ij) = a
312 END DO
313 CALL pression(ip1jmp1, ap, bp, ps, p)
314 CALL massdair(p, masse)
315 END IF
316
317
318 ! q
319 IF (guide_q) THEN
320 DO l = 1, llm
321 DO ij = 1, ip1jmp1
322 a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
323 ! hum relative en % -> hum specif
324 a = qsat(ij, l)*a*0.01
325 q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
326 IF (first .AND. ini_anal) q(ij, l) = a
327 END DO
328 END DO
329 END IF
330
331 ! vcov
332 IF (guide_v) THEN
333 DO l = 1, llm
334 DO ij = 1, ip1jm
335 a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
336 vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
337 IF (first .AND. ini_anal) vcov(ij, l) = a
338 END DO
339 IF (first .AND. ini_anal) vcov(ij, l) = a
340 END DO
341 END IF
342
343 first = .FALSE.
344
345 END SUBROUTINE guide
346
347 !=======================================================================
348 SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)
349 !=======================================================================
350
351 USE dimens_m, ONLY : iim, jjm
352 USE paramet_m, ONLY : iip1, jjp1
353 USE comconst, ONLY : pi
354 USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv
355 USE serre, ONLY : clat, clon, grossismx, grossismy
356 IMPLICIT NONE
357
358 ! arguments :
359 INTEGER :: type
360 INTEGER :: pim, pjm
361 REAL :: factt, taumin, taumax
362 REAL :: dxdy_, alpha(pim, pjm)
363 REAL :: dxdy_min, dxdy_max
364
365 ! local :
366 REAL :: alphamin, alphamax, gamma, xi
367 SAVE gamma
368 INTEGER :: i, j, ilon, ilat
369
370 LOGICAL :: first
371 SAVE first
372 DATA first/ .TRUE./
373
374 REAL :: zdx(iip1, jjp1), zdy(iip1, jjp1)
375
376 REAL :: zlat
377 REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
378 COMMON /comdxdy/dxdys, dxdyu, dxdyv
379
380 IF (first) THEN
381 DO j = 2, jjm
382 DO i = 2, iip1
383 zdx(i, j) = 0.5*(cu_2d(i-1, j)+cu_2d(i, j))/cos(rlatu(j))
384 END DO
385 zdx(1, j) = zdx(iip1, j)
386 END DO
387 DO j = 2, jjm
388 DO i = 1, iip1
389 zdy(i, j) = 0.5*(cv_2d(i, j-1)+cv_2d(i, j))
390 END DO
391 END DO
392 DO i = 1, iip1
393 zdx(i, 1) = zdx(i, 2)
394 zdx(i, jjp1) = zdx(i, jjm)
395 zdy(i, 1) = zdy(i, 2)
396 zdy(i, jjp1) = zdy(i, jjm)
397 END DO
398 DO j = 1, jjp1
399 DO i = 1, iip1
400 dxdys(i, j) = sqrt(zdx(i, j)*zdx(i, j)+zdy(i, j)*zdy(i, j))
401 END DO
402 END DO
403 DO j = 1, jjp1
404 DO i = 1, iim
405 dxdyu(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
406 END DO
407 dxdyu(iip1, j) = dxdyu(1, j)
408 END DO
409 DO j = 1, jjm
410 DO i = 1, iip1
411 dxdyv(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
412 END DO
413 END DO
414
415 CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL ')
416 CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U ')
417 CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v ')
418
419 ! coordonnees du centre du zoom
420 CALL coordij(clon, clat, ilon, ilat)
421 ! aire de la maille au centre du zoom
422 dxdy_min = dxdys(ilon, ilat)
423 ! dxdy maximale de la maille
424 dxdy_max = 0.
425 DO j = 1, jjp1
426 DO i = 1, iip1
427 dxdy_max = max(dxdy_max, dxdys(i, j))
428 END DO
429 END DO
430
431 IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
432 PRINT *, 'ATTENTION modele peu zoome'
433 PRINT *, 'ATTENTION on prend une constante de guidage cste'
434 gamma = 0.
435 ELSE
436 gamma = (dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
437 PRINT *, 'gamma=', gamma
438 IF (gamma<1.E-5) THEN
439 PRINT *, 'gamma =', gamma, '<1e-5'
440 STOP
441 END IF
442 PRINT *, 'gamma=', gamma
443 gamma = log(0.5)/log(gamma)
444 END IF
445 END IF
446
447 alphamin = factt/taumax
448 alphamax = factt/taumin
449
450 DO j = 1, pjm
451 DO i = 1, pim
452 IF (type==1) THEN
453 dxdy_ = dxdys(i, j)
454 zlat = rlatu(j)*180./pi
455 ELSE IF (type==2) THEN
456 dxdy_ = dxdyu(i, j)
457 zlat = rlatu(j)*180./pi
458 ELSE IF (type==3) THEN
459 dxdy_ = dxdyv(i, j)
460 zlat = rlatv(j)*180./pi
461 END IF
462 IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
463 ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
464 alpha(i, j) = alphamin
465 ELSE
466 xi = ((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
467 xi = min(xi, 1.)
468 IF (lat_min_guide<=zlat .AND. zlat<=lat_max_guide) THEN
469 alpha(i, j) = xi*alphamin + (1.-xi)*alphamax
470 ELSE
471 alpha(i, j) = 0.
472 END IF
473 END IF
474 END DO
475 END DO
476
477
478 RETURN
479 END SUBROUTINE tau2alpha
480
481 END MODULE guide_m

  ViewVC Help
Powered by ViewVC 1.1.21