source: trunk/LMDZ.MARS/libf/dynphy_lonlat/calfis.F @ 1547

Last change on this file since 1547 was 1547, checked in by emillour, 9 years ago

Mars GCM: in local dynamics package

  • bug fix in calfis: wrong array (pw) sent to physics: the transfered mass flux should be on the physics grid, not the dynamics grid. Moreover values at the poles needed to be correctly recomputed.

JL+EM

File size: 15.5 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      USE comvert_mod, ONLY: preff
10      USE comconst_mod, ONLY: dtphys,cpp,kappa,pi
11
12      IMPLICIT NONE
13c=======================================================================
14c
15c   1. rearrangement des tableaux et transformation
16c      variables dynamiques  >  variables physiques
17c   2. calcul des termes physiques
18c   3. retransformation des tendances physiques en tendances dynamiques
19c
20c   remarques:
21c   ----------
22c
23c    - les vents sont donnes dans la physique par leurs composantes
24c      naturelles.
25c    - la variable thermodynamique de la physique est une variable
26c      intensive :   T
27c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
28c    - les deux seules variables dependant de la geometrie necessaires
29c      pour la physique sont la latitude pour le rayonnement et
30c      l'aire de la maille quand on veut integrer une grandeur
31c      horizontalement.
32c    - les points de la physique sont les points scalaires de la
33c      la dynamique; numerotation:
34c          1 pour le pole nord
35c          (jjm-1)*iim pour l'interieur du domaine
36c          ngridmx pour le pole sud
37c      ---> ngridmx=2+(jjm-1)*iim
38c
39c     Input :
40c     -------
41c       ecritphy        frequence d'ecriture (en jours)de histphy
42c       pucov           covariant zonal velocity
43c       pvcov           covariant meridional velocity
44c       pteta           potential temperature
45c       pps             surface pressure
46c       pmasse          masse d'air dans chaque maille
47c       pts             surface temperature  (K)
48c       pw              flux vertical (kg/s)
49c
50c    Output :
51c    --------
52c        pdufi          tendency for the natural zonal velocity (ms-1)
53c        pdvfi          tendency for the natural meridional velocity
54c        pdhfi          tendency for the potential temperature
55c        pdtsfi         tendency for the surface temperature
56c
57c        tracer         Call tracer in  gcm.F ? (decided in callphys.def)
58c
59c=======================================================================
60c
61c-----------------------------------------------------------------------
62c
63c    0.  Declarations :
64c    ------------------
65
66#include "dimensions.h"
67#include "paramet.h"
68
69      INTEGER ngridmx,nq
70      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
71
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,nq)
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,nq)
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,nq)
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,nq)
116c
117!      REAL zvervel(ngridmx,llm)
118      REAL flxwfi(ngridmx,llm) ! vertical mass flux (kg/s) on physics grid
119c
120      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
121      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nq)
122      REAL zdpsrf(ngridmx)
123c
124      REAL zsin(iim),zcos(iim),z1(iim)
125      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
126      REAL unskap, pksurcp
127c
128     
129      EXTERNAL gr_dyn_fi,gr_fi_dyn
130      EXTERNAL physiq,multipl
131      REAL SSUM
132      EXTERNAL SSUM
133
134      REAL latfi(ngridmx),lonfi(ngridmx)
135      REAL airefi(ngridmx)
136      SAVE latfi, lonfi, airefi
137
138      LOGICAL firstcal, debut
139      DATA firstcal/.true./
140      SAVE firstcal,debut
141      REAL rdayvrai,rday_ecri
142c
143c-----------------------------------------------------------------------
144c
145c    1. Initialisations :
146c    --------------------
147c
148
149
150      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
151         PRINT*,'STOP dans calfis'
152         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
153         PRINT*,'  ngridmx  jjm   iim   '
154         PRINT*,ngridmx,jjm,iim
155         STOP
156      ENDIF
157
158c-----------------------------------------------------------------------
159c   latitude, longitude et aires des mailles pour la physique:
160c   ----------------------------------------------------------
161
162c
163      IF ( firstcal )  THEN
164          debut = .TRUE.
165      ELSE
166          debut = .FALSE.
167      ENDIF
168
169c
170!      IF (firstcal) THEN
171!         latfi(1)=rlatu(1)
172!         lonfi(1)=0.
173!         DO j=2,jjm
174!            DO i=1,iim
175!               latfi((j-2)*iim+1+i)= rlatu(j)
176!               lonfi((j-2)*iim+1+i)= rlonv(i)
177!            ENDDO
178!         ENDDO
179!         latfi(ngridmx)= rlatu(jjp1)
180!         lonfi(ngridmx)= 0.
181!         
182!         ! build airefi(), mesh area on physics grid
183!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
184!         ! Poles are single points on physics grid
185!         airefi(1)=airefi(1)*iim
186!         airefi(ngridmx)=airefi(ngridmx)*iim
187!
188!         CALL inifis(ngridmx,llm,nq,day_ini,daysec,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         ig0         =  2
246
247         DO j = 2, jjm
248            DO i = 1, iim
249              pksurcp        = ppk(i,j,l) / cpp
250              zplay(ig0,l)   = preff * pksurcp ** unskap
251              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
252              ig0            = ig0 + 1
253            ENDDO
254         ENDDO
255
256         pksurcp       = ppk(1,jjp1,l) / cpp
257         zplay(ig0,l)  = preff * pksurcp ** unskap
258         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
259
260      ENDDO
261
262      DO l=1, llm
263        DO ig=1,ngridmx
264             if (ztfi(ig,l).lt.15) then
265                  write(*,*) 'New Temperature below 15 K !!! '
266                      write(*,*) 'Stop in calfis.F '
267                      write(*,*) 'ig=', ig, ' l=', l
268                      write(*,*) 'ztfi(ig,l)=',ztfi(ig,l)
269                  stop
270             end if
271        ENDDO
272      ENDDO
273
274
275
276c   43.bis Taceurs (en kg/kg)
277c   --------------------------
278      DO iq=1,nq
279         DO l=1,llm
280            zqfi(1,l,iq) = pq(1,1,l,iq)
281            ig0          = 2
282            DO j=2,jjm
283               DO i = 1, iim
284                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
285                  ig0             = ig0 + 1
286               ENDDO
287            ENDDO
288            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
289         ENDDO
290      ENDDO
291
292c   Geopotentiel calcule par rapport a la surface locale:
293c   -----------------------------------------------------
294
295      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
296      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
297      DO l=1,llm
298         DO ig=1,ngridmx
299            zphi(ig,l)=zphi(ig,l)-zphis(ig)
300         ENDDO
301      ENDDO
302
303c   Calcul de la vitesse  verticale (m/s) pour diagnostique
304c   -------------------------------
305c     pw est en kg/s
306c On interpole "lineairement" la temperature entre les couches(FF,10/95)
307
308!      DO ig=1,ngridmx
309!         zvervel(ig,1)=0.
310!      END DO
311!      DO l=2,llm
312!        zvervel(1,l)=(pw(1,1,l)/apoln)
313!     &  * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l)             
314!        ig0=2
315!       DO j=2,jjm
316!           DO i = 1, iim
317!              zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
318!     &        * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
319!              ig0 = ig0 + 1
320!           ENDDO
321!       ENDDO
322!        zvervel(ig0,l)=(pw(1,jjp1,l)/apols)
323!     &  * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
324!      ENDDO
325
326c    .........  Reindexation : calcul de zvervel au MILIEU des couches
327!       DO l=1,llm-1
328!             DO ig=1,ngridmx
329!                    zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1))
330!          END DO
331!       END DO
332c      (dans la couche llm, on garde la valeur à la limite inférieure llm)
333
334! vertical mass flux
335      ! tranfer values from dynamics grid to physics grid:
336      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pw,flxwfi)
337      ! but mass flux is an extensive variable, so take the sum at the poles
338      DO l=1,llm
339        flxwfi(1,l)=sum(pw(1:iim,1,l))
340        flxwfi(ngridmx,l)=sum(pw(1:iim,jjp1,l))
341      ENDDO
342
343c   45. champ u:
344c   ------------
345
346      DO l=1,llm
347         DO j=2,jjm
348            ig0 = 1+(j-2)*iim
349            zufi(ig0+1,l)= 0.5 *
350     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
351            DO i=2,iim
352               zufi(ig0+i,l)= 0.5 *
353     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
354            ENDDO
355        ENDDO
356      ENDDO
357
358
359c   46.champ v:
360c   -----------
361
362      DO l=1,llm
363         DO j=2,jjm
364            ig0=1+(j-2)*iim
365            DO i=1,iim
366               zvfi(ig0+i,l)= 0.5 *
367     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
368            ENDDO
369         ENDDO
370      ENDDO
371
372
373c   47. champs de vents aux pole nord   
374c   ------------------------------
375c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
376c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
377
378      DO l=1,llm
379
380         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
381         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
382         DO i=2,iim
383            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
384            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
385         ENDDO
386
387         DO i=1,iim
388            zcos(i)   = COS(rlonv(i))*z1(i)
389            zcosbis(i)= COS(rlonv(i))*z1bis(i)
390            zsin(i)   = SIN(rlonv(i))*z1(i)
391            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
392         ENDDO
393
394         zufi(1,l)  = SSUM(iim,zcos,1)/pi
395         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
396
397      ENDDO
398
399
400c   48. champs de vents aux pole sud:
401c   ---------------------------------
402c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
403c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
404
405      DO l=1,llm
406
407         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
408         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
409         DO i=2,iim
410            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
411            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
412      ENDDO
413
414         DO i=1,iim
415            zcos(i)    = COS(rlonv(i))*z1(i)
416            zcosbis(i) = COS(rlonv(i))*z1bis(i)
417            zsin(i)    = SIN(rlonv(i))*z1(i)
418            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
419      ENDDO
420
421         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
422         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
423
424      ENDDO
425
426c-----------------------------------------------------------------------
427c   Appel de la physique:
428c   ---------------------
429
430
431      CALL physiq (ngridmx,llm,nq,
432     ,     debut,lafin,
433     ,     rday_ecri,heure,dtphys,
434     ,     zplev,zplay,zphi,
435     ,     zufi, zvfi,ztfi, zqfi, 
436!     ,     zvervel,
437     ,     flxwfi,
438C - sorties
439     s     zdufi, zdvfi, zdtfi, zdqfi,zdpsrf,tracer)
440
441
442c-----------------------------------------------------------------------
443c   transformation des tendances physiques en tendances dynamiques:
444c   ---------------------------------------------------------------
445
446c  tendance sur la pression :
447c  -----------------------------------
448
449      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
450c
451ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
452
453c   62. enthalpie potentielle
454c   ---------------------
455
456      DO l=1,llm
457
458         DO i=1,iip1
459          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
460          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
461         ENDDO
462
463         DO j=2,jjm
464            ig0=1+(j-2)*iim
465            DO i=1,iim
466               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
467            ENDDO
468               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
469         ENDDO
470
471      ENDDO
472
473
474c   62. humidite specifique
475c   ---------------------
476
477      DO iq=1,nq
478         DO l=1,llm
479            DO i=1,iip1
480               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
481               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
482            ENDDO
483            DO j=2,jjm
484               ig0=1+(j-2)*iim
485               DO i=1,iim
486                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
487               ENDDO
488               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
489            ENDDO
490         ENDDO
491      ENDDO
492
493c   65. champ u:
494c   ------------
495
496      DO l=1,llm
497
498         DO i=1,iip1
499            pdufi(i,1,l)    = 0.
500            pdufi(i,jjp1,l) = 0.
501         ENDDO
502
503         DO j=2,jjm
504            ig0=1+(j-2)*iim
505            DO i=1,iim-1
506               pdufi(i,j,l)=
507     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
508            ENDDO
509            pdufi(iim,j,l)=
510     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
511            pdufi(iip1,j,l)=pdufi(1,j,l)
512         ENDDO
513
514      ENDDO
515
516
517c   67. champ v:
518c   ------------
519
520      DO l=1,llm
521
522         DO j=2,jjm-1
523            ig0=1+(j-2)*iim
524            DO i=1,iim
525               pdvfi(i,j,l)=
526     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
527            ENDDO
528            pdvfi(iip1,j,l) = pdvfi(1,j,l)
529         ENDDO
530      ENDDO
531
532
533c   68. champ v pres des poles:
534c   ---------------------------
535c      v = U * cos(long) + V * SIN(long)
536
537      DO l=1,llm
538
539         DO i=1,iim
540            pdvfi(i,1,l)=
541     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
542            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
543     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
544            pdvfi(i,1,l)=
545     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
546            pdvfi(i,jjm,l)=
547     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
548          ENDDO
549
550         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
551         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
552
553      ENDDO
554
555c-----------------------------------------------------------------------
556
557700   CONTINUE
558
559      firstcal = .FALSE.
560
561      RETURN
562      END
Note: See TracBrowser for help on using the repository browser.