source: LMDZ5/trunk/libf/dyn3d/inigeom.F @ 1645

Last change on this file since 1645 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files


Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

  • 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 1403 2010-07-01 09:02:53Z lguez $
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.