/[lmdze]/trunk/dyn3d/calfis.f90
ViewVC logotype

Contents of /trunk/dyn3d/calfis.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/calfis.f90
File size: 13896 byte(s)
Initial import
1 module calfis_m
2
3 ! Clean: no C preprocessor directive, no include line
4
5 IMPLICIT NONE
6
7 contains
8
9 SUBROUTINE calfis(nq, lafin, rdayvrai, heure, pucov, pvcov, pteta, pq, &
10 pmasse, pps, pp, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
11 clesphy0, pdufi, pdvfi, pdhfi, pdqfi, pdpsfi)
12
13 ! From dyn3d/calfis.F,v 1.3 2005/05/25 13:10:09
14
15 ! Auteurs : P. Le Van, F. Hourdin
16
17 ! 1. rearrangement des tableaux et transformation
18 ! variables dynamiques > variables physiques
19 ! 2. calcul des termes physiques
20 ! 3. retransformation des tendances physiques en tendances dynamiques
21
22 ! remarques:
23 ! ----------
24
25 ! - les vents sont donnes dans la physique par leurs composantes
26 ! naturelles.
27 ! - la variable thermodynamique de la physique est une variable
28 ! intensive : T
29 ! pour la dynamique on prend T * (preff / p(l)) **kappa
30 ! - les deux seules variables dependant de la geometrie necessaires
31 ! pour la physique sont la latitude pour le rayonnement et
32 ! l'aire de la maille quand on veut integrer une grandeur
33 ! horizontalement.
34
35 ! Input :
36 ! -------
37 ! pucov covariant zonal velocity
38 ! pvcov covariant meridional velocity
39 ! pteta potential temperature
40 ! pps surface pressure
41 ! pmasse masse d'air dans chaque maille
42 ! pts surface temperature (K)
43 ! callrad clef d'appel au rayonnement
44
45 ! Output :
46 ! --------
47 ! pdufi tendency for the natural zonal velocity (ms-1)
48 ! pdvfi tendency for the natural meridional velocity
49 ! pdhfi tendency for the potential temperature
50 ! pdtsfi tendency for the surface temperature
51
52 ! pdtrad radiative tendencies \ both input
53 ! pfluxrad radiative fluxes / and output
54
55 use dimens_m, only: iim, jjm, llm, nqmx
56 use dimphy, only: klon
57 use comconst, only: kappa, cpp, dtphys, g, pi
58 use comvert, only: preff, presnivs
59 use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
60 use advtrac_m, only: niadv
61 use grid_change, only: dyn_phy, gr_fi_dyn
62 use physiq_m, only: physiq
63
64 ! 0. Declarations :
65
66 INTEGER nq
67
68 ! Arguments :
69
70 LOGICAL, intent(in):: lafin
71 REAL, intent(in):: heure ! heure de la journée en fraction de jour
72
73 REAL pvcov(iim + 1,jjm,llm)
74 REAL pucov(iim + 1,jjm + 1,llm)
75 REAL pteta(iim + 1,jjm + 1,llm)
76 REAL pmasse(iim + 1,jjm + 1,llm)
77
78 REAL, intent(in):: pq(iim + 1,jjm + 1,llm,nqmx)
79 ! (mass fractions of advected fields)
80
81 REAL pphis(iim + 1,jjm + 1)
82 REAL pphi(iim + 1,jjm + 1,llm)
83
84 REAL pdvcov(iim + 1,jjm,llm)
85 REAL pducov(iim + 1,jjm + 1,llm)
86 REAL pdteta(iim + 1,jjm + 1,llm)
87 REAL pdq(iim + 1,jjm + 1,llm,nqmx)
88
89 REAL pw(iim + 1,jjm + 1,llm)
90
91 REAL pps(iim + 1,jjm + 1)
92 REAL pp(iim + 1,jjm + 1,llm + 1)
93 REAL ppk(iim + 1,jjm + 1,llm)
94
95 REAL pdvfi(iim + 1,jjm,llm)
96 REAL pdufi(iim + 1,jjm + 1,llm)
97 REAL pdhfi(iim + 1,jjm + 1,llm)
98 REAL pdqfi(iim + 1,jjm + 1,llm,nqmx)
99 REAL pdpsfi(iim + 1,jjm + 1)
100
101 INTEGER, PARAMETER:: longcles = 20
102 REAL clesphy0(longcles)
103
104 ! Local variables :
105
106 INTEGER i,j,l,ig0,ig,iq,iiq
107 REAL zpsrf(klon)
108 REAL zplev(klon,llm+1),zplay(klon,llm)
109 REAL zphi(klon,llm),zphis(klon)
110
111 REAL zufi(klon,llm), zvfi(klon,llm)
112 REAL ztfi(klon,llm) ! temperature
113 real zqfi(klon,llm,nqmx) ! mass fractions of advected fields
114
115 REAL pcvgu(klon,llm), pcvgv(klon,llm)
116 REAL pcvgt(klon,llm), pcvgq(klon,llm,2)
117
118 REAL pvervel(klon,llm)
119
120 REAL zdufi(klon,llm),zdvfi(klon,llm)
121 REAL zdtfi(klon,llm),zdqfi(klon,llm,nqmx)
122 REAL zdpsrf(klon)
123
124 REAL zsin(iim),zcos(iim),z1(iim)
125 REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
126 REAL unskap, pksurcp
127
128 ! I. Musat: diagnostic PVteta, Amip2
129 INTEGER, PARAMETER:: ntetaSTD=3
130 REAL:: rtetaSTD(ntetaSTD) = (/350., 380., 405./)
131 REAL PVteta(klon,ntetaSTD)
132
133 REAL SSUM
134
135 LOGICAL:: firstcal = .true.
136 REAL rdayvrai
137
138 !-----------------------------------------------------------------------
139
140 !!print *, "Call sequence information: calfis"
141
142 ! 1. Initialisations :
143 ! latitude, longitude et aires des mailles pour la physique:
144
145 ! 40. transformation des variables dynamiques en variables physiques:
146 ! 41. pressions au sol (en Pascals)
147
148 zpsrf(1) = pps(1,1)
149
150 ig0 = 2
151 DO j = 2,jjm
152 CALL SCOPY(iim,pps(1,j),1,zpsrf(ig0), 1)
153 ig0 = ig0+iim
154 ENDDO
155
156 zpsrf(klon) = pps(1,jjm + 1)
157
158 ! 42. pression intercouches :
159
160 ! .... zplev definis aux (llm +1) interfaces des couches ....
161 ! .... zplay definis aux (llm) milieux des couches ....
162
163 ! ... Exner = cp * (p(l) / preff) ** kappa ....
164
165 unskap = 1./ kappa
166
167 DO l = 1, llm + 1
168 zplev(1,l) = pp(1,1,l)
169 ig0 = 2
170 DO j = 2, jjm
171 DO i =1, iim
172 zplev(ig0,l) = pp(i,j,l)
173 ig0 = ig0 +1
174 ENDDO
175 ENDDO
176 zplev(klon,l) = pp(1,jjm + 1,l)
177 ENDDO
178
179 ! 43. temperature naturelle (en K) et pressions milieux couches .
180
181 DO l=1,llm
182
183 pksurcp = ppk(1,1,l) / cpp
184 zplay(1,l) = preff * pksurcp ** unskap
185 ztfi(1,l) = pteta(1,1,l) * pksurcp
186 pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
187 ig0 = 2
188
189 DO j = 2, jjm
190 DO i = 1, iim
191 pksurcp = ppk(i,j,l) / cpp
192 zplay(ig0,l) = preff * pksurcp ** unskap
193 ztfi(ig0,l) = pteta(i,j,l) * pksurcp
194 pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
195 ig0 = ig0 + 1
196 ENDDO
197 ENDDO
198
199 pksurcp = ppk(1,jjm + 1,l) / cpp
200 zplay(ig0,l) = preff * pksurcp ** unskap
201 ztfi (ig0,l) = pteta(1,jjm + 1,l) * pksurcp
202 pcvgt(ig0,l) = pdteta(1,jjm + 1,l) * pksurcp/ pmasse(1,jjm + 1,l)
203
204 ENDDO
205
206 ! 43.bis traceurs
207
208 DO iq=1,nq
209 iiq=niadv(iq)
210 DO l=1,llm
211 zqfi(1,l,iq) = pq(1,1,l,iiq)
212 ig0 = 2
213 DO j=2,jjm
214 DO i = 1, iim
215 zqfi(ig0,l,iq) = pq(i,j,l,iiq)
216 ig0 = ig0 + 1
217 ENDDO
218 ENDDO
219 zqfi(ig0,l,iq) = pq(1,jjm + 1,l,iiq)
220 ENDDO
221 ENDDO
222
223 ! convergence dynamique pour les traceurs "EAU"
224
225 DO iq=1,2
226 DO l=1,llm
227 pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
228 ig0 = 2
229 DO j=2,jjm
230 DO i = 1, iim
231 pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
232 ig0 = ig0 + 1
233 ENDDO
234 ENDDO
235 pcvgq(ig0,l,iq)= pdq(1,jjm + 1,l,iq) / pmasse(1,jjm + 1,l)
236 ENDDO
237 ENDDO
238
239 ! Geopotentiel calcule par rapport a la surface locale:
240
241 forall (l = 1:llm) zphi(:, l) = pack(pphi(:, :, l), dyn_phy)
242 zphis = pack(pphis, dyn_phy)
243 DO l=1,llm
244 DO ig=1,klon
245 zphi(ig,l)=zphi(ig,l)-zphis(ig)
246 ENDDO
247 ENDDO
248
249 ! .... Calcul de la vitesse verticale (en Pa*m*s ou Kg/s) ....
250
251 DO l=1,llm
252 pvervel(1,l)=pw(1,1,l) * g /apoln
253 ig0=2
254 DO j=2,jjm
255 DO i = 1, iim
256 pvervel(ig0,l) = pw(i,j,l) * g * unsaire_2d(i,j)
257 ig0 = ig0 + 1
258 ENDDO
259 ENDDO
260 pvervel(ig0,l)=pw(1,jjm + 1,l) * g /apols
261 ENDDO
262
263 ! 45. champ u:
264
265 DO l=1,llm
266
267 DO j=2,jjm
268 ig0 = 1+(j-2)*iim
269 zufi(ig0+1,l)= 0.5 * &
270 (pucov(iim,j,l)/cu_2d(iim,j) + pucov(1,j,l)/cu_2d(1,j))
271 pcvgu(ig0+1,l)= 0.5 * &
272 (pducov(iim,j,l)/cu_2d(iim,j) + pducov(1,j,l)/cu_2d(1,j))
273 DO i=2,iim
274 zufi(ig0+i,l)= 0.5 * &
275 (pucov(i-1,j,l)/cu_2d(i-1,j) &
276 + pucov(i,j,l)/cu_2d(i,j))
277 pcvgu(ig0+i,l)= 0.5 * &
278 (pducov(i-1,j,l)/cu_2d(i-1,j) &
279 + pducov(i,j,l)/cu_2d(i,j))
280 end DO
281 end DO
282
283 end DO
284
285 ! 46.champ v:
286
287 DO l=1,llm
288 DO j=2,jjm
289 ig0=1+(j-2)*iim
290 DO i=1,iim
291 zvfi(ig0+i,l)= 0.5 * &
292 (pvcov(i,j-1,l)/cv_2d(i,j-1) &
293 + pvcov(i,j,l)/cv_2d(i,j))
294 pcvgv(ig0+i,l)= 0.5 * &
295 (pdvcov(i,j-1,l)/cv_2d(i,j-1) &
296 + pdvcov(i,j,l)/cv_2d(i,j))
297 ENDDO
298 ENDDO
299 ENDDO
300
301 ! 47. champs de vents aux pole nord
302 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
303 ! V = 1 / pi * integrale [ v * sin(long) * d long ]
304
305 DO l=1,llm
306
307 z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv_2d(1,1)
308 z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv_2d(1,1)
309 DO i=2,iim
310 z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv_2d(i,1)
311 z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv_2d(i,1)
312 ENDDO
313
314 DO i=1,iim
315 zcos(i) = COS(rlonv(i))*z1(i)
316 zcosbis(i)= COS(rlonv(i))*z1bis(i)
317 zsin(i) = SIN(rlonv(i))*z1(i)
318 zsinbis(i)= SIN(rlonv(i))*z1bis(i)
319 ENDDO
320
321 zufi(1,l) = SSUM(iim,zcos,1)/pi
322 pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
323 zvfi(1,l) = SSUM(iim,zsin,1)/pi
324 pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
325
326 ENDDO
327
328 ! 48. champs de vents aux pole sud:
329 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
330 ! V = 1 / pi * integrale [ v * sin(long) * d long ]
331
332 DO l=1,llm
333
334 z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l) &
335 /cv_2d(1,jjm)
336 z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l) &
337 /cv_2d(1,jjm)
338 DO i=2,iim
339 z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv_2d(i,jjm)
340 z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv_2d(i,jjm)
341 ENDDO
342
343 DO i=1,iim
344 zcos(i) = COS(rlonv(i))*z1(i)
345 zcosbis(i) = COS(rlonv(i))*z1bis(i)
346 zsin(i) = SIN(rlonv(i))*z1(i)
347 zsinbis(i) = SIN(rlonv(i))*z1bis(i)
348 ENDDO
349
350 zufi(klon,l) = SSUM(iim,zcos,1)/pi
351 pcvgu(klon,l) = SSUM(iim,zcosbis,1)/pi
352 zvfi(klon,l) = SSUM(iim,zsin,1)/pi
353 pcvgv(klon,l) = SSUM(iim,zsinbis,1)/pi
354
355 ENDDO
356
357 !IM calcul PV a teta=350, 380, 405K
358 CALL PVtheta(klon,llm,pucov,pvcov,pteta, &
359 ztfi,zplay,zplev, &
360 ntetaSTD,rtetaSTD,PVteta)
361
362 ! Appel de la physique:
363
364 CALL physiq(nq, firstcal, lafin, rdayvrai, heure, dtphys, &
365 zplev, zplay, zphi, zphis, presnivs, clesphy0, zufi, zvfi, &
366 ztfi, zqfi, pvervel, zdufi, zdvfi, zdtfi, zdqfi, zdpsrf, pducov, &
367 PVteta) ! IM diagnostique PVteta, Amip2
368
369 ! transformation des tendances physiques en tendances dynamiques:
370
371 ! tendance sur la pression :
372
373 pdpsfi = gr_fi_dyn(zdpsrf)
374
375 ! 62. enthalpie potentielle
376
377 DO l=1,llm
378
379 DO i=1,iim + 1
380 pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
381 pdhfi(i,jjm + 1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjm + 1,l)
382 ENDDO
383
384 DO j=2,jjm
385 ig0=1+(j-2)*iim
386 DO i=1,iim
387 pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
388 ENDDO
389 pdhfi(iim + 1,j,l) = pdhfi(1,j,l)
390 ENDDO
391
392 ENDDO
393
394 ! 62. humidite specifique
395
396 DO iq=1,nqmx
397 DO l=1,llm
398 DO i=1,iim + 1
399 pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
400 pdqfi(i,jjm + 1,l,iq) = zdqfi(klon,l,iq)
401 ENDDO
402 DO j=2,jjm
403 ig0=1+(j-2)*iim
404 DO i=1,iim
405 pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
406 ENDDO
407 pdqfi(iim + 1,j,l,iq) = pdqfi(1,j,l,iq)
408 ENDDO
409 ENDDO
410 ENDDO
411
412 ! 63. traceurs
413
414 ! initialisation des tendances
415 pdqfi=0.
416
417 DO iq=1,nq
418 iiq=niadv(iq)
419 DO l=1,llm
420 DO i=1,iim + 1
421 pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
422 pdqfi(i,jjm + 1,l,iiq) = zdqfi(klon,l,iq)
423 ENDDO
424 DO j=2,jjm
425 ig0=1+(j-2)*iim
426 DO i=1,iim
427 pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
428 ENDDO
429 pdqfi(iim + 1,j,l,iiq) = pdqfi(1,j,l,iq)
430 ENDDO
431 ENDDO
432 ENDDO
433
434 ! 65. champ u:
435
436 DO l=1,llm
437
438 DO i=1,iim + 1
439 pdufi(i,1,l) = 0.
440 pdufi(i,jjm + 1,l) = 0.
441 ENDDO
442
443 DO j=2,jjm
444 ig0=1+(j-2)*iim
445 DO i=1,iim-1
446 pdufi(i,j,l)= &
447 0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu_2d(i,j)
448 ENDDO
449 pdufi(iim,j,l)= &
450 0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu_2d(iim,j)
451 pdufi(iim + 1,j,l)=pdufi(1,j,l)
452 ENDDO
453
454 ENDDO
455
456 ! 67. champ v:
457
458 DO l=1,llm
459
460 DO j=2,jjm-1
461 ig0=1+(j-2)*iim
462 DO i=1,iim
463 pdvfi(i,j,l)= &
464 0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv_2d(i,j)
465 ENDDO
466 pdvfi(iim + 1,j,l) = pdvfi(1,j,l)
467 ENDDO
468 ENDDO
469
470 ! 68. champ v pres des poles:
471 ! v = U * cos(long) + V * SIN(long)
472
473 DO l=1,llm
474
475 DO i=1,iim
476 pdvfi(i,1,l)= &
477 zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
478 pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i)) &
479 +zdvfi(klon,l)*SIN(rlonv(i))
480 pdvfi(i,1,l)= &
481 0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv_2d(i,1)
482 pdvfi(i,jjm,l)= &
483 0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iim-1+i,l))*cv_2d(i,jjm)
484 ENDDO
485
486 pdvfi(iim + 1,1,l) = pdvfi(1,1,l)
487 pdvfi(iim + 1,jjm,l)= pdvfi(1,jjm,l)
488
489 ENDDO
490
491 firstcal = .FALSE.
492
493 END SUBROUTINE calfis
494
495 end module calfis_m

  ViewVC Help
Powered by ViewVC 1.1.21