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

Last change on this file since 937 was 697, checked in by emillour, 13 years ago

Mars GCM:

  • Minor cleanup in surfini.F
  • Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)

EM

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