source: LMDZ.3.3/trunk/libf/dyn3d/calfis.F @ 221

Last change on this file since 221 was 221, checked in by lmdz, 23 years ago

Bug sur la vitesse verticale. P LeVan?
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.4 KB
Line 
1      SUBROUTINE calfis(nq, lafin, rdayvrai,rday_ecri, heure,
2     $            pucov,pvcov,pteta,pq,pmasse,pps,pp,ppk,pphis,pphi,
3     $            pducov,pdvcov,pdteta,pdq,pw, clesphy0,
4     $            pdufi,pdvfi,pdhfi,pdqfi,pdpsfi                   )
5c
6c    Auteur :  P. Le Van, F. Hourdin
7c   .........
8
9      IMPLICIT NONE
10c=======================================================================
11c
12c   1. rearrangement des tableaux et transformation
13c      variables dynamiques  >  variables physiques
14c   2. calcul des termes physiques
15c   3. retransformation des tendances physiques en tendances dynamiques
16c
17c   remarques:
18c   ----------
19c
20c    - les vents sont donnes dans la physique par leurs composantes
21c      naturelles.
22c    - la variable thermodynamique de la physique est une variable
23c      intensive :   T
24c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
25c    - les deux seules variables dependant de la geometrie necessaires
26c      pour la physique sont la latitude pour le rayonnement et
27c      l'aire de la maille quand on veut integrer une grandeur
28c      horizontalement.
29c    - les points de la physique sont les points scalaires de la
30c      la dynamique; numerotation:
31c          1 pour le pole nord
32c          (jjm-1)*iim pour l'interieur du domaine
33c          ngridmx pour le pole sud
34c      ---> ngridmx=2+(jjm-1)*iim
35c
36c     Input :
37c     -------
38c       ecritphy        frequence d'ecriture (en jours)de histphy
39c       pucov           covariant zonal velocity
40c       pvcov           covariant meridional velocity
41c       pteta           potential temperature
42c       pps             surface pressure
43c       pmasse          masse d'air dans chaque maille
44c       pts             surface temperature  (K)
45c       callrad         clef d'appel au rayonnement
46c
47c    Output :
48c    --------
49c        pdufi          tendency for the natural zonal velocity (ms-1)
50c        pdvfi          tendency for the natural meridional velocity
51c        pdhfi          tendency for the potential temperature
52c        pdtsfi         tendency for the surface temperature
53c
54c        pdtrad         radiative tendencies  \  both input
55c        pfluxrad       radiative fluxes      /  and output
56c
57c=======================================================================
58c
59c-----------------------------------------------------------------------
60c
61c    0.  Declarations :
62c    ------------------
63
64#include "dimensions.h"
65#include "paramet.h"
66#include "temps.h"
67
68      INTEGER ngridmx,nq
69      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
70
71#include "comconst.h"
72#include "comvert.h"
73#include "comgeom2.h"
74#include "control.h"
75
76c    Arguments :
77c    -----------
78      LOGICAL  lafin
79      REAL heure
80
81      REAL pvcov(iip1,jjm,llm)
82      REAL pucov(iip1,jjp1,llm)
83      REAL pteta(iip1,jjp1,llm)
84      REAL pmasse(iip1,jjp1,llm)
85      REAL pq(iip1,jjp1,llm,nqmx)
86      REAL pphis(iip1,jjp1)
87      REAL pphi(iip1,jjp1,llm)
88c
89      REAL pdvcov(iip1,jjm,llm)
90      REAL pducov(iip1,jjp1,llm)
91      REAL pdteta(iip1,jjp1,llm)
92      REAL pdq(iip1,jjp1,llm,nqmx)
93c
94      REAL pw(iip1,jjp1,llm)
95c
96      REAL pps(iip1,jjp1)
97      REAL pp(iip1,jjp1,llmp1)
98      REAL ppk(iip1,jjp1,llm)
99c
100      REAL pdvfi(iip1,jjm,llm)
101      REAL pdufi(iip1,jjp1,llm)
102      REAL pdhfi(iip1,jjp1,llm)
103      REAL pdqfi(iip1,jjp1,llm,nqmx)
104      REAL pdpsfi(iip1,jjp1)
105
106      INTEGER        longcles
107      PARAMETER    ( longcles = 20 )
108      REAL clesphy0( longcles )
109
110
111c    Local variables :
112c    -----------------
113
114      INTEGER i,j,l,ig0,ig,iq
115      REAL zpsrf(ngridmx)
116      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
117      REAL zphi(ngridmx,llm),zphis(ngridmx)
118c
119      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
120      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
121c
122      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
123      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
124c
125      REAL pvervel(ngridmx,llm)
126c
127      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
128      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
129      REAL zdpsrf(ngridmx)
130c
131      REAL zsin(iim),zcos(iim),z1(iim)
132      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
133      REAL unskap, pksurcp
134c
135     
136      EXTERNAL gr_dyn_fi,gr_fi_dyn
137      EXTERNAL physiq,multipl
138      REAL SSUM
139      EXTERNAL SSUM
140
141      REAL latfi(ngridmx),lonfi(ngridmx)
142      REAL airefi(ngridmx)
143      SAVE latfi, lonfi, airefi
144
145      LOGICAL firstcal, debut
146      DATA firstcal/.true./
147      SAVE firstcal,debut
148      REAL rdayvrai,rday_ecri
149c
150c-----------------------------------------------------------------------
151c
152c    1. Initialisations :
153c    --------------------
154c
155
156      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
157         PRINT*,'STOP dans calfis'
158         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
159         PRINT*,'  ngridmx  jjm   iim   '
160         PRINT*,ngridmx,jjm,iim
161         STOP
162      ENDIF
163
164c-----------------------------------------------------------------------
165c   latitude, longitude et aires des mailles pour la physique:
166c   ----------------------------------------------------------
167
168c
169      IF ( firstcal )  THEN
170          debut = .TRUE.
171      ELSE
172          debut = .FALSE.
173      ENDIF
174
175c
176      IF (firstcal) THEN
177         latfi(1)=rlatu(1)
178         lonfi(1)=0.
179         DO j=2,jjm
180            DO i=1,iim
181               latfi((j-2)*iim+1+i)= rlatu(j)
182               lonfi((j-2)*iim+1+i)= rlonv(i)
183            ENDDO
184         ENDDO
185         latfi(ngridmx)= rlatu(jjp1)
186         lonfi(ngridmx)= 0.
187         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
188          PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
189         CALL inifis(ngridmx,llm,daysec,day_ini,dtphys ,
190     ,                latfi,lonfi,airefi,rad,g,r,cpp     )
191          ENDIF
192
193c
194c-----------------------------------------------------------------------
195c   40. transformation des variables dynamiques en variables physiques:
196c   ---------------------------------------------------------------
197
198c   41. pressions au sol (en Pascals)
199c   ----------------------------------
200
201       
202      zpsrf(1) = pps(1,1)
203
204      ig0  = 2
205      DO j = 2,jjm
206         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
207         ig0 = ig0+iim
208      ENDDO
209
210      zpsrf(ngridmx) = pps(1,jjp1)
211
212
213c   42. pression intercouches :
214c
215c   -----------------------------------------------------------------
216c     .... zplev  definis aux (llm +1) interfaces des couches  ....
217c     .... zplay  definis aux (  llm )    milieux des couches  ....
218c   -----------------------------------------------------------------
219
220c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
221c
222       unskap   = 1./ kappa
223c
224      DO l = 1, llmp1
225        zplev( 1,l ) = pp(1,1,l)
226        ig0 = 2
227          DO j = 2, jjm
228             DO i =1, iim
229              zplev( ig0,l ) = pp(i,j,l)
230              ig0 = ig0 +1
231             ENDDO
232          ENDDO
233        zplev( ngridmx,l ) = pp(1,jjp1,l)
234      ENDDO
235c
236c
237
238c   43. temperature naturelle (en K) et pressions milieux couches .
239c   ---------------------------------------------------------------
240
241      DO l=1,llm
242
243         pksurcp     =  ppk(1,1,l) / cpp
244         zplay(1,l)  =  preff * pksurcp ** unskap
245         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
246         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
247         ig0         = 2
248
249         DO j = 2, jjm
250            DO i = 1, iim
251              pksurcp        = ppk(i,j,l) / cpp
252              zplay(ig0,l)   = preff * pksurcp ** unskap
253              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
254              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
255              ig0            = ig0 + 1
256            ENDDO
257         ENDDO
258
259         pksurcp       = ppk(1,jjp1,l) / cpp
260         zplay(ig0,l)  = preff * pksurcp ** unskap
261         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
262         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
263
264      ENDDO
265
266c   43.bis humidite specifique (en kg/kg)
267c   -------------------------------------
268
269      DO iq=1,nqmx
270         DO l=1,llm
271            zqfi(1,l,iq) = pq(1,1,l,iq)
272            ig0          = 2
273            DO j=2,jjm
274               DO i = 1, iim
275                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
276                  ig0             = ig0 + 1
277               ENDDO
278            ENDDO
279            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
280         ENDDO
281      ENDDO
282
283c   convergence dynamique pour les traceurs "EAU"
284
285      DO iq=1,2
286         DO l=1,llm
287            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
288            ig0          = 2
289            DO j=2,jjm
290               DO i = 1, iim
291                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
292                  ig0             = ig0 + 1
293               ENDDO
294            ENDDO
295            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
296         ENDDO
297      ENDDO
298
299
300c   Geopotentiel calcule par rapport a la surface locale:
301c   -----------------------------------------------------
302
303      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
304      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
305      DO l=1,llm
306         DO ig=1,ngridmx
307           zphi(ig,l)=zphi(ig,l)-zphis(ig)
308         ENDDO
309      ENDDO
310
311c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
312c
313      DO l=1,llm
314        pvervel(1,l)=pw(1,1,l) * g /apoln
315        ig0=2
316       DO j=2,jjm
317           DO i = 1, iim
318              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
319              ig0 = ig0 + 1
320           ENDDO
321       ENDDO
322        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
323      ENDDO
324
325c      CALL initial0( ngridmx*llm,pvervel )
326c
327c
328c   45. champ u:
329c   ------------
330
331      DO 50 l=1,llm
332
333         DO 25 j=2,jjm
334            ig0 = 1+(j-2)*iim
335            zufi(ig0+1,l)= 0.5 *
336     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
337            pcvgu(ig0+1,l)= 0.5 *
338     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
339            DO 10 i=2,iim
340               zufi(ig0+i,l)= 0.5 *
341     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
342               pcvgu(ig0+i,l)= 0.5 *
343     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
34410         CONTINUE
34525      CONTINUE
346
34750    CONTINUE
348
349
350c   46.champ v:
351c   -----------
352
353      DO l=1,llm
354         DO j=2,jjm
355            ig0=1+(j-2)*iim
356            DO i=1,iim
357               zvfi(ig0+i,l)= 0.5 *
358     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
359               pcvgv(ig0+i,l)= 0.5 *
360     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
361            ENDDO
362         ENDDO
363      ENDDO
364
365
366c   47. champs de vents aux pole nord   
367c   ------------------------------
368c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
369c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
370
371      DO l=1,llm
372
373         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
374         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
375         DO i=2,iim
376            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
377            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
378         ENDDO
379
380         DO i=1,iim
381            zcos(i)   = COS(rlonv(i))*z1(i)
382            zcosbis(i)= COS(rlonv(i))*z1bis(i)
383            zsin(i)   = SIN(rlonv(i))*z1(i)
384            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
385         ENDDO
386
387         zufi(1,l)  = SSUM(iim,zcos,1)/pi
388         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
389         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
390         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
391
392      ENDDO
393
394
395c   48. champs de vents aux pole sud:
396c   ---------------------------------
397c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
398c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
399
400      DO l=1,llm
401
402         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
403         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
404         DO i=2,iim
405            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
406            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
407         ENDDO
408
409         DO i=1,iim
410            zcos(i)    = COS(rlonv(i))*z1(i)
411            zcosbis(i) = COS(rlonv(i))*z1bis(i)
412            zsin(i)    = SIN(rlonv(i))*z1(i)
413            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
414         ENDDO
415
416         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
417         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
418         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
419         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
420
421      ENDDO
422
423c-----------------------------------------------------------------------
424c   Appel de la physique:
425c   ---------------------
426
427
428      CALL physiq (ngridmx,llm,nq,debut,lafin,
429     ,   rdayvrai,rday_ecri,heure,dtphys,zplev,zplay,zphi,zphis,airefi,
430     ,             presnivs,clesphy0,  zufi, zvfi,ztfi, zqfi, 
431ccc     ,             pcvgu, pcvgv, pcvgt, pcvgq,
432     ,             pvervel,
433C - sorties
434     s             zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
435
436500   CONTINUE
437
438c-----------------------------------------------------------------------
439c   transformation des tendances physiques en tendances dynamiques:
440c   ---------------------------------------------------------------
441
442c  tendance sur la pression :
443c  -----------------------------------
444
445      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
446c
447ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
448
449c   62. enthalpie potentielle
450c   ---------------------
451
452      DO l=1,llm
453
454         DO i=1,iip1
455          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
456          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
457         ENDDO
458
459         DO j=2,jjm
460            ig0=1+(j-2)*iim
461            DO i=1,iim
462               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
463            ENDDO
464               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
465         ENDDO
466
467      ENDDO
468
469
470c   62. humidite specifique
471c   ---------------------
472
473      DO iq=1,nqmx
474         DO l=1,llm
475            DO i=1,iip1
476               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
477               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
478            ENDDO
479            DO j=2,jjm
480               ig0=1+(j-2)*iim
481               DO i=1,iim
482                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
483               ENDDO
484               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
485            ENDDO
486         ENDDO
487      ENDDO
488
489c   65. champ u:
490c   ------------
491
492      DO l=1,llm
493
494         DO i=1,iip1
495            pdufi(i,1,l)    = 0.
496            pdufi(i,jjp1,l) = 0.
497         ENDDO
498
499         DO j=2,jjm
500            ig0=1+(j-2)*iim
501            DO i=1,iim-1
502               pdufi(i,j,l)=
503     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
504            ENDDO
505            pdufi(iim,j,l)=
506     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
507            pdufi(iip1,j,l)=pdufi(1,j,l)
508         ENDDO
509
510      ENDDO
511
512
513c   67. champ v:
514c   ------------
515
516      DO l=1,llm
517
518         DO j=2,jjm-1
519            ig0=1+(j-2)*iim
520            DO i=1,iim
521               pdvfi(i,j,l)=
522     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
523            ENDDO
524            pdvfi(iip1,j,l) = pdvfi(1,j,l)
525         ENDDO
526      ENDDO
527
528
529c   68. champ v pres des poles:
530c   ---------------------------
531c      v = U * cos(long) + V * SIN(long)
532
533      DO l=1,llm
534
535         DO i=1,iim
536            pdvfi(i,1,l)=
537     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
538            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
539     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
540            pdvfi(i,1,l)=
541     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
542            pdvfi(i,jjm,l)=
543     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
544          ENDDO
545
546         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
547         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
548
549      ENDDO
550
551c-----------------------------------------------------------------------
552
553700   CONTINUE
554
555      firstcal = .FALSE.
556
557      RETURN
558      END
Note: See TracBrowser for help on using the repository browser.