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

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

Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA
LF

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