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

Last change on this file since 18 was 18, checked in by lmdz, 25 years ago

correction erreur initialisation airefi L. Li

  • 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)/apoln
315        ig0=2
316       DO j=2,jjm
317           DO i = 1, iim
318              pvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
319              ig0 = ig0 + 1
320           ENDDO
321       ENDDO
322        pvervel(ig0,l)=pw(1,jjp1,l)/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, 
431     ,             pcvgu, pcvgv, pcvgt, pcvgq, pvervel,
432C - sorties
433     s             zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
434
435500   CONTINUE
436
437c-----------------------------------------------------------------------
438c   transformation des tendances physiques en tendances dynamiques:
439c   ---------------------------------------------------------------
440
441c  tendance sur la pression :
442c  -----------------------------------
443
444      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
445c
446ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
447
448c   62. enthalpie potentielle
449c   ---------------------
450
451      DO l=1,llm
452
453         DO i=1,iip1
454          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
455          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
456         ENDDO
457
458         DO j=2,jjm
459            ig0=1+(j-2)*iim
460            DO i=1,iim
461               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
462            ENDDO
463               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
464         ENDDO
465
466      ENDDO
467
468
469c   62. humidite specifique
470c   ---------------------
471
472      DO iq=1,nqmx
473         DO l=1,llm
474            DO i=1,iip1
475               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
476               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
477            ENDDO
478            DO j=2,jjm
479               ig0=1+(j-2)*iim
480               DO i=1,iim
481                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
482               ENDDO
483               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
484            ENDDO
485         ENDDO
486      ENDDO
487
488c   65. champ u:
489c   ------------
490
491      DO l=1,llm
492
493         DO i=1,iip1
494            pdufi(i,1,l)    = 0.
495            pdufi(i,jjp1,l) = 0.
496         ENDDO
497
498         DO j=2,jjm
499            ig0=1+(j-2)*iim
500            DO i=1,iim-1
501               pdufi(i,j,l)=
502     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
503            ENDDO
504            pdufi(iim,j,l)=
505     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
506            pdufi(iip1,j,l)=pdufi(1,j,l)
507         ENDDO
508
509      ENDDO
510
511
512c   67. champ v:
513c   ------------
514
515      DO l=1,llm
516
517         DO j=2,jjm-1
518            ig0=1+(j-2)*iim
519            DO i=1,iim
520               pdvfi(i,j,l)=
521     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
522            ENDDO
523            pdvfi(iip1,j,l) = pdvfi(1,j,l)
524         ENDDO
525      ENDDO
526
527
528c   68. champ v pres des poles:
529c   ---------------------------
530c      v = U * cos(long) + V * SIN(long)
531
532      DO l=1,llm
533
534         DO i=1,iim
535            pdvfi(i,1,l)=
536     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
537            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
538     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
539            pdvfi(i,1,l)=
540     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
541            pdvfi(i,jjm,l)=
542     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
543          ENDDO
544
545         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
546         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
547
548      ENDDO
549
550c-----------------------------------------------------------------------
551
552700   CONTINUE
553
554      firstcal = .FALSE.
555
556      RETURN
557      END
Note: See TracBrowser for help on using the repository browser.