source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F @ 265

Last change on this file since 265 was 257, checked in by lmdzadmin, 23 years ago

cosmetique
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.7 KB
Line 
1C
2C $Header$
3C
4      SUBROUTINE calfis(nq, lafin, rdayvrai,rday_ecri, heure,
5     $            pucov,pvcov,pteta,pq,pmasse,pps,pp,ppk,pphis,pphi,
6     $            pducov,pdvcov,pdteta,pdq,pw, clesphy0,
7     $            pdufi,pdvfi,pdhfi,pdqfi,pdpsfi                   )
8c
9c    Auteur :  P. Le Van, F. Hourdin
10c   .........
11
12      IMPLICIT NONE
13c=======================================================================
14c
15c   1. rearrangement des tableaux et transformation
16c      variables dynamiques  >  variables physiques
17c   2. calcul des termes physiques
18c   3. retransformation des tendances physiques en tendances dynamiques
19c
20c   remarques:
21c   ----------
22c
23c    - les vents sont donnes dans la physique par leurs composantes
24c      naturelles.
25c    - la variable thermodynamique de la physique est une variable
26c      intensive :   T
27c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
28c    - les deux seules variables dependant de la geometrie necessaires
29c      pour la physique sont la latitude pour le rayonnement et
30c      l'aire de la maille quand on veut integrer une grandeur
31c      horizontalement.
32c    - les points de la physique sont les points scalaires de la
33c      la dynamique; numerotation:
34c          1 pour le pole nord
35c          (jjm-1)*iim pour l'interieur du domaine
36c          ngridmx pour le pole sud
37c      ---> ngridmx=2+(jjm-1)*iim
38c
39c     Input :
40c     -------
41c       ecritphy        frequence d'ecriture (en jours)de histphy
42c       pucov           covariant zonal velocity
43c       pvcov           covariant meridional velocity
44c       pteta           potential temperature
45c       pps             surface pressure
46c       pmasse          masse d'air dans chaque maille
47c       pts             surface temperature  (K)
48c       callrad         clef d'appel au rayonnement
49c
50c    Output :
51c    --------
52c        pdufi          tendency for the natural zonal velocity (ms-1)
53c        pdvfi          tendency for the natural meridional velocity
54c        pdhfi          tendency for the potential temperature
55c        pdtsfi         tendency for the surface temperature
56c
57c        pdtrad         radiative tendencies  \  both input
58c        pfluxrad       radiative fluxes      /  and output
59c
60c=======================================================================
61c
62c-----------------------------------------------------------------------
63c
64c    0.  Declarations :
65c    ------------------
66
67#include "dimensions.h"
68#include "paramet.h"
69#include "temps.h"
70
71      INTEGER ngridmx,nq
72      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
73
74#include "comconst.h"
75#include "comvert.h"
76#include "comgeom2.h"
77#include "control.h"
78
79c    Arguments :
80c    -----------
81      LOGICAL  lafin
82      REAL heure
83
84      REAL pvcov(iip1,jjm,llm)
85      REAL pucov(iip1,jjp1,llm)
86      REAL pteta(iip1,jjp1,llm)
87      REAL pmasse(iip1,jjp1,llm)
88      REAL pq(iip1,jjp1,llm,nqmx)
89      REAL pphis(iip1,jjp1)
90      REAL pphi(iip1,jjp1,llm)
91c
92      REAL pdvcov(iip1,jjm,llm)
93      REAL pducov(iip1,jjp1,llm)
94      REAL pdteta(iip1,jjp1,llm)
95      REAL pdq(iip1,jjp1,llm,nqmx)
96c
97      REAL pw(iip1,jjp1,llm)
98c
99      REAL pps(iip1,jjp1)
100      REAL pp(iip1,jjp1,llmp1)
101      REAL ppk(iip1,jjp1,llm)
102c
103      REAL pdvfi(iip1,jjm,llm)
104      REAL pdufi(iip1,jjp1,llm)
105      REAL pdhfi(iip1,jjp1,llm)
106      REAL pdqfi(iip1,jjp1,llm,nqmx)
107      REAL pdpsfi(iip1,jjp1)
108
109      INTEGER        longcles
110      PARAMETER    ( longcles = 20 )
111      REAL clesphy0( longcles )
112
113
114c    Local variables :
115c    -----------------
116
117      INTEGER i,j,l,ig0,ig,iq
118      REAL zpsrf(ngridmx)
119      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
120      REAL zphi(ngridmx,llm),zphis(ngridmx)
121c
122      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
123      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
124c
125      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
126      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
127c
128      REAL pvervel(ngridmx,llm)
129c
130      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
131      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
132      REAL zdpsrf(ngridmx)
133      REAL zcufi(ngridmx),zcvfi(ngridmx)
134c
135      REAL zsin(iim),zcos(iim),z1(iim)
136      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
137      REAL unskap, pksurcp
138c
139     
140      EXTERNAL gr_dyn_fi,gr_fi_dyn
141      EXTERNAL physiq,multipl
142      REAL SSUM
143      EXTERNAL SSUM
144
145      REAL latfi(ngridmx),lonfi(ngridmx)
146      REAL airefi(ngridmx)
147      SAVE latfi, lonfi, airefi
148
149      LOGICAL firstcal, debut
150      DATA firstcal/.true./
151      SAVE firstcal,debut
152      REAL rdayvrai,rday_ecri
153c
154c-----------------------------------------------------------------------
155c
156c    1. Initialisations :
157c    --------------------
158c
159
160      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
161         PRINT*,'STOP dans calfis'
162         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
163         PRINT*,'  ngridmx  jjm   iim   '
164         PRINT*,ngridmx,jjm,iim
165         STOP
166      ENDIF
167
168c-----------------------------------------------------------------------
169c   latitude, longitude et aires des mailles pour la physique:
170c   ----------------------------------------------------------
171
172c
173      IF ( firstcal )  THEN
174          debut = .TRUE.
175      ELSE
176          debut = .FALSE.
177      ENDIF
178
179c
180      IF (firstcal) THEN
181         latfi(1)=rlatu(1)
182         lonfi(1)=0.
183         zcufi(1) = cu(1,1)
184         zcvfi(1) = cv(1,1)
185         DO j=2,jjm
186            DO i=1,iim
187               latfi((j-2)*iim+1+i)= rlatu(j)
188               lonfi((j-2)*iim+1+i)= rlonv(i)
189               zcufi((j-2)*iim+1+i) = cu(i,j)
190               zcvfi((j-2)*iim+1+i) = cv(i,j)
191            ENDDO
192         ENDDO
193         latfi(ngridmx)= rlatu(jjp1)
194         lonfi(ngridmx)= 0.
195         zcufi(ngridmx) = cu(1,jjp1)
196         zcvfi(ngridmx) = cv(1,jjm)
197         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
198          PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
199         CALL inifis(ngridmx,llm,daysec,day_ini,dtphys ,
200     ,                latfi,lonfi,airefi,rad,g,r,cpp     )
201          ENDIF
202
203c
204c-----------------------------------------------------------------------
205c   40. transformation des variables dynamiques en variables physiques:
206c   ---------------------------------------------------------------
207
208c   41. pressions au sol (en Pascals)
209c   ----------------------------------
210
211       
212      zpsrf(1) = pps(1,1)
213
214      ig0  = 2
215      DO j = 2,jjm
216         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
217         ig0 = ig0+iim
218      ENDDO
219
220      zpsrf(ngridmx) = pps(1,jjp1)
221
222
223c   42. pression intercouches :
224c
225c   -----------------------------------------------------------------
226c     .... zplev  definis aux (llm +1) interfaces des couches  ....
227c     .... zplay  definis aux (  llm )    milieux des couches  ....
228c   -----------------------------------------------------------------
229
230c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
231c
232       unskap   = 1./ kappa
233c
234      DO l = 1, llmp1
235        zplev( 1,l ) = pp(1,1,l)
236        ig0 = 2
237          DO j = 2, jjm
238             DO i =1, iim
239              zplev( ig0,l ) = pp(i,j,l)
240              ig0 = ig0 +1
241             ENDDO
242          ENDDO
243        zplev( ngridmx,l ) = pp(1,jjp1,l)
244      ENDDO
245c
246c
247
248c   43. temperature naturelle (en K) et pressions milieux couches .
249c   ---------------------------------------------------------------
250
251      DO l=1,llm
252
253         pksurcp     =  ppk(1,1,l) / cpp
254         zplay(1,l)  =  preff * pksurcp ** unskap
255         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
256         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
257         ig0         = 2
258
259         DO j = 2, jjm
260            DO i = 1, iim
261              pksurcp        = ppk(i,j,l) / cpp
262              zplay(ig0,l)   = preff * pksurcp ** unskap
263              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
264              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
265              ig0            = ig0 + 1
266            ENDDO
267         ENDDO
268
269         pksurcp       = ppk(1,jjp1,l) / cpp
270         zplay(ig0,l)  = preff * pksurcp ** unskap
271         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
272         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
273
274      ENDDO
275
276c   43.bis humidite specifique (en kg/kg)
277c   -------------------------------------
278
279      DO iq=1,nqmx
280         DO l=1,llm
281            zqfi(1,l,iq) = pq(1,1,l,iq)
282            ig0          = 2
283            DO j=2,jjm
284               DO i = 1, iim
285                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
286                  ig0             = ig0 + 1
287               ENDDO
288            ENDDO
289            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
290         ENDDO
291      ENDDO
292
293c   convergence dynamique pour les traceurs "EAU"
294
295      DO iq=1,2
296         DO l=1,llm
297            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
298            ig0          = 2
299            DO j=2,jjm
300               DO i = 1, iim
301                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
302                  ig0             = ig0 + 1
303               ENDDO
304            ENDDO
305            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
306         ENDDO
307      ENDDO
308
309
310c   Geopotentiel calcule par rapport a la surface locale:
311c   -----------------------------------------------------
312
313      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
314      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
315      DO l=1,llm
316         DO ig=1,ngridmx
317           zphi(ig,l)=zphi(ig,l)-zphis(ig)
318         ENDDO
319      ENDDO
320
321c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
322c
323      DO l=1,llm
324        pvervel(1,l)=pw(1,1,l) * g /apoln
325        ig0=2
326       DO j=2,jjm
327           DO i = 1, iim
328              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
329              ig0 = ig0 + 1
330           ENDDO
331       ENDDO
332        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
333      ENDDO
334
335c      CALL initial0( ngridmx*llm,pvervel )
336c
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
433c-----------------------------------------------------------------------
434c   Appel de la physique:
435c   ---------------------
436
437
438      CALL physiq (ngridmx, llm, nq,
439     ,    debut, lafin, rdayvrai, rday_ecri, heure, dtphys,
440     ,    zplev, zplay, zphi, zphis, airefi, presnivs, clesphy0,
441     ,    zufi, zvfi, ztfi, zqfi, 
442     ,    pvervel, zcufi, zcvfi,
443C - sorties
444     s    zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
445
446500   CONTINUE
447
448c-----------------------------------------------------------------------
449c   transformation des tendances physiques en tendances dynamiques:
450c   ---------------------------------------------------------------
451
452c  tendance sur la pression :
453c  -----------------------------------
454
455      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
456c
457ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
458
459c   62. enthalpie potentielle
460c   ---------------------
461
462      DO l=1,llm
463
464         DO i=1,iip1
465          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
466          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
467         ENDDO
468
469         DO j=2,jjm
470            ig0=1+(j-2)*iim
471            DO i=1,iim
472               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
473            ENDDO
474               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
475         ENDDO
476
477      ENDDO
478
479
480c   62. humidite specifique
481c   ---------------------
482
483      DO iq=1,nqmx
484         DO l=1,llm
485            DO i=1,iip1
486               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
487               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
488            ENDDO
489            DO j=2,jjm
490               ig0=1+(j-2)*iim
491               DO i=1,iim
492                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
493               ENDDO
494               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
495            ENDDO
496         ENDDO
497      ENDDO
498
499c   65. champ u:
500c   ------------
501
502      DO l=1,llm
503
504         DO i=1,iip1
505            pdufi(i,1,l)    = 0.
506            pdufi(i,jjp1,l) = 0.
507         ENDDO
508
509         DO j=2,jjm
510            ig0=1+(j-2)*iim
511            DO i=1,iim-1
512               pdufi(i,j,l)=
513     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
514            ENDDO
515            pdufi(iim,j,l)=
516     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
517            pdufi(iip1,j,l)=pdufi(1,j,l)
518         ENDDO
519
520      ENDDO
521
522
523c   67. champ v:
524c   ------------
525
526      DO l=1,llm
527
528         DO j=2,jjm-1
529            ig0=1+(j-2)*iim
530            DO i=1,iim
531               pdvfi(i,j,l)=
532     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
533            ENDDO
534            pdvfi(iip1,j,l) = pdvfi(1,j,l)
535         ENDDO
536      ENDDO
537
538
539c   68. champ v pres des poles:
540c   ---------------------------
541c      v = U * cos(long) + V * SIN(long)
542
543      DO l=1,llm
544
545         DO i=1,iim
546            pdvfi(i,1,l)=
547     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
548            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
549     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
550            pdvfi(i,1,l)=
551     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
552            pdvfi(i,jjm,l)=
553     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
554          ENDDO
555
556         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
557         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
558
559      ENDDO
560
561c-----------------------------------------------------------------------
562
563700   CONTINUE
564
565      firstcal = .FALSE.
566
567      RETURN
568      END
Note: See TracBrowser for help on using the repository browser.