source: LMDZ4/trunk/libf/dyn3dpar/calfis.F @ 757

Last change on this file since 757 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.2 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_CH4
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
164
165#ifdef INCA_CH4
166      REAL flxw(iip1,jjp1,llm)
167      REAL flxwfi(ngridmx,llm)
168#endif
169c
170     
171      REAL SSUM
172
173      LOGICAL firstcal, debut
174      DATA firstcal/.true./
175      SAVE firstcal,debut
176      REAL rdayvrai
177c
178c-----------------------------------------------------------------------
179c
180c    1. Initialisations :
181c    --------------------
182c
183
184      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
185         PRINT*,'STOP dans calfis'
186         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
187         PRINT*,'  ngridmx  jjm   iim   '
188         PRINT*,ngridmx,jjm,iim
189         STOP
190      ENDIF
191
192c-----------------------------------------------------------------------
193c   latitude, longitude et aires des mailles pour la physique:
194c   ----------------------------------------------------------
195
196c
197      IF ( firstcal )  THEN
198          debut = .TRUE.
199      ELSE
200          debut = .FALSE.
201      ENDIF
202
203c
204c
205c-----------------------------------------------------------------------
206c   40. transformation des variables dynamiques en variables physiques:
207c   ---------------------------------------------------------------
208
209c   41. pressions au sol (en Pascals)
210c   ----------------------------------
211
212       
213      zpsrf(1) = pps(1,1)
214
215      ig0  = 2
216      DO j = 2,jjm
217         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
218         ig0 = ig0+iim
219      ENDDO
220
221      zpsrf(ngridmx) = pps(1,jjp1)
222
223
224c   42. pression intercouches :
225c
226c   -----------------------------------------------------------------
227c     .... zplev  definis aux (llm +1) interfaces des couches  ....
228c     .... zplay  definis aux (  llm )    milieux des couches  ....
229c   -----------------------------------------------------------------
230
231c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
232c
233       unskap   = 1./ kappa
234c
235      DO l = 1, llmp1
236        zplev( 1,l ) = pp(1,1,l)
237        ig0 = 2
238          DO j = 2, jjm
239             DO i =1, iim
240              zplev( ig0,l ) = pp(i,j,l)
241              ig0 = ig0 +1
242             ENDDO
243          ENDDO
244        zplev( ngridmx,l ) = pp(1,jjp1,l)
245      ENDDO
246c
247c
248
249c   43. temperature naturelle (en K) et pressions milieux couches .
250c   ---------------------------------------------------------------
251
252      DO l=1,llm
253
254         pksurcp     =  ppk(1,1,l) / cpp
255         zplay(1,l)  =  preff * pksurcp ** unskap
256         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
257         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
258         ig0         = 2
259
260         DO j = 2, jjm
261            DO i = 1, iim
262              pksurcp        = ppk(i,j,l) / cpp
263              zplay(ig0,l)   = preff * pksurcp ** unskap
264              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
265              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
266              ig0            = ig0 + 1
267            ENDDO
268         ENDDO
269
270         pksurcp       = ppk(1,jjp1,l) / cpp
271         zplay(ig0,l)  = preff * pksurcp ** unskap
272         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
273         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
274
275      ENDDO
276
277c   43.bis traceurs
278c   ---------------
279c
280      DO iq=1,nq
281          iiq=niadv(iq)
282         DO l=1,llm
283            zqfi(1,l,iq) = pq(1,1,l,iiq)
284            ig0          = 2
285            DO j=2,jjm
286               DO i = 1, iim
287                  zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
288                  ig0             = ig0 + 1
289               ENDDO
290            ENDDO
291            zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
292         ENDDO
293      ENDDO
294
295c   convergence dynamique pour les traceurs "EAU"
296
297      DO iq=1,2
298         DO l=1,llm
299            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
300            ig0          = 2
301            DO j=2,jjm
302               DO i = 1, iim
303                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
304                  ig0             = ig0 + 1
305               ENDDO
306            ENDDO
307            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
308         ENDDO
309      ENDDO
310
311
312c   Geopotentiel calcule par rapport a la surface locale:
313c   -----------------------------------------------------
314
315      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
316      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
317      DO l=1,llm
318         DO ig=1,ngridmx
319           zphi(ig,l)=zphi(ig,l)-zphis(ig)
320         ENDDO
321      ENDDO
322
323c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
324c
325      DO l=1,llm
326        pvervel(1,l)=pw(1,1,l) * g /apoln
327        ig0=2
328       DO j=2,jjm
329           DO i = 1, iim
330              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
331              ig0 = ig0 + 1
332           ENDDO
333       ENDDO
334        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
335      ENDDO
336
337c
338c   45. champ u:
339c   ------------
340
341      DO 50 l=1,llm
342
343         DO 25 j=2,jjm
344            ig0 = 1+(j-2)*iim
345            zufi(ig0+1,l)= 0.5 *
346     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
347            pcvgu(ig0+1,l)= 0.5 *
348     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
349            DO 10 i=2,iim
350               zufi(ig0+i,l)= 0.5 *
351     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
352               pcvgu(ig0+i,l)= 0.5 *
353     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
35410         CONTINUE
35525      CONTINUE
356
35750    CONTINUE
358
359
360c   46.champ v:
361c   -----------
362
363      DO l=1,llm
364         DO j=2,jjm
365            ig0=1+(j-2)*iim
366            DO i=1,iim
367               zvfi(ig0+i,l)= 0.5 *
368     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
369               pcvgv(ig0+i,l)= 0.5 *
370     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
371            ENDDO
372         ENDDO
373      ENDDO
374
375
376c   47. champs de vents aux pole nord   
377c   ------------------------------
378c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
379c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
380
381      DO l=1,llm
382
383         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
384         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
385         DO i=2,iim
386            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
387            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
388         ENDDO
389
390         DO i=1,iim
391            zcos(i)   = COS(rlonv(i))*z1(i)
392            zcosbis(i)= COS(rlonv(i))*z1bis(i)
393            zsin(i)   = SIN(rlonv(i))*z1(i)
394            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
395         ENDDO
396
397         zufi(1,l)  = SSUM(iim,zcos,1)/pi
398         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
399         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
400         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
401
402      ENDDO
403
404
405c   48. champs de vents aux pole sud:
406c   ---------------------------------
407c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
408c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
409
410      DO l=1,llm
411
412         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
413         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
414         DO i=2,iim
415            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
416            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
417         ENDDO
418
419         DO i=1,iim
420            zcos(i)    = COS(rlonv(i))*z1(i)
421            zcosbis(i) = COS(rlonv(i))*z1bis(i)
422            zsin(i)    = SIN(rlonv(i))*z1(i)
423            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
424         ENDDO
425
426         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
427         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
428         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
429         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
430
431      ENDDO
432
433
434#ifdef INCA_CH4
435      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
436#endif
437
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     .             pvervel,
463#ifdef INCA_CH4
464     .             flxwfi,
465#endif
466     .             zdufi,
467     .             zdvfi,
468     .             zdtfi,
469     .             zdqfi,
470     .             zdpsrf)
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.