source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3d/calfis.F @ 1440

Last change on this file since 1440 was 1225, checked in by lguez, 15 years ago

Added some "intent" attributes in declarations.

In "phyredem", "dtime" is not declared. It is not in any included
file. Probably accepted by compilers as an intrinsic non-standard
function. Removed this element of "tab_cntrl".

Added some "only" clauses in "use" statements.

If the ozone field is read from a file, it is now updated every
360th of the length of the current year, regardless of that length.

In "physiq", "omega" was output before it was defined. Moved the
output instruction after the definition.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 KB
Line 
1!
2! $Id: calfis.F 1225 2009-08-07 16:29:23Z acozic $
3!
4C
5C
6      SUBROUTINE calfis(lafin,
7     $                  jD_cur, jH_cur,
8     $                  pucov,
9     $                  pvcov,
10     $                  pteta,
11     $                  pq,
12     $                  pmasse,
13     $                  pps,
14     $                  pp,
15     $                  ppk,
16     $                  pphis,
17     $                  pphi,
18     $                  pducov,
19     $                  pdvcov,
20     $                  pdteta,
21     $                  pdq,
22     $                  flxw,
23     $                  clesphy0,
24     $                  pdufi,
25     $                  pdvfi,
26     $                  pdhfi,
27     $                  pdqfi,
28     $                  pdpsfi)
29c
30c    Auteur :  P. Le Van, F. Hourdin
31c   .........
32      USE infotrac
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       pucov           covariant zonal velocity
64c       pvcov           covariant meridional velocity
65c       pteta           potential temperature
66c       pps             surface pressure
67c       pmasse          masse d'air dans chaque maille
68c       pts             surface temperature  (K)
69c       callrad         clef d'appel au rayonnement
70c
71c    Output :
72c    --------
73c        pdufi          tendency for the natural zonal velocity (ms-1)
74c        pdvfi          tendency for the natural meridional velocity
75c        pdhfi          tendency for the potential temperature
76c        pdtsfi         tendency for the surface temperature
77c
78c        pdtrad         radiative tendencies  \  both input
79c        pfluxrad       radiative fluxes      /  and output
80c
81c=======================================================================
82c
83c-----------------------------------------------------------------------
84c
85c    0.  Declarations :
86c    ------------------
87
88#include "dimensions.h"
89#include "paramet.h"
90#include "temps.h"
91
92      INTEGER ngridmx
93      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
94
95#include "comconst.h"
96#include "comvert.h"
97#include "comgeom2.h"
98#include "control.h"
99
100c    Arguments :
101c    -----------
102      LOGICAL  lafin
103
104
105      REAL pvcov(iip1,jjm,llm)
106      REAL pucov(iip1,jjp1,llm)
107      REAL pteta(iip1,jjp1,llm)
108      REAL pmasse(iip1,jjp1,llm)
109      REAL pq(iip1,jjp1,llm,nqtot)
110      REAL pphis(iip1,jjp1)
111      REAL pphi(iip1,jjp1,llm)
112c
113      REAL pdvcov(iip1,jjm,llm)
114      REAL pducov(iip1,jjp1,llm)
115      REAL pdteta(iip1,jjp1,llm)
116      REAL pdq(iip1,jjp1,llm,nqtot)
117c
118      REAL pps(iip1,jjp1)
119      REAL pp(iip1,jjp1,llmp1)
120      REAL ppk(iip1,jjp1,llm)
121c
122      REAL pdvfi(iip1,jjm,llm)
123      REAL pdufi(iip1,jjp1,llm)
124      REAL pdhfi(iip1,jjp1,llm)
125      REAL pdqfi(iip1,jjp1,llm,nqtot)
126      REAL pdpsfi(iip1,jjp1)
127
128      INTEGER        longcles
129      PARAMETER    ( longcles = 20 )
130      REAL clesphy0( longcles )
131
132
133c    Local variables :
134c    -----------------
135
136      INTEGER i,j,l,ig0,ig,iq,iiq
137      REAL zpsrf(ngridmx)
138      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
139      REAL zphi(ngridmx,llm),zphis(ngridmx)
140c
141      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
142      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
143c
144      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
145      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
146c
147      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
148      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
149      REAL zdpsrf(ngridmx)
150c
151      REAL zsin(iim),zcos(iim),z1(iim)
152      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
153      REAL unskap, pksurcp
154c
155cIM diagnostique PVteta, Amip2
156      INTEGER ntetaSTD
157      PARAMETER(ntetaSTD=3)
158      REAL rtetaSTD(ntetaSTD)
159      DATA rtetaSTD/350., 380., 405./
160      REAL PVteta(ngridmx,ntetaSTD)
161c
162      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
163      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
164c
165     
166      REAL SSUM
167
168      LOGICAL firstcal, debut
169      DATA firstcal/.true./
170      SAVE firstcal,debut
171!      REAL rdayvrai
172      REAL, intent(in):: jD_cur, jH_cur
173c
174c-----------------------------------------------------------------------
175c
176c    1. Initialisations :
177c    --------------------
178c
179c
180      IF ( firstcal )  THEN
181        debut = .TRUE.
182        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
183         PRINT*,'STOP dans calfis'
184         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
185         PRINT*,'  ngridmx  jjm   iim   '
186         PRINT*,ngridmx,jjm,iim
187         STOP
188        ENDIF
189      ELSE
190        debut = .FALSE.
191      ENDIF ! of IF (firstcal)
192
193c
194c
195c-----------------------------------------------------------------------
196c   40. transformation des variables dynamiques en variables physiques:
197c   ---------------------------------------------------------------
198
199c   41. pressions au sol (en Pascals)
200c   ----------------------------------
201
202       
203      zpsrf(1) = pps(1,1)
204
205      ig0  = 2
206      DO j = 2,jjm
207         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
208         ig0 = ig0+iim
209      ENDDO
210
211      zpsrf(ngridmx) = pps(1,jjp1)
212
213
214c   42. pression intercouches :
215c
216c   -----------------------------------------------------------------
217c     .... zplev  definis aux (llm +1) interfaces des couches  ....
218c     .... zplay  definis aux (  llm )    milieux des couches  ....
219c   -----------------------------------------------------------------
220
221c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
222c
223       unskap   = 1./ kappa
224c
225      DO l = 1, llmp1
226        zplev( 1,l ) = pp(1,1,l)
227        ig0 = 2
228          DO j = 2, jjm
229             DO i =1, iim
230              zplev( ig0,l ) = pp(i,j,l)
231              ig0 = ig0 +1
232             ENDDO
233          ENDDO
234        zplev( ngridmx,l ) = pp(1,jjp1,l)
235      ENDDO
236c
237c
238
239c   43. temperature naturelle (en K) et pressions milieux couches .
240c   ---------------------------------------------------------------
241
242      DO l=1,llm
243
244         pksurcp     =  ppk(1,1,l) / cpp
245         zplay(1,l)  =  preff * pksurcp ** unskap
246         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
247         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
248         ig0         = 2
249
250         DO j = 2, jjm
251            DO i = 1, iim
252              pksurcp        = ppk(i,j,l) / cpp
253              zplay(ig0,l)   = preff * pksurcp ** unskap
254              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
255              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
256              ig0            = ig0 + 1
257            ENDDO
258         ENDDO
259
260         pksurcp       = ppk(1,jjp1,l) / cpp
261         zplay(ig0,l)  = preff * pksurcp ** unskap
262         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
263         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
264
265      ENDDO
266
267c   43.bis traceurs
268c   ---------------
269c
270      DO iq=1,nqtot
271          iiq=niadv(iq)
272         DO l=1,llm
273            zqfi(1,l,iq) = pq(1,1,l,iiq)
274            ig0          = 2
275            DO j=2,jjm
276               DO i = 1, iim
277                  zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
278                  ig0             = ig0 + 1
279               ENDDO
280            ENDDO
281            zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
282         ENDDO
283      ENDDO
284
285c   convergence dynamique pour les traceurs "EAU"
286! Earth-specific treatment of first 2 tracers (water)
287       if (planet_type=="earth") then
288        DO iq=1,2
289         DO l=1,llm
290            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
291            ig0          = 2
292            DO j=2,jjm
293               DO i = 1, iim
294                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
295                  ig0             = ig0 + 1
296               ENDDO
297            ENDDO
298            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
299         ENDDO
300        ENDDO
301       endif ! of if (planet_type=="earth")
302
303
304c   Geopotentiel calcule par rapport a la surface locale:
305c   -----------------------------------------------------
306
307      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
308      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
309      DO l=1,llm
310         DO ig=1,ngridmx
311           zphi(ig,l)=zphi(ig,l)-zphis(ig)
312         ENDDO
313      ENDDO
314
315c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
316c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
317c    de masse est calclue dans advtrac.F 
318c      DO l=1,llm
319c        pvervel(1,l)=pw(1,1,l) * g /apoln
320c        ig0=2
321c       DO j=2,jjm
322c           DO i = 1, iim
323c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
324c              ig0 = ig0 + 1
325c           ENDDO
326c       ENDDO
327c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
328c      ENDDO
329
330c
331c   45. champ u:
332c   ------------
333
334      DO 50 l=1,llm
335
336         DO 25 j=2,jjm
337            ig0 = 1+(j-2)*iim
338            zufi(ig0+1,l)= 0.5 *
339     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
340            pcvgu(ig0+1,l)= 0.5 *
341     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
342            DO 10 i=2,iim
343               zufi(ig0+i,l)= 0.5 *
344     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
345               pcvgu(ig0+i,l)= 0.5 *
346     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
34710         CONTINUE
34825      CONTINUE
349
35050    CONTINUE
351
352
353c   46.champ v:
354c   -----------
355
356      DO l=1,llm
357         DO j=2,jjm
358            ig0=1+(j-2)*iim
359            DO i=1,iim
360               zvfi(ig0+i,l)= 0.5 *
361     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
362               pcvgv(ig0+i,l)= 0.5 *
363     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
364            ENDDO
365         ENDDO
366      ENDDO
367
368
369c   47. champs de vents aux pole nord   
370c   ------------------------------
371c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
372c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
373
374      DO l=1,llm
375
376         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
377         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
378         DO i=2,iim
379            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
380            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
381         ENDDO
382
383         DO i=1,iim
384            zcos(i)   = COS(rlonv(i))*z1(i)
385            zcosbis(i)= COS(rlonv(i))*z1bis(i)
386            zsin(i)   = SIN(rlonv(i))*z1(i)
387            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
388         ENDDO
389
390         zufi(1,l)  = SSUM(iim,zcos,1)/pi
391         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
392         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
393         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
394
395      ENDDO
396
397
398c   48. champs de vents aux pole sud:
399c   ---------------------------------
400c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
401c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
402
403      DO l=1,llm
404
405         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
406         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
407         DO i=2,iim
408            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
409            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
410         ENDDO
411
412         DO i=1,iim
413            zcos(i)    = COS(rlonv(i))*z1(i)
414            zcosbis(i) = COS(rlonv(i))*z1bis(i)
415            zsin(i)    = SIN(rlonv(i))*z1(i)
416            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
417         ENDDO
418
419         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
420         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
421         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
422         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
423
424      ENDDO
425c
426      if (planet_type=="earth") then
427#ifdef CPP_EARTH
428cIM calcul PV a teta=350, 380, 405K
429      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
430     $           ztfi,zplay,zplev,
431     $           ntetaSTD,rtetaSTD,PVteta)
432#endif
433      endif
434c
435c On change de grille, dynamique vers physiq, pour le flux de masse verticale
436      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
437
438c-----------------------------------------------------------------------
439c   Appel de la physique:
440c   ---------------------
441
442
443      if (planet_type=="earth") then
444#ifdef CPP_EARTH
445      CALL physiq (ngridmx,
446     .             llm,
447     .             debut,
448     .             lafin,
449     .             jD_cur,
450     .             jH_cur,
451     .             dtphys,
452     .             zplev,
453     .             zplay,
454     .             zphi,
455     .             zphis,
456     .             presnivs,
457     .             clesphy0,
458     .             zufi,
459     .             zvfi,
460     .             ztfi,
461     .             zqfi,
462     .             flxwfi,
463     .             zdufi,
464     .             zdvfi,
465     .             zdtfi,
466     .             zdqfi,
467     .             zdpsrf,
468cIM diagnostique PVteta, Amip2         
469     .             pducov,
470     .             PVteta)
471#endif
472      endif !of if (planet_type=="earth")
473
474500   CONTINUE
475
476c-----------------------------------------------------------------------
477c   transformation des tendances physiques en tendances dynamiques:
478c   ---------------------------------------------------------------
479
480c  tendance sur la pression :
481c  -----------------------------------
482
483      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
484c
485c   62. enthalpie potentielle
486c   ---------------------
487
488      DO l=1,llm
489
490         DO i=1,iip1
491          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
492          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
493         ENDDO
494
495         DO j=2,jjm
496            ig0=1+(j-2)*iim
497            DO i=1,iim
498               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
499            ENDDO
500               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
501         ENDDO
502
503      ENDDO
504
505
506c   62. humidite specifique
507c   ---------------------
508! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
509!      DO iq=1,nqtot
510!         DO l=1,llm
511!            DO i=1,iip1
512!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
513!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
514!            ENDDO
515!            DO j=2,jjm
516!               ig0=1+(j-2)*iim
517!               DO i=1,iim
518!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
519!               ENDDO
520!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
521!            ENDDO
522!         ENDDO
523!      ENDDO
524
525c   63. traceurs
526c   ------------
527C     initialisation des tendances
528      pdqfi(:,:,:,:)=0.
529C
530      DO iq=1,nqtot
531         iiq=niadv(iq)
532         DO l=1,llm
533            DO i=1,iip1
534               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
535               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
536            ENDDO
537            DO j=2,jjm
538               ig0=1+(j-2)*iim
539               DO i=1,iim
540                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
541               ENDDO
542               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
543            ENDDO
544         ENDDO
545      ENDDO
546
547c   65. champ u:
548c   ------------
549
550      DO l=1,llm
551
552         DO i=1,iip1
553            pdufi(i,1,l)    = 0.
554            pdufi(i,jjp1,l) = 0.
555         ENDDO
556
557         DO j=2,jjm
558            ig0=1+(j-2)*iim
559            DO i=1,iim-1
560               pdufi(i,j,l)=
561     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
562            ENDDO
563            pdufi(iim,j,l)=
564     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
565            pdufi(iip1,j,l)=pdufi(1,j,l)
566         ENDDO
567
568      ENDDO
569
570
571c   67. champ v:
572c   ------------
573
574      DO l=1,llm
575
576         DO j=2,jjm-1
577            ig0=1+(j-2)*iim
578            DO i=1,iim
579               pdvfi(i,j,l)=
580     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
581            ENDDO
582            pdvfi(iip1,j,l) = pdvfi(1,j,l)
583         ENDDO
584      ENDDO
585
586
587c   68. champ v pres des poles:
588c   ---------------------------
589c      v = U * cos(long) + V * SIN(long)
590
591      DO l=1,llm
592
593         DO i=1,iim
594            pdvfi(i,1,l)=
595     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
596            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
597     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
598            pdvfi(i,1,l)=
599     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
600            pdvfi(i,jjm,l)=
601     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
602          ENDDO
603
604         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
605         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
606
607      ENDDO
608
609c-----------------------------------------------------------------------
610
611700   CONTINUE
612 
613      firstcal = .FALSE.
614
615      RETURN
616      END
Note: See TracBrowser for help on using the repository browser.