source: trunk/LMDZ.PLUTO.old/libf/dyn3d/calfis.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 14.8 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,
4     $            pdufi,pdvfi,pdhfi,pdqfi,pdpsfi,tracer)
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       pw              flux vertical (kg m-2)
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        tracer         Call tracer in  gcm.F ? (decided in callphys.def)
55c
56c=======================================================================
57c
58c-----------------------------------------------------------------------
59c
60c    0.  Declarations :
61c    ------------------
62
63#include "dimensions.h"
64#include "paramet.h"
65#include "temps.h"
66
67      INTEGER ngridmx,nq
68      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
69
70#include "comconst.h"
71#include "comvert.h"
72#include "comgeom2.h"
73#include "control.h"
74
75c    Arguments :
76c    -----------
77      LOGICAL  lafin
78      REAL heure
79
80      REAL pvcov(iip1,jjm,llm)
81      REAL pucov(iip1,jjp1,llm)
82      REAL pteta(iip1,jjp1,llm)
83      REAL pmasse(iip1,jjp1,llm)
84      REAL pq(iip1,jjp1,llm,nqmx)
85      REAL pphis(iip1,jjp1)
86      REAL pphi(iip1,jjp1,llm)
87c
88      REAL pdvcov(iip1,jjm,llm)
89      REAL pducov(iip1,jjp1,llm)
90      REAL pdteta(iip1,jjp1,llm)
91      REAL pdq(iip1,jjp1,llm,nqmx)
92c
93      REAL pw(iip1,jjp1,llm)
94c
95      REAL pps(iip1,jjp1)
96      REAL pp(iip1,jjp1,llmp1)
97      REAL ppk(iip1,jjp1,llm)
98c
99      REAL pdvfi(iip1,jjm,llm)
100      REAL pdufi(iip1,jjp1,llm)
101      REAL pdhfi(iip1,jjp1,llm)
102      REAL pdqfi(iip1,jjp1,llm,nqmx)
103      REAL pdpsfi(iip1,jjp1)
104      logical tracer
105
106c    Local variables :
107c    -----------------
108
109      INTEGER i,j,l,ig0,ig,iq
110      REAL zpsrf(ngridmx)
111      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
112      REAL zphi(ngridmx,llm),zphis(ngridmx)
113c
114      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
115      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
116c
117      REAL zvervel(ngridmx,llm)
118c
119      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
120      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
121      REAL zdpsrf(ngridmx)
122c
123      REAL zsin(iim),zcos(iim),z1(iim)
124      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
125      REAL unskap, pksurcp
126c
127     
128      EXTERNAL gr_dyn_fi,gr_fi_dyn
129      EXTERNAL physiq,multipl
130      REAL SSUM
131      EXTERNAL SSUM
132
133      REAL latfi(ngridmx),lonfi(ngridmx)
134      REAL airefi(ngridmx)
135      SAVE latfi, lonfi, airefi
136
137      LOGICAL firstcal, debut
138      DATA firstcal/.true./
139      SAVE firstcal,debut
140      REAL rdayvrai,rday_ecri
141c
142c-----------------------------------------------------------------------
143c
144c    1. Initialisations :
145c    --------------------
146c
147
148
149      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
150         PRINT*,'STOP dans calfis'
151         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
152         PRINT*,'  ngridmx  jjm   iim   '
153         PRINT*,ngridmx,jjm,iim
154         STOP
155      ENDIF
156
157c-----------------------------------------------------------------------
158c   latitude, longitude et aires des mailles pour la physique:
159c   ----------------------------------------------------------
160
161c
162      IF ( firstcal )  THEN
163          debut = .TRUE.
164      ELSE
165          debut = .FALSE.
166      ENDIF
167
168c
169      IF (firstcal) THEN
170         latfi(1)=rlatu(1)
171         lonfi(1)=0.
172         DO j=2,jjm
173            DO i=1,iim
174               latfi((j-2)*iim+1+i)= rlatu(j)
175               lonfi((j-2)*iim+1+i)= rlonv(i)
176            ENDDO
177         ENDDO
178         latfi(ngridmx)= rlatu(jjp1)
179         lonfi(ngridmx)= 0.
180         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
181
182         CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
183     .                latfi,lonfi,airefi,rad,g,r,cpp)
184      ENDIF
185c
186c-----------------------------------------------------------------------
187c   40. transformation des variables dynamiques en variables physiques:
188c   ---------------------------------------------------------------
189
190c   41. pressions au sol (en Pascals)
191c   ----------------------------------
192
193      zpsrf(1) = pps(1,1)
194
195      ig0  = 2
196      DO j = 2,jjm
197         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
198         ig0 = ig0+iim
199      ENDDO
200
201      zpsrf(ngridmx) = pps(1,jjp1)
202
203
204c   42. pression intercouches :
205c
206c   -----------------------------------------------------------------
207c     .... zplev  definis aux (llm +1) interfaces des couches  ....
208c     .... zplay  definis aux (  llm )    milieux des couches  ....
209c   -----------------------------------------------------------------
210
211c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
212c
213       unskap   = 1./ kappa
214c
215      DO l = 1, llmp1
216        zplev( 1,l ) = pp(1,1,l)
217        ig0 = 2
218          DO j = 2, jjm
219             DO i =1, iim
220              zplev( ig0,l ) = pp(i,j,l)
221              ig0 = ig0 +1
222             ENDDO
223          ENDDO
224        zplev( ngridmx,l ) = pp(1,jjp1,l)
225      ENDDO
226
227c   43. temperature naturelle (en K) et pressions milieux couches .
228c   ---------------------------------------------------------------
229
230      DO l=1,llm
231
232         pksurcp     =  ppk(1,1,l) / cpp
233         zplay(1,l)  =  preff * pksurcp ** unskap
234         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
235         ig0         =  2
236
237         DO j = 2, jjm
238            DO i = 1, iim
239              pksurcp        = ppk(i,j,l) / cpp
240              zplay(ig0,l)   = preff * pksurcp ** unskap
241              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
242              ig0            = ig0 + 1
243            ENDDO
244         ENDDO
245
246         pksurcp       = ppk(1,jjp1,l) / cpp
247         zplay(ig0,l)  = preff * pksurcp ** unskap
248         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
249
250      ENDDO
251
252!      DO l=1, llm
253!        DO ig=1,ngridmx
254!             if (ztfi(ig,l).lt.10) then
255!                  write(*,*) 'New Temperature below 10 K !!! '
256!                 write(*,*) 'Stop in calfis.F '
257!                  write(*,*) 'ig=', ig, ' l=', l
258!                  write(*,*) 'ztfi(ig,l)=',ztfi(ig,l)
259!                  stop
260!             end if
261!        ENDDO
262!      ENDDO
263
264
265
266c   43.bis Taceurs (en kg/kg)
267c   --------------------------
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   Geopotentiel calcule par rapport a la surface locale:
283c   -----------------------------------------------------
284
285      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
286      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
287      DO l=1,llm
288         DO ig=1,ngridmx
289            zphi(ig,l)=zphi(ig,l)-zphis(ig)
290         ENDDO
291      ENDDO
292
293c   Calcul de la vitesse  verticale (m/s) pour diagnostique
294c   -------------------------------
295c     pw est en kg/s
296c On interpole "lineairement" la temperature entre les couches(FF,10/95)
297
298      DO ig=1,ngridmx
299         zvervel(ig,1)=0.
300      END DO
301      DO l=2,llm
302        zvervel(1,l)=(pw(1,1,l)/apoln)
303     &  * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l)             
304        ig0=2
305       DO j=2,jjm
306           DO i = 1, iim
307              zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
308     &        * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
309              ig0 = ig0 + 1
310           ENDDO
311       ENDDO
312        zvervel(ig0,l)=(pw(1,jjp1,l)/apols)
313     &  * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
314      ENDDO
315
316c    .........  Reindexation : calcul de zvervel au MILIEU des couches
317       DO l=1,llm-1
318              DO ig=1,ngridmx
319                     zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1))
320          END DO
321       END DO
322c      (dans la couche llm, on garde la valeur à la limite inférieure llm)
323
324c   45. champ u:
325c   ------------
326
327      DO 50 l=1,llm
328
329         DO 25 j=2,jjm
330            ig0 = 1+(j-2)*iim
331            zufi(ig0+1,l)= 0.5 *
332     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
333            DO 10 i=2,iim
334               zufi(ig0+i,l)= 0.5 *
335     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
33610         CONTINUE
33725      CONTINUE
338
33950    CONTINUE
340
341
342c   46.champ v:
343c   -----------
344
345      DO l=1,llm
346         DO j=2,jjm
347            ig0=1+(j-2)*iim
348            DO i=1,iim
349               zvfi(ig0+i,l)= 0.5 *
350     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
351            ENDDO
352         ENDDO
353      ENDDO
354
355
356c   47. champs de vents aux pole nord   
357c   ------------------------------
358c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
359c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
360
361      DO l=1,llm
362
363         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
364         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
365         DO i=2,iim
366            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
367            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
368         ENDDO
369
370         DO i=1,iim
371            zcos(i)   = COS(rlonv(i))*z1(i)
372            zcosbis(i)= COS(rlonv(i))*z1bis(i)
373            zsin(i)   = SIN(rlonv(i))*z1(i)
374            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
375         ENDDO
376
377         zufi(1,l)  = SSUM(iim,zcos,1)/pi
378         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
379
380      ENDDO
381
382
383c   48. champs de vents aux pole sud:
384c   ---------------------------------
385c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
386c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
387
388      DO l=1,llm
389
390         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
391         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
392         DO i=2,iim
393            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
394            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
395      ENDDO
396
397         DO i=1,iim
398            zcos(i)    = COS(rlonv(i))*z1(i)
399            zcosbis(i) = COS(rlonv(i))*z1bis(i)
400            zsin(i)    = SIN(rlonv(i))*z1(i)
401            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
402      ENDDO
403
404         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
405         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
406
407      ENDDO
408
409c-----------------------------------------------------------------------
410c   Appel de la physique:
411c   ---------------------
412
413      !print*, 'TB18 ztfi=',ztfi
414      CALL physiq (ngridmx,llm,nq,
415     ,     debut,lafin,
416     ,     rday_ecri,heure,dtphys,
417     ,     zplev,zplay,zphi,
418     ,     zufi, zvfi,ztfi, zqfi, 
419     ,     zvervel,
420C - sorties
421     s     zdufi, zdvfi, zdtfi, zdqfi,zdpsrf,tracer)
422
423
424c-----------------------------------------------------------------------
425c   transformation des tendances physiques en tendances dynamiques:
426c   ---------------------------------------------------------------
427
428c  tendance sur la pression :
429c  -----------------------------------
430
431      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
432c
433ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
434
435c   62. enthalpie potentielle
436c   ---------------------
437
438      DO l=1,llm
439
440         DO i=1,iip1
441          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
442          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
443         ENDDO
444
445         DO j=2,jjm
446            ig0=1+(j-2)*iim
447            DO i=1,iim
448               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
449            ENDDO
450               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
451         ENDDO
452
453      ENDDO
454
455
456c   62. humidite specifique
457c   ---------------------
458
459      DO iq=1,nqmx
460         DO l=1,llm
461            DO i=1,iip1
462               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
463               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
464            ENDDO
465            DO j=2,jjm
466               ig0=1+(j-2)*iim
467               DO i=1,iim
468                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
469               ENDDO
470               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
471            ENDDO
472         ENDDO
473      ENDDO
474
475c   65. champ u:
476c   ------------
477
478      DO l=1,llm
479
480         DO i=1,iip1
481            pdufi(i,1,l)    = 0.
482            pdufi(i,jjp1,l) = 0.
483         ENDDO
484
485         DO j=2,jjm
486            ig0=1+(j-2)*iim
487            DO i=1,iim-1
488               pdufi(i,j,l)=
489     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
490            ENDDO
491            pdufi(iim,j,l)=
492     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
493            pdufi(iip1,j,l)=pdufi(1,j,l)
494         ENDDO
495
496      ENDDO
497
498
499c   67. champ v:
500c   ------------
501
502      DO l=1,llm
503
504         DO j=2,jjm-1
505            ig0=1+(j-2)*iim
506            DO i=1,iim
507               pdvfi(i,j,l)=
508     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
509            ENDDO
510            pdvfi(iip1,j,l) = pdvfi(1,j,l)
511         ENDDO
512      ENDDO
513
514
515c   68. champ v pres des poles:
516c   ---------------------------
517c      v = U * cos(long) + V * SIN(long)
518
519      DO l=1,llm
520
521         DO i=1,iim
522            pdvfi(i,1,l)=
523     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
524            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
525     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
526            pdvfi(i,1,l)=
527     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
528            pdvfi(i,jjm,l)=
529     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
530          ENDDO
531
532         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
533         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
534
535      ENDDO
536
537c-----------------------------------------------------------------------
538
539700   CONTINUE
540
541      firstcal = .FALSE.
542
543      RETURN
544      END
Note: See TracBrowser for help on using the repository browser.