source: LMDZ.3.3/branches/LF/libf/dyn3d/calfis.F @ 5007

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

Initial revision

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