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

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

Merge par rapport a la branche principale
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 
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)
130      REAL zcufi(ngridmx),zcvfi(ngridmx)
131c
132      REAL zsin(iim),zcos(iim),z1(iim)
133      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
134      REAL unskap, pksurcp
135c
136     
137      EXTERNAL gr_dyn_fi,gr_fi_dyn
138      EXTERNAL physiq,multipl
139      REAL SSUM
140      EXTERNAL SSUM
141
142      REAL latfi(ngridmx),lonfi(ngridmx)
143      REAL airefi(ngridmx)
144      SAVE latfi, lonfi, airefi
145
146      LOGICAL firstcal, debut
147      DATA firstcal/.true./
148      SAVE firstcal,debut
149      REAL rdayvrai,rday_ecri
150c
151c-----------------------------------------------------------------------
152c
153c    1. Initialisations :
154c    --------------------
155c
156
157      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
158         PRINT*,'STOP dans calfis'
159         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
160         PRINT*,'  ngridmx  jjm   iim   '
161         PRINT*,ngridmx,jjm,iim
162         STOP
163      ENDIF
164
165c-----------------------------------------------------------------------
166c   latitude, longitude et aires des mailles pour la physique:
167c   ----------------------------------------------------------
168
169c
170      IF ( firstcal )  THEN
171          debut = .TRUE.
172      ELSE
173          debut = .FALSE.
174      ENDIF
175
176c
177      IF (firstcal) THEN
178         latfi(1)=rlatu(1)
179         lonfi(1)=0.
180         zcufi(1) = cu(1,1)
181         zcvfi(1) = cv(1,1)
182         DO j=2,jjm
183            DO i=1,iim
184               latfi((j-2)*iim+1+i)= rlatu(j)
185               lonfi((j-2)*iim+1+i)= rlonv(i)
186               zcufi((j-2)*iim+1+i) = cu(i,j)
187               zcvfi((j-2)*iim+1+i) = cv(i,j)
188            ENDDO
189         ENDDO
190         latfi(ngridmx)= rlatu(jjp1)
191         lonfi(ngridmx)= 0.
192         zcufi(ngridmx) = cu(1,jjp1)
193         zcvfi(ngridmx) = cv(1,jjm)
194         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
195          PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
196         CALL inifis(ngridmx,llm,daysec,day_ini,dtphys ,
197     ,                latfi,lonfi,airefi,rad,g,r,cpp     )
198          ENDIF
199
200c
201c-----------------------------------------------------------------------
202c   40. transformation des variables dynamiques en variables physiques:
203c   ---------------------------------------------------------------
204
205c   41. pressions au sol (en Pascals)
206c   ----------------------------------
207
208       
209      zpsrf(1) = pps(1,1)
210
211      ig0  = 2
212      DO j = 2,jjm
213         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
214         ig0 = ig0+iim
215      ENDDO
216
217      zpsrf(ngridmx) = pps(1,jjp1)
218
219
220c   42. pression intercouches :
221c
222c   -----------------------------------------------------------------
223c     .... zplev  definis aux (llm +1) interfaces des couches  ....
224c     .... zplay  definis aux (  llm )    milieux des couches  ....
225c   -----------------------------------------------------------------
226
227c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
228c
229       unskap   = 1./ kappa
230c
231      DO l = 1, llmp1
232        zplev( 1,l ) = pp(1,1,l)
233        ig0 = 2
234          DO j = 2, jjm
235             DO i =1, iim
236              zplev( ig0,l ) = pp(i,j,l)
237              ig0 = ig0 +1
238             ENDDO
239          ENDDO
240        zplev( ngridmx,l ) = pp(1,jjp1,l)
241      ENDDO
242c
243c
244
245c   43. temperature naturelle (en K) et pressions milieux couches .
246c   ---------------------------------------------------------------
247
248      DO l=1,llm
249
250         pksurcp     =  ppk(1,1,l) / cpp
251         zplay(1,l)  =  preff * pksurcp ** unskap
252         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
253         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
254         ig0         = 2
255
256         DO j = 2, jjm
257            DO i = 1, iim
258              pksurcp        = ppk(i,j,l) / cpp
259              zplay(ig0,l)   = preff * pksurcp ** unskap
260              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
261              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
262              ig0            = ig0 + 1
263            ENDDO
264         ENDDO
265
266         pksurcp       = ppk(1,jjp1,l) / cpp
267         zplay(ig0,l)  = preff * pksurcp ** unskap
268         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
269         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
270
271      ENDDO
272
273c   43.bis humidite specifique (en kg/kg)
274c   -------------------------------------
275
276      DO iq=1,nqmx
277         DO l=1,llm
278            zqfi(1,l,iq) = pq(1,1,l,iq)
279            ig0          = 2
280            DO j=2,jjm
281               DO i = 1, iim
282                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
283                  ig0             = ig0 + 1
284               ENDDO
285            ENDDO
286            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
287         ENDDO
288      ENDDO
289
290c   convergence dynamique pour les traceurs "EAU"
291
292      DO iq=1,2
293         DO l=1,llm
294            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
295            ig0          = 2
296            DO j=2,jjm
297               DO i = 1, iim
298                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
299                  ig0             = ig0 + 1
300               ENDDO
301            ENDDO
302            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
303         ENDDO
304      ENDDO
305
306
307c   Geopotentiel calcule par rapport a la surface locale:
308c   -----------------------------------------------------
309
310      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
311      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
312      DO l=1,llm
313         DO ig=1,ngridmx
314           zphi(ig,l)=zphi(ig,l)-zphis(ig)
315         ENDDO
316      ENDDO
317
318c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
319c
320      DO l=1,llm
321        pvervel(1,l)=pw(1,1,l) * g /apoln
322        ig0=2
323       DO j=2,jjm
324           DO i = 1, iim
325              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
326              ig0 = ig0 + 1
327           ENDDO
328       ENDDO
329        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
330      ENDDO
331
332c      CALL initial0( ngridmx*llm,pvervel )
333c
334c
335c   45. champ u:
336c   ------------
337
338      DO 50 l=1,llm
339
340         DO 25 j=2,jjm
341            ig0 = 1+(j-2)*iim
342            zufi(ig0+1,l)= 0.5 *
343     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
344            pcvgu(ig0+1,l)= 0.5 *
345     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
346            DO 10 i=2,iim
347               zufi(ig0+i,l)= 0.5 *
348     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
349               pcvgu(ig0+i,l)= 0.5 *
350     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
35110         CONTINUE
35225      CONTINUE
353
35450    CONTINUE
355
356
357c   46.champ v:
358c   -----------
359
360      DO l=1,llm
361         DO j=2,jjm
362            ig0=1+(j-2)*iim
363            DO i=1,iim
364               zvfi(ig0+i,l)= 0.5 *
365     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
366               pcvgv(ig0+i,l)= 0.5 *
367     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
368            ENDDO
369         ENDDO
370      ENDDO
371
372
373c   47. champs de vents aux pole nord   
374c   ------------------------------
375c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
376c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
377
378      DO l=1,llm
379
380         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
381         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
382         DO i=2,iim
383            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
384            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
385         ENDDO
386
387         DO i=1,iim
388            zcos(i)   = COS(rlonv(i))*z1(i)
389            zcosbis(i)= COS(rlonv(i))*z1bis(i)
390            zsin(i)   = SIN(rlonv(i))*z1(i)
391            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
392         ENDDO
393
394         zufi(1,l)  = SSUM(iim,zcos,1)/pi
395         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
396         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
397         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
398
399      ENDDO
400
401
402c   48. champs de vents aux pole sud:
403c   ---------------------------------
404c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
405c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
406
407      DO l=1,llm
408
409         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
410         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
411         DO i=2,iim
412            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
413            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
414         ENDDO
415
416         DO i=1,iim
417            zcos(i)    = COS(rlonv(i))*z1(i)
418            zcosbis(i) = COS(rlonv(i))*z1bis(i)
419            zsin(i)    = SIN(rlonv(i))*z1(i)
420            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
421         ENDDO
422
423         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
424         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
425         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
426         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
427
428      ENDDO
429
430c-----------------------------------------------------------------------
431c   Appel de la physique:
432c   ---------------------
433
434
435      CALL physiq (ngridmx,llm,nq,debut,lafin,
436     ,   rdayvrai,rday_ecri,heure,dtphys,zplev,zplay,zphi,zphis,airefi,
437     ,             presnivs,clesphy0,  zufi, zvfi,ztfi, zqfi, 
438ccc     ,             pcvgu, pcvgv, pcvgt, pcvgq,
439     ,             pvervel, zcufi, zcvfi,
440C - sorties
441     s             zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
442
443500   CONTINUE
444
445c-----------------------------------------------------------------------
446c   transformation des tendances physiques en tendances dynamiques:
447c   ---------------------------------------------------------------
448
449c  tendance sur la pression :
450c  -----------------------------------
451
452      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
453c
454ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
455
456c   62. enthalpie potentielle
457c   ---------------------
458
459      DO l=1,llm
460
461         DO i=1,iip1
462          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
463          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
464         ENDDO
465
466         DO j=2,jjm
467            ig0=1+(j-2)*iim
468            DO i=1,iim
469               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
470            ENDDO
471               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
472         ENDDO
473
474      ENDDO
475
476
477c   62. humidite specifique
478c   ---------------------
479
480      DO iq=1,nqmx
481         DO l=1,llm
482            DO i=1,iip1
483               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
484               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
485            ENDDO
486            DO j=2,jjm
487               ig0=1+(j-2)*iim
488               DO i=1,iim
489                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
490               ENDDO
491               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
492            ENDDO
493         ENDDO
494      ENDDO
495
496c   65. champ u:
497c   ------------
498
499      DO l=1,llm
500
501         DO i=1,iip1
502            pdufi(i,1,l)    = 0.
503            pdufi(i,jjp1,l) = 0.
504         ENDDO
505
506         DO j=2,jjm
507            ig0=1+(j-2)*iim
508            DO i=1,iim-1
509               pdufi(i,j,l)=
510     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
511            ENDDO
512            pdufi(iim,j,l)=
513     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
514            pdufi(iip1,j,l)=pdufi(1,j,l)
515         ENDDO
516
517      ENDDO
518
519
520c   67. champ v:
521c   ------------
522
523      DO l=1,llm
524
525         DO j=2,jjm-1
526            ig0=1+(j-2)*iim
527            DO i=1,iim
528               pdvfi(i,j,l)=
529     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
530            ENDDO
531            pdvfi(iip1,j,l) = pdvfi(1,j,l)
532         ENDDO
533      ENDDO
534
535
536c   68. champ v pres des poles:
537c   ---------------------------
538c      v = U * cos(long) + V * SIN(long)
539
540      DO l=1,llm
541
542         DO i=1,iim
543            pdvfi(i,1,l)=
544     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
545            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
546     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
547            pdvfi(i,1,l)=
548     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
549            pdvfi(i,jjm,l)=
550     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
551          ENDDO
552
553         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
554         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
555
556      ENDDO
557
558c-----------------------------------------------------------------------
559
560700   CONTINUE
561
562      firstcal = .FALSE.
563
564      RETURN
565      END
Note: See TracBrowser for help on using the repository browser.