source: LMDZ4/trunk/libf/dyn3d/calfis.F @ 1117

Last change on this file since 1117 was 960, checked in by lsce, 16 years ago
  • Ajoute du parametre config_inca dans conf_gcm.F config_inca='none'(sans INCA, par defaut) config_inca='chem'(avec INCA config chemie) config_inca='aero'(avec INCA config aerosol)
  • Menage parmis les cles CPP INCA
  • Enleve le calcul d'omega dans calfis.F et active le calcul correspondant dans advtrac.F(avant uniquement pour INCA).

JG

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