source: LMDZ.3.3/trunk/libf/dyn3d/inigeom.F @ 2

Last change on this file since 2 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 21.8 KB
Line 
1      SUBROUTINE inigeom
2c
3c     Auteur :  P. Le Van
4c    .....................
5c
6c   ............      Version  du 20/12/98     ........................
7c
8c  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-
9c     endroits que les aires aireij1,..aireij4 .
10c  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
11c
12c
13      IMPLICIT NONE
14c
15      REAL        deltay, tauy, taux
16      PARAMETER ( deltay =  0. , tauy = 1. , taux = 1. )
17c
18c     deltay  est ( en degres ) le deplacement eventuel en Y du zoom
19cc    taux  et  tauy  sont  les raideurs   du  zoom
20c
21c
22#include "dimensions.h"
23#include "paramet.h"
24#include "comconst.h"
25#include "comgeom2.h"
26#include "serre.h"
27#include "logic.h"
28#include "comdissnew.h"
29
30c-----------------------------------------------------------------------
31c   ....  Variables  locales   ....
32c
33      INTEGER  i,j,itmax,itmay,iter
34      REAL cvu(iip1,jjp1),cuv(iip1,jjm)
35      REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
36      REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
37      REAL coslatm,coslatp,radclatm,radclatp
38      REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
39     *     cuij4(iip1,jjp1)
40      REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
41     *     cvij4(iip1,jjp1)
42      REAL rlonvv(iip1),rlatuu(jjp1)
43      REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
44     *     yprimv(jjm),yprimu(jjp1)
45      REAL gamdi_gdiv, gamdi_grot, gamdi_h
46 
47      REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
48     ,  xprimp025(iip1)
49      SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
50      SAVE rlonm025,xprimm025,rlonp025,xprimp025
51
52
53
54c----------------------------------------------------------------------
55      REAL      SSUM
56      EXTERNAL  SSUM
57c
58c
59c
60c   ------------------------------------------------------------------
61c   -                                                                -
62c   -    calcul des coeff. ( cu, cv , 1./cu**2,  1./cv**2  )         -
63c   -                                                                -
64c   ------------------------------------------------------------------
65c   ------------------------------------------------------------------
66c
67c      les coef. ( cu, cv ) permettent de passer des vitesses naturelles
68c      aux vitesses covariantes et contravariantes , ou vice-versa ...
69c
70c
71c     on a :  u (covariant) = cu * u (naturel)   , u(contrav)= u(nat)/cu
72c             v (covariant) = cv * v (naturel)   , v(contrav)= v(nat)/cv
73c
74c       on en tire :  u(covariant) = cu * cu * u(contravariant)
75c                     v(covariant) = cv * cv * v(contravariant)
76c
77c
78c     on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2
79c                                                          =     =
80c                                           et   - jm/2    <  Y  < jm/2
81c                                                          =     =
82c
83c      ...................................................
84c      ...................................................
85c      .  x  est la longitude du point  en radians       .
86c      .  y  est la  latitude du point  en radians       .
87c      .                                                 .
88c      .  on a :  cu(i,j) = rad * COS(y) * dx/dX         .
89c      .          cv( j ) = rad          * dy/dY         .
90c      .        aire(i,j) =  cu(i,j) * cv(j)             .
91c      .                                                 .
92c      . y, dx/dX, dy/dY calcules aux points concernes   .
93c      .                                                 .
94c      ...................................................
95c      ...................................................
96c
97c
98c
99c                                                           ,
100c    cv , bien que dependant de j uniquement,sera ici indice aussi en i
101c          pour un adressage plus facile en  ij  .
102c
103c
104c
105c  **************  aux points  u  et  v ,           *****************
106c      xprimu et xprimv sont respectivement les valeurs de  dx/dX
107c      yprimu et yprimv    .  .  .  .  .  .  .  .  .  .  .  dy/dY
108c      rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude
109c      cvu    et   cv      .  .  .  .  .  .  .  .  .  .  .    cv
110c
111c  **************  aux points u, v, scalaires, et z  ****************
112c      cu, cuv, cuscal, cuz sont respectiv. les valeurs de    cu
113c
114c
115c
116c         Exemple de distribution de variables sur la grille dans le
117c             domaine de travail ( X,Y ) .
118c     ................................................................
119c                  DX=DY= 1
120c
121c   
122c        +     represente  un  point scalaire ( p.exp  la pression )
123c        >     represente  la composante zonale du  vent
124c        V     represente  la composante meridienne du vent
125c        o     represente  la  vorticite
126c
127c       ----  , car aux poles , les comp.zonales covariantes sont nulles
128c
129c
130c
131c         i ->
132c
133c         1      2      3      4      5      6      7      8
134c  j
135c  v  1   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
136c
137c         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
138c
139c     2   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
140c
141c         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
142c
143c     3   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
144c
145c         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
146c
147c     4   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
148c
149c         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
150c
151c     5   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
152c
153c
154c      Ci-dessus,  on voit que le nombre de pts.en longitude est egal
155c                 a   IM = 8
156c      De meme ,   le nombre d'intervalles entre les 2 poles est egal
157c                 a   JM = 4
158c
159c      Les points scalaires ( + ) correspondent donc a des valeurs
160c       entieres  de  i ( 1 a IM )   et  de  j ( 1 a  JM +1 )   .
161c
162c      Les vents    U       ( > ) correspondent a des valeurs  semi-
163c       entieres  de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
164c
165c      Les vents    V       ( V ) correspondent a des valeurs entieres
166c       de     i ( 1 a  IM ) et semi-entieres de  j ( 1 +0.5  a JM +0.5)
167c
168c
169c
170      PRINT 3
171 3    FORMAT( // 10x,' ....  INIGEOM  date du 01/06/98   ..... ',
172     * //5x,'   Calcul des elongations cu et cv  comme sommes des 4 ' /
173     *  5x,' elong. cuij1, .. 4  , cvij1,.. 4  qui les entourent , aux
174     * '/ 5x,' memes endroits que les aires aireij1,...j4   . ' / )
175c
176c
177      IF( nitergdiv.NE.2 ) THEN
178        gamdi_gdiv = coefdis/ ( float(nitergdiv) -2. )
179      ELSE
180        gamdi_gdiv = 0.
181      ENDIF
182      IF( nitergrot.NE.2 ) THEN
183        gamdi_grot = coefdis/ ( float(nitergrot) -2. )
184      ELSE
185        gamdi_grot = 0.
186      ENDIF
187      IF( niterh.NE.2 ) THEN
188        gamdi_h = coefdis/ ( float(niterh) -2. )
189      ELSE
190        gamdi_h = 0.
191      ENDIF
192
193      PRINT *,' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
194     *  nitergdiv,nitergrot,niterh
195c
196      pi    = 2.* ASIN(1.)
197c
198      PRINT 990
199
200c     ----------------------------------------------------------------
201c
202      IF( .NOT.fxyhypb )   THEN
203c
204c
205       IF( ysinus )  THEN
206c
207        PRINT *,' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
208c
209c   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
210
211        CALL  fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
212     ,                    rlatu2,yprimu2,
213     ,  rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
214
215       ELSE
216c
217        PRINT *,' *** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
218
219c  .... utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ...
220c
221 
222        pxo   = clon *pi /180.
223        pyo   = 2.* clat* pi /180.
224c
225c  ....  determination de  transx ( pour le zoom ) par Newton-Raphson ...
226c
227        itmax = 10
228        eps   = .1e-7
229c
230        xo1 = 0.
231        DO 10 iter = 1, itmax
232        x1  = xo1
233        f   = x1+ alphax *SIN(x1-pxo)
234        df  = 1.+ alphax *COS(x1-pxo)
235        x1  = x1 - f/df
236        xdm = ABS( x1- xo1 )
237        IF( xdm.LE.eps )GO TO 11
238        xo1 = x1
239 10     CONTINUE
240 11     CONTINUE
241c
242        transx = xo1
243
244        itmay = 10
245        eps   = .1e-7
246C
247        yo1  = 0.
248        DO 15 iter = 1,itmay
249        y1   = yo1
250        f    = y1 + alphay* SIN(y1-pyo)
251        df   = 1. + alphay* COS(y1-pyo)
252        y1   = y1 -f/df
253        ydm  = ABS(y1-yo1)
254        IF(ydm.LE.eps) GO TO 17
255        yo1  = y1
256 15     CONTINUE
257c
258 17     CONTINUE
259        transy = yo1
260
261        CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
262     ,              rlatu2,yprimu2,
263     ,  rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
264
265       ENDIF
266c
267      ELSE
268c
269c   ....  utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
270c   .....................................................................
271
272       PRINT *,'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
273
274       CALL fxyhyper( clat, grossismy, dzoomy, tauy, deltay  ,
275     ,               clon, grossismx, dzoomx, taux          ,
276     , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2  ,
277     , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
278
279 
280      ENDIF
281c
282c  -------------------------------------------------------------------
283
284c
285      rlatu(1)    =     ASIN(1.)
286      rlatu(jjp1) =  - rlatu(1)
287c
288c
289c   ....  calcul  aux  poles  ....
290c
291      yprimu(1)      = 0.
292      yprimu(jjp1)   = 0.
293c
294c
295      un4rad2 = 0.25 * rad * rad
296c
297c   --------------------------------------------------------------------
298c   --------------------------------------------------------------------
299c   -                                                                  -
300c   -  calcul  des aires ( aire,aireu,airev, 1./aire, 1./airez  )      -
301c   -      et de   fext ,  force de coriolis  extensive  .             -
302c   -                                                                  -
303c   --------------------------------------------------------------------
304c   --------------------------------------------------------------------
305c
306c
307c
308c   A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
309c   affectees 4 aires entourant P , calculees respectivement aux points
310c            ( i + 1/4, j - 1/4 )    :    aireij1 (i,j)
311c            ( i + 1/4, j + 1/4 )    :    aireij2 (i,j)
312c            ( i - 1/4, j + 1/4 )    :    aireij3 (i,j)
313c            ( i - 1/4, j - 1/4 )    :    aireij4 (i,j)
314c
315c           ,
316c   Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
317c   Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
318c   des 4 aires  aireij1,aireij2,aireij3,aireij4 qui sont affectees au
319c   point (i,j) .
320c   On definit en outre les coefficients  alpha comme etant egaux a
321c    (aireij / aire), c.a.d par exp.  alpha1(i,j)=aireij1(i,j)/aire(i,j)
322c
323c   De meme, toute aire centree en 1 point U est egale a la somme des
324c   4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
325c    Idem pour  airev, airez .
326c
327c       On a ,pour chaque maille :    dX = dY = 1
328c
329c
330c                             . V
331c
332c                 aireij4 .        . aireij1
333c
334c                   U .       . P      . U
335c
336c                 aireij3 .        . aireij2
337c
338c                             . V
339c
340c
341c
342c
343c
344c   ....................................................................
345c
346c    Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
347c   qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
348c   taires cuij et les 4 elongat. cvij qui sont calculees aux memes
349c     endroits  que les aireij   .   
350c
351c   ....................................................................
352c
353c     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....
354c
355c
356      DO 35 j = 1, jjp1
357c
358      IF ( j. eq. 1 )  THEN
359c
360      yprm           = yprimu1(j)
361      rlatm          = rlatu1(j)
362c
363      coslatm        = COS( rlatm )
364      radclatm       = 0.5* rad * coslatm
365c
366      DO 30 i = 1, iim
367      xprp           = xprimp025( i )
368      xprm           = xprimm025( i )
369      aireij2( i,1 ) = un4rad2 * coslatm  * xprp * yprm
370      aireij3( i,1 ) = un4rad2 * coslatm  * xprm * yprm
371      cuij2  ( i,1 ) = radclatm * xprp
372      cuij3  ( i,1 ) = radclatm * xprm
373      cvij2  ( i,1 ) = 0.5* rad * yprm
374      cvij3  ( i,1 ) = cvij2(i,1)
375  30  CONTINUE
376c
377      DO  i = 1, iim
378      aireij1( i,1 ) = 0.
379      aireij4( i,1 ) = 0.
380      cuij1  ( i,1 ) = 0.
381      cuij4  ( i,1 ) = 0.
382      cvij1  ( i,1 ) = 0.
383      cvij4  ( i,1 ) = 0.
384      ENDDO
385c
386      END IF
387c
388      IF ( j. eq. jjp1 )  THEN
389
390       yprp               = yprimu2(j-1)
391       rlatp              = rlatu2 (j-1)
392cc       yprp             = fyprim( FLOAT(j) - 0.25 )
393cc       rlatp            = fy    ( FLOAT(j) - 0.25 )
394c
395      coslatp             = COS( rlatp )
396      radclatp            = 0.5* rad * coslatp
397c
398      DO 31 i = 1,iim
399        xprp              = xprimp025( i )
400        xprm              = xprimm025( i )
401        aireij1( i,jjp1 ) = un4rad2 * coslatp  * xprp * yprp
402        aireij4( i,jjp1 ) = un4rad2 * coslatp  * xprm * yprp
403        cuij1(i,jjp1)     = radclatp * xprp
404        cuij4(i,jjp1)     = radclatp * xprm
405        cvij1(i,jjp1)     = 0.5 * rad* yprp
406        cvij4(i,jjp1)     = cvij1(i,jjp1)
407 31   CONTINUE
408c
409       DO   i    = 1, iim
410        aireij2( i,jjp1 ) = 0.
411        aireij3( i,jjp1 ) = 0.
412        cvij2  ( i,jjp1 ) = 0.
413        cvij3  ( i,jjp1 ) = 0.
414        cuij2  ( i,jjp1 ) = 0.
415        cuij3  ( i,jjp1 ) = 0.
416       ENDDO
417c
418      END IF
419c
420
421      IF ( j .gt. 1 .AND. j .lt. jjp1 )  THEN
422c
423        rlatp    = rlatu2 ( j-1 )
424        yprp     = yprimu2( j-1 )
425        rlatm    = rlatu1 (  j  )
426        yprm     = yprimu1(  j  )
427c         rlatp    = fy    ( FLOAT(j) - 0.25 )
428c         yprp     = fyprim( FLOAT(j) - 0.25 )
429c         rlatm    = fy    ( FLOAT(j) + 0.25 )
430c         yprm     = fyprim( FLOAT(j) + 0.25 )
431
432         coslatm  = COS( rlatm )
433         coslatp  = COS( rlatp )
434         radclatp = 0.5* rad * coslatp
435         radclatm = 0.5* rad * coslatm
436c
437         DO 32 i = 1,iim
438         xprp            = xprimp025( i )
439         xprm            = xprimm025( i )
440     
441         ai14            = un4rad2 * coslatp * yprp
442         ai23            = un4rad2 * coslatm * yprm
443         aireij1 ( i,j ) = ai14 * xprp
444         aireij2 ( i,j ) = ai23 * xprp
445         aireij3 ( i,j ) = ai23 * xprm
446         aireij4 ( i,j ) = ai14 * xprm
447         cuij1   ( i,j ) = radclatp * xprp
448         cuij2   ( i,j ) = radclatm * xprp
449         cuij3   ( i,j ) = radclatm * xprm
450         cuij4   ( i,j ) = radclatp * xprm
451         cvij1   ( i,j ) = 0.5* rad * yprp
452         cvij2   ( i,j ) = 0.5* rad * yprm
453         cvij3   ( i,j ) = cvij2(i,j)
454         cvij4   ( i,j ) = cvij1(i,j)
455  32     CONTINUE
456c
457      END IF
458c
459c    ........       periodicite   ............
460c
461         cvij1   (iip1,j) = cvij1   (1,j)
462         cvij2   (iip1,j) = cvij2   (1,j)
463         cvij3   (iip1,j) = cvij3   (1,j)
464         cvij4   (iip1,j) = cvij4   (1,j)
465         cuij1   (iip1,j) = cuij1   (1,j)
466         cuij2   (iip1,j) = cuij2   (1,j)
467         cuij3   (iip1,j) = cuij3   (1,j)
468         cuij4   (iip1,j) = cuij4   (1,j)
469         aireij1 (iip1,j) = aireij1 (1,j )
470         aireij2 (iip1,j) = aireij2 (1,j )
471         aireij3 (iip1,j) = aireij3 (1,j )
472         aireij4 (iip1,j) = aireij4 (1,j )
473       
474  35  CONTINUE
475c
476c    ..............................................................
477c
478      DO 37 j = 1, jjp1
479      DO 36 i = 1, iim
480      aire    ( i,j )  = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
481     *                          aireij4(i,j)
482      alpha1  ( i,j )  = aireij1(i,j) / aire(i,j)
483      alpha2  ( i,j )  = aireij2(i,j) / aire(i,j)
484      alpha3  ( i,j )  = aireij3(i,j) / aire(i,j)
485      alpha4  ( i,j )  = aireij4(i,j) / aire(i,j)
486      alpha1p2( i,j )  = alpha1 (i,j) + alpha2 (i,j)
487      alpha1p4( i,j )  = alpha1 (i,j) + alpha4 (i,j)
488      alpha2p3( i,j )  = alpha2 (i,j) + alpha3 (i,j)
489      alpha3p4( i,j )  = alpha3 (i,j) + alpha4 (i,j)
490  36  CONTINUE
491c
492c   ....  Modif  P. Le Van  ( 4/07/96 )  .....
493c
494      aire    (iip1,j) = aire    (1,j)
495      alpha1  (iip1,j) = alpha1  (1,j)
496      alpha2  (iip1,j) = alpha2  (1,j)
497      alpha3  (iip1,j) = alpha3  (1,j)
498      alpha4  (iip1,j) = alpha4  (1,j)
499      alpha1p2(iip1,j) = alpha1p2(1,j)
500      alpha1p4(iip1,j) = alpha1p4(1,j)
501      alpha2p3(iip1,j) = alpha2p3(1,j)
502      alpha3p4(iip1,j) = alpha3p4(1,j)
503  37  CONTINUE
504c
505
506      DO 42 j = 1,jjp1
507      DO 41 i = 1,iim
508      aireu       (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
509     *                                aireij3(i+1,j)
510      unsaire    ( i,j)= 1./ aire(i,j)
511      unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv )
512      unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h    )
513      airesurg   ( i,j)= aire(i,j)/ g
514  41  CONTINUE
515      aireu     (iip1,j)  = aireu  (1,j)
516      unsaire   (iip1,j)  = unsaire(1,j)
517      unsair_gam1(iip1,j) = unsair_gam1(1,j)
518      unsair_gam2(iip1,j) = unsair_gam2(1,j)
519      airesurg   (iip1,j) = airesurg(1,j)
520  42  CONTINUE
521c
522c
523      DO 48 j = 1,jjm
524c
525        DO i=1,iim
526         airev     (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
527     *                           aireij4(i,j+1)
528        ENDDO
529         DO i=1,iim
530          airez         = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) +
531     *                           aireij4(i+1,j+1)
532          unsairez(i,j) = 1./ airez
533          unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot )
534          fext    (i,j)   = airez * SIN(rlatv(j))* 2.* omeg
535         ENDDO
536        airev     (iip1,j)  = airev(1,j)
537        unsairez  (iip1,j)  = unsairez(1,j)
538        fext      (iip1,j)  = fext(1,j)
539        unsairz_gam(iip1,j) = unsairz_gam(1,j)
540c
541  48  CONTINUE
542c
543c
544c    .....      Calcul  des elongations cu,cv, cvu     .........
545c
546      DO    j   = 1, jjm
547       DO   i  = 1, iim
548       cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1))
549       cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j)  +cvij3(i,j) )
550       cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1))
551       unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) )
552       ENDDO
553       DO   i  = 1, iim
554       cuvsurcv (i,j)    = airev(i,j)  * unscv2(i,j)
555       cvsurcuv (i,j)    = 1./cuvsurcv(i,j)
556       cuvscvgam1(i,j)   = cuvsurcv (i,j) ** ( - gamdi_gdiv )
557       cuvscvgam2(i,j)   = cuvsurcv (i,j) ** ( - gamdi_h )
558       cvscuvgam(i,j)    = cvsurcuv (i,j) ** ( - gamdi_grot )
559       ENDDO
560       cv       (iip1,j)  = cv       (1,j)
561       cvu      (iip1,j)  = cvu      (1,j)
562       unscv2   (iip1,j)  = unscv2   (1,j)
563       cuv      (iip1,j)  = cuv      (1,j)
564       cuvsurcv (iip1,j)  = cuvsurcv (1,j)
565       cvsurcuv (iip1,j)  = cvsurcuv (1,j)
566       cuvscvgam1(iip1,j) = cuvscvgam1(1,j)
567       cuvscvgam2(iip1,j) = cuvscvgam2(1,j)
568       cvscuvgam(iip1,j)  = cvscuvgam(1,j)
569      ENDDO
570
571      DO  j     = 2, jjm
572        DO   i  = 1, iim
573        cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j))
574        unscu2    (i,j)  = 1./ ( cu(i,j) * cu(i,j) )
575        cvusurcu  (i,j)  =  aireu(i,j) * unscu2(i,j)
576        cusurcvu  (i,j)  = 1./ cvusurcu(i,j)
577        cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv )
578        cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h    )
579        cuscvugam (i,j)  = cusurcvu(i,j) ** ( - gamdi_grot )
580        ENDDO
581        cu       (iip1,j)  = cu(1,j)
582        unscu2   (iip1,j)  = unscu2(1,j)
583        cvusurcu (iip1,j)  = cvusurcu(1,j)
584        cusurcvu (iip1,j)  = cusurcvu(1,j)
585        cvuscugam1(iip1,j) = cvuscugam1(1,j)
586        cvuscugam2(iip1,j) = cvuscugam2(1,j)
587        cuscvugam (iip1,j) = cuscvugam(1,j)
588      ENDDO
589
590c
591c   ....  calcul aux  poles  ....
592c
593      DO    i      =  1, iip1
594        cu    ( i, 1 )  =   0.
595        unscu2( i, 1 )  =   0.
596        cvu   ( i, 1 )  =   0.
597c
598        cu    (i, jjp1) =   0.
599        unscu2(i, jjp1) =   0.
600        cvu   (i, jjp1) =   0.
601      ENDDO
602c
603c    ..............................................................
604c
605      DO j = 1, jjm
606        DO i= 1, iim
607         airvscu2  (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) )
608         aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot )
609        ENDDO
610         airvscu2  (iip1,j)  = airvscu2(1,j)
611         aivscu2gam(iip1,j)  = aivscu2gam(1,j)
612      ENDDO
613
614      DO j=2,jjm
615        DO i=1,iim
616         airuscv2   (i,j)    = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) )
617         aiuscv2gam (i,j)    = airuscv2(i,j)** ( - gamdi_grot )
618        ENDDO
619         airuscv2  (iip1,j)  = airuscv2  (1,j)
620         aiuscv2gam(iip1,j)  = aiuscv2gam(1,j)
621      ENDDO
622c
623c
624c   calcul des aires aux  poles :
625c   -----------------------------
626c
627      apoln       = SSUM(iim,aire(1,1),1)
628      apols       = SSUM(iim,aire(1,jjp1),1)
629      unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
630      unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
631      unsapolnga2 = 1./ ( apoln ** ( - gamdi_h    ) )
632      unsapolsga2 = 1./ ( apols ** ( - gamdi_h    ) )
633c
634c-----------------------------------------------------------------------
635c     gtitre='Coriolis version ancienne'
636c     gfichier='fext1'
637c     CALL writestd(fext,iip1*jjm)
638c
639c   changement F. Hourdin calcul conservatif pour fext
640c   constang contient le produit a * cos ( latitude ) * omega
641c
642      DO i=1,iim
643         constang(i,1) = 0.
644      ENDDO
645      DO j=1,jjm-1
646        DO i=1,iim
647         constang(i,j+1) = rad*omeg*cu(i,j+1)*COS(rlatu(j+1))
648        ENDDO
649      ENDDO
650      DO i=1,iim
651         constang(i,jjp1) = 0.
652      ENDDO
653c
654c   periodicite en longitude
655c
656      DO j=1,jjm
657        fext(iip1,j)     = fext(1,j)
658      ENDDO
659      DO j=1,jjp1
660        constang(iip1,j) = constang(1,j)
661      ENDDO
662
663c fin du changement
664
665c
666c-----------------------------------------------------------------------
667c
668       PRINT *,' INIGEOM  RLONV '
669        DO i=1,iip1
670         rlonvv(i) = rlonv(i)*180./pi
671        ENDDO
672       PRINT 400,rlonvv
673c
674       PRINT *,' RLATV '
675        DO i=1,jjm
676         rlatuu(i)=rlatv(i)*180./pi
677        ENDDO
678       PRINT 400,(rlatuu(i),i=1,jjm)
679c
680        DO i=1,iip1
681          rlonvv(i)=rlonu(i)*180./pi
682        ENDDO
683       PRINT *,' RLONU '
684       PRINT 400,rlonvv
685c
686       PRINT *,' RLATU '
687        DO i=1,jjp1
688         rlatuu(i)=rlatu(i)*180./pi
689        ENDDO
690       PRINT 400,(rlatuu(i),i=1,jjp1)
691c
692400    FORMAT(1x,8f8.2)
693990    FORMAT(//)
694c
695      RETURN
696      END
Note: See TracBrowser for help on using the repository browser.