source: LMDZ4/branches/LMDZ4V5.0-LF/libf/dyn3d/calfis.F @ 5214

Last change on this file since 5214 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper


General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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