source: LMDZ4/branches/LMDZ4_V2_patch/libf/dyn3d/calfis.F @ 4990

Last change on this file since 4990 was 665, checked in by (none), 19 years ago

This commit was manufactured by cvs2svn to create branch
'LMDZ4_V2_patch'.

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