/[lmdze]/trunk/libf/dyn3d/inigeom.f90
ViewVC logotype

Diff of /trunk/libf/dyn3d/inigeom.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.3  
changed lines
  Added in v.25

  ViewVC Help
Powered by ViewVC 1.1.21