source: LMDZ5/trunk/libf/dyn3d_common/inigeom.F @ 2306

Last change on this file since 2306 was 2218, checked in by lguez, 10 years ago

Bug fix in fxhyp. There was some bad logic for the shifting of
longitude grids rlonuv, rlonv, rlonm025 and rlonp025 toward [- pi,
pi]. In some cases, one of the four grids was not shifted and the
others were. For example, you could see the bug with: iim = 48, clon =
20, grossismx = 3, dzoomx = 0.15. The bad logic involved the variable
is2 in the loop on ik = 1, 4. The loop included some tests depending
on ik. The whole thing was quite convoluted. Rewrote fxhyp. In
particular, replaced the loop on ik by a call to a new procedure,
invert_zoom_x. fxhyp.F was in dyn3d, it is replaced by fxhyp_m.F90
which is in dyn3d_common. Removed several arguments of fxhyp: zoom
parameters are accessed through serre.h; rlonm025 and rlonp025 were
not used in inigeom; min and max of longitude steps are written inside
fxhyp.

Some simplifications and modernizations in fyhyp. In particular,
several arguments are removed: zoom parameters are accessed through
serre.h; yprimv was not used in inigeom; min and max of latitude steps
are written inside fyhyp.

Removed now useless intermediary procedure fxyhyper.

Changed default value of dzoomx and dzoomy from 0 to 0.2. dzoom[xy] is
only needed when grossism[xy] > 1 and then it should be > 0.

dzoom[xy] can no longer be the extent of the zoom in degrees: it must
be expressed as the fraction of the total domain.

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