source: trunk/LMDZ.MARS/libf/dyn3d/calfis.F @ 620

Last change on this file since 620 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
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
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                      write(*,*) 'ztfi(ig,l)=',ztfi(ig,l)
263                  stop
264             end if
265        ENDDO
266      ENDDO
267
268
269
270c   43.bis Taceurs (en kg/kg)
271c   --------------------------
272      DO iq=1,nqmx
273         DO l=1,llm
274            zqfi(1,l,iq) = pq(1,1,l,iq)
275            ig0          = 2
276            DO j=2,jjm
277               DO i = 1, iim
278                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
279                  ig0             = ig0 + 1
280               ENDDO
281            ENDDO
282            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
283         ENDDO
284      ENDDO
285
286c   Geopotentiel calcule par rapport a la surface locale:
287c   -----------------------------------------------------
288
289      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
290      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
291      DO l=1,llm
292         DO ig=1,ngridmx
293            zphi(ig,l)=zphi(ig,l)-zphis(ig)
294         ENDDO
295      ENDDO
296
297c   Calcul de la vitesse  verticale (m/s) pour diagnostique
298c   -------------------------------
299c     pw est en kg/s
300c On interpole "lineairement" la temperature entre les couches(FF,10/95)
301
302      DO ig=1,ngridmx
303         zvervel(ig,1)=0.
304      END DO
305      DO l=2,llm
306        zvervel(1,l)=(pw(1,1,l)/apoln)
307     &  * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l)             
308        ig0=2
309       DO j=2,jjm
310           DO i = 1, iim
311              zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
312     &        * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
313              ig0 = ig0 + 1
314           ENDDO
315       ENDDO
316        zvervel(ig0,l)=(pw(1,jjp1,l)/apols)
317     &  * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
318      ENDDO
319
320c    .........  Reindexation : calcul de zvervel au MILIEU des couches
321       DO l=1,llm-1
322              DO ig=1,ngridmx
323                     zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1))
324          END DO
325       END DO
326c      (dans la couche llm, on garde la valeur à la limite inférieure llm)
327
328c   45. champ u:
329c   ------------
330
331      DO 50 l=1,llm
332
333         DO 25 j=2,jjm
334            ig0 = 1+(j-2)*iim
335            zufi(ig0+1,l)= 0.5 *
336     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
337            DO 10 i=2,iim
338               zufi(ig0+i,l)= 0.5 *
339     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
34010         CONTINUE
34125      CONTINUE
342
34350    CONTINUE
344
345
346c   46.champ v:
347c   -----------
348
349      DO l=1,llm
350         DO j=2,jjm
351            ig0=1+(j-2)*iim
352            DO i=1,iim
353               zvfi(ig0+i,l)= 0.5 *
354     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
355            ENDDO
356         ENDDO
357      ENDDO
358
359
360c   47. champs de vents aux pole nord   
361c   ------------------------------
362c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
363c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
364
365      DO l=1,llm
366
367         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
368         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
369         DO i=2,iim
370            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
371            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
372         ENDDO
373
374         DO i=1,iim
375            zcos(i)   = COS(rlonv(i))*z1(i)
376            zcosbis(i)= COS(rlonv(i))*z1bis(i)
377            zsin(i)   = SIN(rlonv(i))*z1(i)
378            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
379         ENDDO
380
381         zufi(1,l)  = SSUM(iim,zcos,1)/pi
382         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
383
384      ENDDO
385
386
387c   48. champs de vents aux pole sud:
388c   ---------------------------------
389c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
390c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
391
392      DO l=1,llm
393
394         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
395         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
396         DO i=2,iim
397            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
398            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
399      ENDDO
400
401         DO i=1,iim
402            zcos(i)    = COS(rlonv(i))*z1(i)
403            zcosbis(i) = COS(rlonv(i))*z1bis(i)
404            zsin(i)    = SIN(rlonv(i))*z1(i)
405            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
406      ENDDO
407
408         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
409         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
410
411      ENDDO
412
413c-----------------------------------------------------------------------
414c   Appel de la physique:
415c   ---------------------
416
417
418      CALL physiq (ngridmx,llm,nq,
419     ,     debut,lafin,
420     ,     rday_ecri,heure,dtphys,
421     ,     zplev,zplay,zphi,
422     ,     zufi, zvfi,ztfi, zqfi, 
423     ,     zvervel,
424C - sorties
425     s     zdufi, zdvfi, zdtfi, zdqfi,zdpsrf,tracer)
426
427
428c-----------------------------------------------------------------------
429c   transformation des tendances physiques en tendances dynamiques:
430c   ---------------------------------------------------------------
431
432c  tendance sur la pression :
433c  -----------------------------------
434
435      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
436c
437ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
438
439c   62. enthalpie potentielle
440c   ---------------------
441
442      DO l=1,llm
443
444         DO i=1,iip1
445          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
446          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
447         ENDDO
448
449         DO j=2,jjm
450            ig0=1+(j-2)*iim
451            DO i=1,iim
452               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
453            ENDDO
454               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
455         ENDDO
456
457      ENDDO
458
459
460c   62. humidite specifique
461c   ---------------------
462
463      DO iq=1,nqmx
464         DO l=1,llm
465            DO i=1,iip1
466               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
467               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
468            ENDDO
469            DO j=2,jjm
470               ig0=1+(j-2)*iim
471               DO i=1,iim
472                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
473               ENDDO
474               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
475            ENDDO
476         ENDDO
477      ENDDO
478
479c   65. champ u:
480c   ------------
481
482      DO l=1,llm
483
484         DO i=1,iip1
485            pdufi(i,1,l)    = 0.
486            pdufi(i,jjp1,l) = 0.
487         ENDDO
488
489         DO j=2,jjm
490            ig0=1+(j-2)*iim
491            DO i=1,iim-1
492               pdufi(i,j,l)=
493     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
494            ENDDO
495            pdufi(iim,j,l)=
496     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
497            pdufi(iip1,j,l)=pdufi(1,j,l)
498         ENDDO
499
500      ENDDO
501
502
503c   67. champ v:
504c   ------------
505
506      DO l=1,llm
507
508         DO j=2,jjm-1
509            ig0=1+(j-2)*iim
510            DO i=1,iim
511               pdvfi(i,j,l)=
512     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
513            ENDDO
514            pdvfi(iip1,j,l) = pdvfi(1,j,l)
515         ENDDO
516      ENDDO
517
518
519c   68. champ v pres des poles:
520c   ---------------------------
521c      v = U * cos(long) + V * SIN(long)
522
523      DO l=1,llm
524
525         DO i=1,iim
526            pdvfi(i,1,l)=
527     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
528            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
529     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
530            pdvfi(i,1,l)=
531     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
532            pdvfi(i,jjm,l)=
533     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
534          ENDDO
535
536         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
537         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
538
539      ENDDO
540
541c-----------------------------------------------------------------------
542
543700   CONTINUE
544
545      firstcal = .FALSE.
546
547      RETURN
548      END
Note: See TracBrowser for help on using the repository browser.