source: LMDZ5/trunk/libf/dyn3dmem/inigeom.F @ 1663

Last change on this file since 1663 was 1658, checked in by Laurent Fairhead, 12 years ago

Phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ4 (HEAD)
Validation effectuee par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modele donnee par Y. Meurdesoif et la version phasee avec la r1428 (fin du tronc LMDZ4)


Phasing of the localised (low memory) parallel dynamics package with the LMDZ4 trunk version of LMDZ
Validation consisted in comparing output debug files (u, v, t, q, masse, etc... ) of a no physics simulation
run with the version of the code given by Y. Meurdesoif and this version phased with r1428 (HEAD of the LMDZ4 trunk)

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