source: trunk/LMDZ.GENERIC/libf/dyn3d/inigeom.F @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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