source: trunk/MESOSCALE/LMDZ.MARS/libf_gcm/dyn3d/calfis.F @ 815

Last change on this file since 815 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

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
185
186c
187c-----------------------------------------------------------------------
188c   40. transformation des variables dynamiques en variables physiques:
189c   ---------------------------------------------------------------
190
191c   41. pressions au sol (en Pascals)
192c   ----------------------------------
193
194       
195      zpsrf(1) = pps(1,1)
196
197      ig0  = 2
198      DO j = 2,jjm
199         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
200         ig0 = ig0+iim
201      ENDDO
202
203      zpsrf(ngridmx) = pps(1,jjp1)
204
205
206c   42. pression intercouches :
207c
208c   -----------------------------------------------------------------
209c     .... zplev  definis aux (llm +1) interfaces des couches  ....
210c     .... zplay  definis aux (  llm )    milieux des couches  ....
211c   -----------------------------------------------------------------
212
213c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
214c
215       unskap   = 1./ kappa
216c
217      DO l = 1, llmp1
218        zplev( 1,l ) = pp(1,1,l)
219        ig0 = 2
220          DO j = 2, jjm
221             DO i =1, iim
222              zplev( ig0,l ) = pp(i,j,l)
223              ig0 = ig0 +1
224             ENDDO
225          ENDDO
226        zplev( ngridmx,l ) = pp(1,jjp1,l)
227      ENDDO
228c
229c
230
231c   43. temperature naturelle (en K) et pressions milieux couches .
232c   ---------------------------------------------------------------
233
234      DO l=1,llm
235
236         pksurcp     =  ppk(1,1,l) / cpp
237         zplay(1,l)  =  preff * pksurcp ** unskap
238         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
239         ig0         =  2
240
241         DO j = 2, jjm
242            DO i = 1, iim
243              pksurcp        = ppk(i,j,l) / cpp
244              zplay(ig0,l)   = preff * pksurcp ** unskap
245              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
246              ig0            = ig0 + 1
247            ENDDO
248         ENDDO
249
250         pksurcp       = ppk(1,jjp1,l) / cpp
251         zplay(ig0,l)  = preff * pksurcp ** unskap
252         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
253
254      ENDDO
255
256      DO l=1, llm
257        DO ig=1,ngridmx
258             if (ztfi(ig,l).lt.15) then
259                  write(*,*) 'New Temperature below 15 K !!! '
260                      write(*,*) 'Stop in calfis.F '
261                      write(*,*) 'ig=', ig, 'l= ', l
262                  stop
263             end if
264        ENDDO
265      ENDDO
266
267
268
269c   43.bis Taceurs (en kg/kg)
270c   --------------------------
271      DO iq=1,nqmx
272         DO l=1,llm
273            zqfi(1,l,iq) = pq(1,1,l,iq)
274            ig0          = 2
275            DO j=2,jjm
276               DO i = 1, iim
277                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
278                  ig0             = ig0 + 1
279               ENDDO
280            ENDDO
281            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
282         ENDDO
283      ENDDO
284
285c   Geopotentiel calcule par rapport a la surface locale:
286c   -----------------------------------------------------
287
288      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
289      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
290      DO l=1,llm
291         DO ig=1,ngridmx
292            zphi(ig,l)=zphi(ig,l)-zphis(ig)
293         ENDDO
294      ENDDO
295
296c   Calcul de la vitesse  verticale (m/s) pour diagnostique
297c   -------------------------------
298c     pw est en kg/s
299c On interpole "lineairement" la temperature entre les couches(FF,10/95)
300
301      DO ig=1,ngridmx
302         zvervel(ig,1)=0.
303      END DO
304      DO l=2,llm
305        zvervel(1,l)=(pw(1,1,l)/apoln)
306     &  * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l)             
307        ig0=2
308       DO j=2,jjm
309           DO i = 1, iim
310              zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
311     &        * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
312              ig0 = ig0 + 1
313           ENDDO
314       ENDDO
315        zvervel(ig0,l)=(pw(1,jjp1,l)/apols)
316     &  * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
317      ENDDO
318
319c    .........  Reindexation : calcul de zvervel au MILIEU des couches
320       DO l=1,llm-1
321              DO ig=1,ngridmx
322                     zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1))
323          END DO
324       END DO
325c      (dans la couche llm, on garde la valeur à la limite inférieure llm)
326
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            DO 10 i=2,iim
337               zufi(ig0+i,l)= 0.5 *
338     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
33910         CONTINUE
34025      CONTINUE
341
34250    CONTINUE
343
344
345c   46.champ v:
346c   -----------
347
348      DO l=1,llm
349         DO j=2,jjm
350            ig0=1+(j-2)*iim
351            DO i=1,iim
352               zvfi(ig0+i,l)= 0.5 *
353     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
354            ENDDO
355         ENDDO
356      ENDDO
357
358
359c   47. champs de vents aux pole nord   
360c   ------------------------------
361c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
362c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
363
364      DO l=1,llm
365
366         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
367         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
368         DO i=2,iim
369            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
370            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
371         ENDDO
372
373         DO i=1,iim
374            zcos(i)   = COS(rlonv(i))*z1(i)
375            zcosbis(i)= COS(rlonv(i))*z1bis(i)
376            zsin(i)   = SIN(rlonv(i))*z1(i)
377            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
378         ENDDO
379
380         zufi(1,l)  = SSUM(iim,zcos,1)/pi
381         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
382
383      ENDDO
384
385
386c   48. champs de vents aux pole sud:
387c   ---------------------------------
388c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
389c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
390
391      DO l=1,llm
392
393         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
394         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
395         DO i=2,iim
396            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
397            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
398      ENDDO
399
400         DO i=1,iim
401            zcos(i)    = COS(rlonv(i))*z1(i)
402            zcosbis(i) = COS(rlonv(i))*z1bis(i)
403            zsin(i)    = SIN(rlonv(i))*z1(i)
404            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
405      ENDDO
406
407         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
408         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
409
410      ENDDO
411
412c-----------------------------------------------------------------------
413c   Appel de la physique:
414c   ---------------------
415
416
417      CALL physiq (ngridmx,llm,nq,
418     ,     debut,lafin,
419     ,     rday_ecri,heure,dtphys,
420     ,     zplev,zplay,zphi,
421     ,     zufi, zvfi,ztfi, zqfi, 
422     ,     zvervel,
423C - sorties
424     s     zdufi, zdvfi, zdtfi, zdqfi,zdpsrf,tracer)
425
426
427c-----------------------------------------------------------------------
428c   transformation des tendances physiques en tendances dynamiques:
429c   ---------------------------------------------------------------
430
431c  tendance sur la pression :
432c  -----------------------------------
433
434      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
435c
436ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
437
438c   62. enthalpie potentielle
439c   ---------------------
440
441      DO l=1,llm
442
443         DO i=1,iip1
444          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
445          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
446         ENDDO
447
448         DO j=2,jjm
449            ig0=1+(j-2)*iim
450            DO i=1,iim
451               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
452            ENDDO
453               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
454         ENDDO
455
456      ENDDO
457
458
459c   62. humidite specifique
460c   ---------------------
461
462      DO iq=1,nqmx
463         DO l=1,llm
464            DO i=1,iip1
465               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
466               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
467            ENDDO
468            DO j=2,jjm
469               ig0=1+(j-2)*iim
470               DO i=1,iim
471                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
472               ENDDO
473               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
474            ENDDO
475         ENDDO
476      ENDDO
477
478c   65. champ u:
479c   ------------
480
481      DO l=1,llm
482
483         DO i=1,iip1
484            pdufi(i,1,l)    = 0.
485            pdufi(i,jjp1,l) = 0.
486         ENDDO
487
488         DO j=2,jjm
489            ig0=1+(j-2)*iim
490            DO i=1,iim-1
491               pdufi(i,j,l)=
492     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
493            ENDDO
494            pdufi(iim,j,l)=
495     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
496            pdufi(iip1,j,l)=pdufi(1,j,l)
497         ENDDO
498
499      ENDDO
500
501
502c   67. champ v:
503c   ------------
504
505      DO l=1,llm
506
507         DO j=2,jjm-1
508            ig0=1+(j-2)*iim
509            DO i=1,iim
510               pdvfi(i,j,l)=
511     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
512            ENDDO
513            pdvfi(iip1,j,l) = pdvfi(1,j,l)
514         ENDDO
515      ENDDO
516
517
518c   68. champ v pres des poles:
519c   ---------------------------
520c      v = U * cos(long) + V * SIN(long)
521
522      DO l=1,llm
523
524         DO i=1,iim
525            pdvfi(i,1,l)=
526     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
527            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
528     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
529            pdvfi(i,1,l)=
530     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
531            pdvfi(i,jjm,l)=
532     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
533          ENDDO
534
535         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
536         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
537
538      ENDDO
539
540c-----------------------------------------------------------------------
541
542700   CONTINUE
543
544      firstcal = .FALSE.
545
546      RETURN
547      END
Note: See TracBrowser for help on using the repository browser.