source: LMDZ5/branches/testing/libf/dynlonlat_phylonlat/calfis.F @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

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