source: trunk/libf/dyn3d/calfis.F @ 1

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

Import initial LMDZ5

File size: 18.1 KB
RevLine 
[1]1!
2! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
3!
4C
5C
6      SUBROUTINE calfis(lafin,
7     $                  jD_cur, jH_cur,
8     $                  pucov,
9     $                  pvcov,
10     $                  pteta,
11     $                  pq,
12     $                  pmasse,
13     $                  pps,
14     $                  pp,
15     $                  ppk,
16     $                  pphis,
17     $                  pphi,
18     $                  pducov,
19     $                  pdvcov,
20     $                  pdteta,
21     $                  pdq,
22     $                  flxw,
23     $                  clesphy0,
24     $                  pdufi,
25     $                  pdvfi,
26     $                  pdhfi,
27     $                  pdqfi,
28     $                  pdpsfi)
29c
30c    Auteur :  P. Le Van, F. Hourdin
31c   .........
32      USE infotrac
33      USE control_mod
34 
35
36      IMPLICIT NONE
37c=======================================================================
38c
39c   1. rearrangement des tableaux et transformation
40c      variables dynamiques  >  variables physiques
41c   2. calcul des termes physiques
42c   3. retransformation des tendances physiques en tendances dynamiques
43c
44c   remarques:
45c   ----------
46c
47c    - les vents sont donnes dans la physique par leurs composantes
48c      naturelles.
49c    - la variable thermodynamique de la physique est une variable
50c      intensive :   T
51c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
52c    - les deux seules variables dependant de la geometrie necessaires
53c      pour la physique sont la latitude pour le rayonnement et
54c      l'aire de la maille quand on veut integrer une grandeur
55c      horizontalement.
56c    - les points de la physique sont les points scalaires de la
57c      la dynamique; numerotation:
58c          1 pour le pole nord
59c          (jjm-1)*iim pour l'interieur du domaine
60c          ngridmx pour le pole sud
61c      ---> ngridmx=2+(jjm-1)*iim
62c
63c     Input :
64c     -------
65c       pucov           covariant zonal velocity
66c       pvcov           covariant meridional velocity
67c       pteta           potential temperature
68c       pps             surface pressure
69c       pmasse          masse d'air dans chaque maille
70c       pts             surface temperature  (K)
71c       callrad         clef d'appel au rayonnement
72c
73c    Output :
74c    --------
75c        pdufi          tendency for the natural zonal velocity (ms-1)
76c        pdvfi          tendency for the natural meridional velocity
77c        pdhfi          tendency for the potential temperature
78c        pdtsfi         tendency for the surface temperature
79c
80c        pdtrad         radiative tendencies  \  both input
81c        pfluxrad       radiative fluxes      /  and output
82c
83c=======================================================================
84c
85c-----------------------------------------------------------------------
86c
87c    0.  Declarations :
88c    ------------------
89
90#include "dimensions.h"
91#include "paramet.h"
92#include "temps.h"
93
94      INTEGER ngridmx
95      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
96
97#include "comconst.h"
98#include "comvert.h"
99#include "comgeom2.h"
100#include "iniprint.h"
101
102c    Arguments :
103c    -----------
104      LOGICAL  lafin
105
106
107      REAL pvcov(iip1,jjm,llm)
108      REAL pucov(iip1,jjp1,llm)
109      REAL pteta(iip1,jjp1,llm)
110      REAL pmasse(iip1,jjp1,llm)
111      REAL pq(iip1,jjp1,llm,nqtot)
112      REAL pphis(iip1,jjp1)
113      REAL pphi(iip1,jjp1,llm)
114c
115      REAL pdvcov(iip1,jjm,llm)
116      REAL pducov(iip1,jjp1,llm)
117      REAL pdteta(iip1,jjp1,llm)
118      REAL pdq(iip1,jjp1,llm,nqtot)
119c
120      REAL pps(iip1,jjp1)
121      REAL pp(iip1,jjp1,llmp1)
122      REAL ppk(iip1,jjp1,llm)
123c
124      REAL pdvfi(iip1,jjm,llm)
125      REAL pdufi(iip1,jjp1,llm)
126      REAL pdhfi(iip1,jjp1,llm)
127      REAL pdqfi(iip1,jjp1,llm,nqtot)
128      REAL pdpsfi(iip1,jjp1)
129
130      INTEGER        longcles
131      PARAMETER    ( longcles = 20 )
132      REAL clesphy0( longcles )
133
134
135c    Local variables :
136c    -----------------
137
138      INTEGER i,j,l,ig0,ig,iq,iiq
139      REAL zpsrf(ngridmx)
140      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
141      REAL zphi(ngridmx,llm),zphis(ngridmx)
142c
143      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
144      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
145c
146      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
147      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
148c
149      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
150      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
151      REAL zdpsrf(ngridmx)
152c
153      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
154      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
155      REAL jH_cur_split,zdt_split
156      LOGICAL debut_split,lafin_split
157      INTEGER isplit
158
159      REAL zsin(iim),zcos(iim),z1(iim)
160      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
161      REAL unskap, pksurcp
162c
163cIM diagnostique PVteta, Amip2
164      INTEGER ntetaSTD
165      PARAMETER(ntetaSTD=3)
166      REAL rtetaSTD(ntetaSTD)
167      DATA rtetaSTD/350., 380., 405./
168      REAL PVteta(ngridmx,ntetaSTD)
169c
170      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
171      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
172c
173     
174      REAL SSUM
175
176      LOGICAL firstcal, debut
177      DATA firstcal/.true./
178      SAVE firstcal,debut
179!      REAL rdayvrai
180      REAL, intent(in):: jD_cur, jH_cur
181c
182c-----------------------------------------------------------------------
183c
184c    1. Initialisations :
185c    --------------------
186c
187c
188      IF ( firstcal )  THEN
189        debut = .TRUE.
190        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
191         write(lunout,*) 'STOP dans calfis'
192         write(lunout,*)
193     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
194         write(lunout,*) '  ngridmx  jjm   iim   '
195         write(lunout,*) ngridmx,jjm,iim
196         STOP
197        ENDIF
198      ELSE
199        debut = .FALSE.
200      ENDIF ! of IF (firstcal)
201
202c
203c
204c-----------------------------------------------------------------------
205c   40. transformation des variables dynamiques en variables physiques:
206c   ---------------------------------------------------------------
207
208c   41. pressions au sol (en Pascals)
209c   ----------------------------------
210
211       
212      zpsrf(1) = pps(1,1)
213
214      ig0  = 2
215      DO j = 2,jjm
216         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
217         ig0 = ig0+iim
218      ENDDO
219
220      zpsrf(ngridmx) = pps(1,jjp1)
221
222
223c   42. pression intercouches :
224c
225c   -----------------------------------------------------------------
226c     .... zplev  definis aux (llm +1) interfaces des couches  ....
227c     .... zplay  definis aux (  llm )    milieux des couches  ....
228c   -----------------------------------------------------------------
229
230c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
231c
232       unskap   = 1./ kappa
233c
234      DO l = 1, llmp1
235        zplev( 1,l ) = pp(1,1,l)
236        ig0 = 2
237          DO j = 2, jjm
238             DO i =1, iim
239              zplev( ig0,l ) = pp(i,j,l)
240              ig0 = ig0 +1
241             ENDDO
242          ENDDO
243        zplev( ngridmx,l ) = pp(1,jjp1,l)
244      ENDDO
245c
246c
247
248c   43. temperature naturelle (en K) et pressions milieux couches .
249c   ---------------------------------------------------------------
250
251      DO l=1,llm
252
253         pksurcp     =  ppk(1,1,l) / cpp
254         zplay(1,l)  =  preff * pksurcp ** unskap
255         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
256         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
257         ig0         = 2
258
259         DO j = 2, jjm
260            DO i = 1, iim
261              pksurcp        = ppk(i,j,l) / cpp
262              zplay(ig0,l)   = preff * pksurcp ** unskap
263              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
264              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
265              ig0            = ig0 + 1
266            ENDDO
267         ENDDO
268
269         pksurcp       = ppk(1,jjp1,l) / cpp
270         zplay(ig0,l)  = preff * pksurcp ** unskap
271         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
272         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
273
274      ENDDO
275
276c   43.bis traceurs
277c   ---------------
278c
279      DO iq=1,nqtot
280          iiq=niadv(iq)
281         DO l=1,llm
282            zqfi(1,l,iq) = pq(1,1,l,iiq)
283            ig0          = 2
284            DO j=2,jjm
285               DO i = 1, iim
286                  zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
287                  ig0             = ig0 + 1
288               ENDDO
289            ENDDO
290            zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
291         ENDDO
292      ENDDO
293
294c   convergence dynamique pour les traceurs "EAU"
295! Earth-specific treatment of first 2 tracers (water)
296       if (planet_type=="earth") then
297        DO iq=1,2
298         DO l=1,llm
299            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
300            ig0          = 2
301            DO j=2,jjm
302               DO i = 1, iim
303                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
304                  ig0             = ig0 + 1
305               ENDDO
306            ENDDO
307            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
308         ENDDO
309        ENDDO
310       endif ! of if (planet_type=="earth")
311
312
313c   Geopotentiel calcule par rapport a la surface locale:
314c   -----------------------------------------------------
315
316      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
317      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
318      DO l=1,llm
319         DO ig=1,ngridmx
320           zphi(ig,l)=zphi(ig,l)-zphis(ig)
321         ENDDO
322      ENDDO
323
324c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
325c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
326c    de masse est calclue dans advtrac.F 
327c      DO l=1,llm
328c        pvervel(1,l)=pw(1,1,l) * g /apoln
329c        ig0=2
330c       DO j=2,jjm
331c           DO i = 1, iim
332c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
333c              ig0 = ig0 + 1
334c           ENDDO
335c       ENDDO
336c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
337c      ENDDO
338
339c
340c   45. champ u:
341c   ------------
342
343      DO 50 l=1,llm
344
345         DO 25 j=2,jjm
346            ig0 = 1+(j-2)*iim
347            zufi(ig0+1,l)= 0.5 *
348     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
349            pcvgu(ig0+1,l)= 0.5 *
350     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
351            DO 10 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               pcvgu(ig0+i,l)= 0.5 *
355     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
35610         CONTINUE
35725      CONTINUE
358
35950    CONTINUE
360
361
362c   46.champ v:
363c   -----------
364
365      DO l=1,llm
366         DO j=2,jjm
367            ig0=1+(j-2)*iim
368            DO i=1,iim
369               zvfi(ig0+i,l)= 0.5 *
370     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
371               pcvgv(ig0+i,l)= 0.5 *
372     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
373            ENDDO
374         ENDDO
375      ENDDO
376
377
378c   47. champs de vents aux pole nord   
379c   ------------------------------
380c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
381c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
382
383      DO l=1,llm
384
385         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
386         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
387         DO i=2,iim
388            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
389            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
390         ENDDO
391
392         DO i=1,iim
393            zcos(i)   = COS(rlonv(i))*z1(i)
394            zcosbis(i)= COS(rlonv(i))*z1bis(i)
395            zsin(i)   = SIN(rlonv(i))*z1(i)
396            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
397         ENDDO
398
399         zufi(1,l)  = SSUM(iim,zcos,1)/pi
400         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
401         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
402         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
403
404      ENDDO
405
406
407c   48. champs de vents aux pole sud:
408c   ---------------------------------
409c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
410c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
411
412      DO l=1,llm
413
414         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
415         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
416         DO i=2,iim
417            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
418            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
419         ENDDO
420
421         DO i=1,iim
422            zcos(i)    = COS(rlonv(i))*z1(i)
423            zcosbis(i) = COS(rlonv(i))*z1bis(i)
424            zsin(i)    = SIN(rlonv(i))*z1(i)
425            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
426         ENDDO
427
428         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
429         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
430         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
431         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
432
433      ENDDO
434c
435      if (planet_type=="earth") then
436#ifdef CPP_EARTH
437cIM calcul PV a teta=350, 380, 405K
438      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
439     $           ztfi,zplay,zplev,
440     $           ntetaSTD,rtetaSTD,PVteta)
441#endif
442      endif
443c
444c On change de grille, dynamique vers physiq, pour le flux de masse verticale
445      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
446
447c-----------------------------------------------------------------------
448c   Appel de la physique:
449c   ---------------------
450
451
452      if (planet_type=="earth") then
453#ifdef CPP_EARTH
454
455!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
456      zdt_split=dtphys/nsplit_phys
457      zdufic(:,:)=0.
458      zdvfic(:,:)=0.
459      zdtfic(:,:)=0.
460      zdqfic(:,:,:)=0.
461
462      do isplit=1,nsplit_phys
463
464         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
465         debut_split=debut.and.isplit==1
466         lafin_split=lafin.and.isplit==nsplit_phys
467
468         CALL physiq (ngridmx,
469     .             llm,
470     .             debut_split,
471     .             lafin_split,
472     .             jD_cur,
473     .             jH_cur_split,
474     .             zdt_split,
475     .             zplev,
476     .             zplay,
477     .             zphi,
478     .             zphis,
479     .             presnivs,
480     .             clesphy0,
481     .             zufi,
482     .             zvfi,
483     .             ztfi,
484     .             zqfi,
485     .             flxwfi,
486     .             zdufi,
487     .             zdvfi,
488     .             zdtfi,
489     .             zdqfi,
490     .             zdpsrf,
491cIM diagnostique PVteta, Amip2         
492     .             pducov,
493     .             PVteta)
494
495         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
496         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
497         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
498         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
499
500         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
501         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
502         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
503         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
504
505      enddo
506      zdufi(:,:)=zdufic(:,:)/nsplit_phys
507      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
508      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
509      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
510
511#endif
512      endif !of if (planet_type=="earth")
513
514500   CONTINUE
515
516c-----------------------------------------------------------------------
517c   transformation des tendances physiques en tendances dynamiques:
518c   ---------------------------------------------------------------
519
520c  tendance sur la pression :
521c  -----------------------------------
522
523      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
524c
525c   62. enthalpie potentielle
526c   ---------------------
527
528      DO l=1,llm
529
530         DO i=1,iip1
531          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
532          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
533         ENDDO
534
535         DO j=2,jjm
536            ig0=1+(j-2)*iim
537            DO i=1,iim
538               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
539            ENDDO
540               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
541         ENDDO
542
543      ENDDO
544
545
546c   62. humidite specifique
547c   ---------------------
548! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
549!      DO iq=1,nqtot
550!         DO l=1,llm
551!            DO i=1,iip1
552!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
553!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
554!            ENDDO
555!            DO j=2,jjm
556!               ig0=1+(j-2)*iim
557!               DO i=1,iim
558!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
559!               ENDDO
560!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
561!            ENDDO
562!         ENDDO
563!      ENDDO
564
565c   63. traceurs
566c   ------------
567C     initialisation des tendances
568      pdqfi(:,:,:,:)=0.
569C
570      DO iq=1,nqtot
571         iiq=niadv(iq)
572         DO l=1,llm
573            DO i=1,iip1
574               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
575               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
576            ENDDO
577            DO j=2,jjm
578               ig0=1+(j-2)*iim
579               DO i=1,iim
580                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
581               ENDDO
582               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
583            ENDDO
584         ENDDO
585      ENDDO
586
587c   65. champ u:
588c   ------------
589
590      DO l=1,llm
591
592         DO i=1,iip1
593            pdufi(i,1,l)    = 0.
594            pdufi(i,jjp1,l) = 0.
595         ENDDO
596
597         DO j=2,jjm
598            ig0=1+(j-2)*iim
599            DO i=1,iim-1
600               pdufi(i,j,l)=
601     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
602            ENDDO
603            pdufi(iim,j,l)=
604     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
605            pdufi(iip1,j,l)=pdufi(1,j,l)
606         ENDDO
607
608      ENDDO
609
610
611c   67. champ v:
612c   ------------
613
614      DO l=1,llm
615
616         DO j=2,jjm-1
617            ig0=1+(j-2)*iim
618            DO i=1,iim
619               pdvfi(i,j,l)=
620     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
621            ENDDO
622            pdvfi(iip1,j,l) = pdvfi(1,j,l)
623         ENDDO
624      ENDDO
625
626
627c   68. champ v pres des poles:
628c   ---------------------------
629c      v = U * cos(long) + V * SIN(long)
630
631      DO l=1,llm
632
633         DO i=1,iim
634            pdvfi(i,1,l)=
635     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
636            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
637     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
638            pdvfi(i,1,l)=
639     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
640            pdvfi(i,jjm,l)=
641     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
642          ENDDO
643
644         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
645         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
646
647      ENDDO
648
649c-----------------------------------------------------------------------
650
651700   CONTINUE
652 
653      firstcal = .FALSE.
654
655      RETURN
656      END
Note: See TracBrowser for help on using the repository browser.