/[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

revision 24 by guez, Wed Mar 3 13:23:49 2010 UTC revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC
# Line 1  Line 1 
1  SUBROUTINE inigeom  module inigeom_m
   
   !     Auteur :  P. Le Van  
   
   !   ............      Version  du 01/04/2001     ...................  
   
   !  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-  
   !     endroits que les aires aireij1_2d,..aireij4_2d .  
   
   !  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.  
   ! Possibilité d'appeler une fonction "f(y)" à  
   ! dérivée tangente hyperbolique à la place de la fonction à dérivée  
   ! sinusoïdale.  
   
   
   USE dimens_m, ONLY : iim, jjm  
   USE paramet_m, ONLY : iip1, jjp1  
   USE comconst, ONLY : g, omeg, pi, rad  
   USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh  
   USE logic, ONLY : fxyhypb, ysinus  
   USE comgeom, ONLY : aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d, &  
        airesurg_2d, aireu_2d, airev_2d, aire_2d, airuscv2_2d, airvscu2_2d, &  
        aiuscv2gam_2d, aivscu2gam_2d, alpha1p2_2d, alpha1p4_2d, alpha1_2d, &  
        alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &  
        apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &  
        cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &  
        cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &  
        rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &  
        unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, unsapolsga1, &  
        unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv  
   USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &  
        grossismy, pxo, pyo, taux, tauy, transx, transy  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    !   ....  Variables  locales   ....  contains
   
   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  
   
   !   calcul des coeff. ( cu_2d, cv_2d , 1./cu_2d**2,  1./cv_2d**2  )  
   !   -                                                                -  
   !   ------------------------------------------------------------------  
   
   ! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles  
   !      aux vitesses covariantes et contravariantes , ou vice-versa ...  
   
   
   ! on a :  u (covariant) = cu_2d * u (naturel)   , u(contrav)= u(nat)/cu_2d  
   !         v (covariant) = cv_2d * v (naturel)   , v(contrav)= v(nat)/cv_2d  
   
   !       on en tire :  u(covariant) = cu_2d * cu_2d * u(contravariant)  
   !                     v(covariant) = cv_2d * cv_2d * v(contravariant)  
   
   
   !     on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2  
   !                                                          =     =  
   !                                           et   - jm/2    <  Y  < jm/2  
   !                                                          =     =  
   
   !      ...................................................  
   !      ...................................................  
   !      .  x  est la longitude du point  en radians       .  
   !      .  y  est la  latitude du point  en radians       .  
   !      .                                                 .  
   !      .  on a :  cu_2d(i,j) = rad * COS(y) * dx/dX         .  
   !      .          cv( j ) = rad          * dy/dY         .  
   !      .        aire_2d(i,j) =  cu_2d(i,j) * cv(j)             .  
   !      .                                                 .  
   !      . y, dx/dX, dy/dY calcules aux points concernes   .  
   !      .                                                 .  
   !      ...................................................  
   !      ...................................................  
   
   
   
   !                                                           ,  
   !    cv , bien que dependant de j uniquement,sera ici indice aussi en i  
   !          pour un adressage plus facile en  ij  .  
   
   
   
   !  **************  aux points  u  et  v ,           *****************  
   !      xprimu et xprimv sont respectivement les valeurs de  dx/dX  
   !      yprimu et yprimv    .  .  .  .  .  .  .  .  .  .  .  dy/dY  
   !      rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude  
   !      cvu    et   cv_2d      .  .  .  .  .  .  .  .  .  .  .    cv_2d  
   
   !  **************  aux points u, v, scalaires, et z  ****************  
   !      cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de    cu_2d  
   
   
   
   !         Exemple de distribution de variables sur la grille dans le  
   !             domaine de travail ( X,Y ) .  
   !     ................................................................  
   !                  DX=DY= 1  
   
   
   !        +     represente  un  point scalaire ( p.exp  la pression )  
   !        >     represente  la composante zonale du  vent  
   !        V     represente  la composante meridienne du vent  
   !        o     represente  la  vorticite  
   
   !     ----  , car aux poles , les comp.zonales covariantes sont nulles  
   
   
   
   !         i ->  
   
   !         1      2      3      4      5      6      7      8  
   !  j  
   !  v  1   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --  
   
   !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
   !     2   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
   !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
   !     3   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
   !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
   !     4   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
   !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
   !     5   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --  
   
   
   !      Ci-dessus,  on voit que le nombre de pts.en longitude est egal  
   !                 a   IM = 8  
   !      De meme ,   le nombre d'intervalles entre les 2 poles est egal  
   !                 a   JM = 4  
   
   !      Les points scalaires ( + ) correspondent donc a des valeurs  
   !       entieres  de  i ( 1 a IM )   et  de  j ( 1 a  JM +1 )   .  
   
   !      Les vents    U       ( > ) correspondent a des valeurs  semi-  
   !       entieres  de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)  
   
   !      Les vents    V       ( V ) correspondent a des valeurs entieres  
   !     de     i ( 1 a  IM ) et semi-entieres de  j ( 1 +0.5  a JM +0.5)  
   
   
   
   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   . '/)  
   
   
   IF (nitergdiv/=2) THEN  
      gamdi_gdiv = coefdis/(float(nitergdiv)-2.)  
   ELSE  
      gamdi_gdiv = 0.  
   END IF  
   IF (nitergrot/=2) THEN  
      gamdi_grot = coefdis/(float(nitergrot)-2.)  
   ELSE  
      gamdi_grot = 0.  
   END IF  
   IF (niterh/=2) THEN  
      gamdi_h = coefdis/(float(niterh)-2.)  
   ELSE  
      gamdi_h = 0.  
   END IF  
   
   WRITE (6,*) ' gamdi_gd ', gamdi_gdiv, gamdi_grot, gamdi_h, coefdis, &  
        nitergdiv, nitergrot, niterh  
   
   pi = 2.*asin(1.)  
   
   WRITE (6,990)  
   
   !     ----------------------------------------------------------------  
   
   IF ( .NOT. fxyhypb) THEN  
   
   
      IF (ysinus) THEN  
   
         WRITE (6,*) ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '  
   
         !   .... 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)  
   
      ELSE  
   
         WRITE (6,*) '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'  
   
         ! utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ..  
   
   
         pxo = clon*pi/180.  
         pyo = 2.*clat*pi/180.  
   
         !  ....  determination de  transx ( pour le zoom ) par Newton-Raphson .  
   
         itmax = 10  
         eps = .1E-7  
   
         xo1 = 0.  
         DO 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<=eps) EXIT  
            xo1 = x1  
         END DO  
   
         transx = xo1  
   
         itmay = 10  
         eps = .1E-7  
   
         yo1 = 0.  
         DO 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<=eps) EXIT  
            yo1 = y1  
         END DO  
   
         transy = yo1  
   
         CALL fxy(rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,rlatu2,yprimu2, &  
              rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)  
   
      END IF  
   
   ELSE  
   
      ! ....  Utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.  
      !   ..................................................................  
   
      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)  
   
   
   END IF  
   
   !  -------------------------------------------------------------------  
   
   
   rlatu(1) = asin(1.)  
   rlatu(jjp1) = -rlatu(1)  
   
6    
7    !   ....  calcul  aux  poles  ....    SUBROUTINE inigeom
8    
9    yprimu(1) = 0.      ! Auteur : P. Le Van
10    yprimu(jjp1) = 0.      ! 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
181              print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
182              ! 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    un4rad2 = 0.25*rad*rad  400 FORMAT (1X, 8F8.2)
   
   !   -------------------------------------------------------------  
   !   -------------------------------------------------------------  
   !                                                                    -  
   ! calcul  des aires ( aire_2d,aireu_2d,airev_2d, 1./aire_2d, 1./airez )  
   !   -      et de   fext_2d ,  force de coriolis  extensive  .  
   !   -  
   !   -------------------------------------------------------------  
   !   -------------------------------------------------------------  
   
   
   
   !   A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont  
   !   affectees 4 aires entourant P , calculees respectivement aux points  
   !            ( i + 1/4, j - 1/4 )    :    aireij1_2d (i,j)  
   !            ( i + 1/4, j + 1/4 )    :    aireij2_2d (i,j)  
   !            ( i - 1/4, j + 1/4 )    :    aireij3_2d (i,j)  
   !            ( i - 1/4, j - 1/4 )    :    aireij4_2d (i,j)  
   
   !           ,  
   ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).  
   !   Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme  
   ! des 4 aires  aireij1_2d,aireij2_2d,aireij3_2d,aireij4_2d qui sont affectees au  
   !   point (i,j) .  
   !   On definit en outre les coefficients  alpha comme etant egaux a  
   ! (aireij / aire_2d), c.a.d par exp.  alpha1_2d(i,j)=aireij1_2d(i,j)/aire_2d(i,j)  
   
   !   De meme, toute aire centree en 1 point U est egale a la somme des  
   ! 4 aires aireij1_2d,aireij2_2d,aireij3_2d,aireij4_2d entourant le point U.  
   !    Idem pour  airev_2d, airez .  
   
   !       On a ,pour chaque maille :    dX = dY = 1  
   
   
   !                             . V  
   
   !                 aireij4_2d .        . aireij1_2d  
   
   !                   U .       . P      . U  
   
   !                 aireij3_2d .        . aireij2_2d  
   
   !                             . V  
   
   
   
   
   
   ! ....................................................................  
   
   ! Calcul des 4 aires elementaires aireij1_2d,aireij2_2d,aireij3_2d,aireij4_2d  
   ! qui entourent chaque aire_2d(i,j) , ainsi que les 4 elongations elementaires  
   ! cuij et les 4 elongat. cvij qui sont calculees aux memes  
   !     endroits  que les aireij   .  
   
   !  ....................................................................  
   
   !     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....  
   
   
   DO j = 1, jjp1  
   
      IF (j==1) THEN  
   
         yprm = yprimu1(j)  
         rlatm = rlatu1(j)  
   
         coslatm = cos(rlatm)  
         radclatm = 0.5*rad*coslatm  
   
         DO 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)  
         END DO  
   
         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.  
         END DO  
   
      END IF  
   
      IF (j==jjp1) THEN  
         yprp = yprimu2(j-1)  
         rlatp = rlatu2(j-1)  
         !cc       yprp             = fyprim( FLOAT(j) - 0.25 )  
         !cc       rlatp            = fy    ( FLOAT(j) - 0.25 )  
   
         coslatp = cos(rlatp)  
         radclatp = 0.5*rad*coslatp  
   
         DO 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)  
         END DO  
   
         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.  
         END DO  
   
      END IF  
   
   
      IF (j>1 .AND. j<jjp1) THEN  
   
         rlatp = rlatu2(j-1)  
         yprp = yprimu2(j-1)  
         rlatm = rlatu1(j)  
         yprm = yprimu1(j)  
         !c         rlatp    = fy    ( FLOAT(j) - 0.25 )  
         !c         yprp     = fyprim( FLOAT(j) - 0.25 )  
         !c         rlatm    = fy    ( FLOAT(j) + 0.25 )  
         !c         yprm     = fyprim( FLOAT(j) + 0.25 )  
   
         coslatm = cos(rlatm)  
         coslatp = cos(rlatp)  
         radclatp = 0.5*rad*coslatp  
         radclatm = 0.5*rad*coslatm  
   
         DO 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)  
         END DO  
   
      END IF  
   
      !    ........       periodicite   ............  
   
      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)  
   
   END DO  
   
   !    ..............................................................  
   
   DO j = 1, jjp1  
      DO 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)  
      END DO  
   
   
      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)  
   END DO  
   
   
   DO j = 1, jjp1  
      DO 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  
      END DO  
      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)  
   END DO  
   
   
   DO j = 1, jjm  
   
      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)  
      END DO  
      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  
      END DO  
      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)  
   
   END DO  
   
   
   !    .....      Calcul  des elongations cu_2d,cv_2d, cvu     .........  
   
   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))  
      END DO  
      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)  
      END DO  
      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)  
   END DO  
   
   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)  
      END DO  
      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)  
   END DO  
   
   
   !   ....  calcul aux  poles  ....  
   
   DO i = 1, iip1  
      cu_2d(i,1) = 0.  
      unscu2_2d(i,1) = 0.  
      cvu(i,1) = 0.  
   
      cu_2d(i,jjp1) = 0.  
      unscu2_2d(i,jjp1) = 0.  
      cvu(i,jjp1) = 0.  
   END DO  
   
   !    ..............................................................  
   
   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)  
      END DO  
      airvscu2_2d(iip1,j) = airvscu2_2d(1,j)  
      aivscu2gam_2d(iip1,j) = aivscu2gam_2d(1,j)  
   END DO  
   
   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)  
      END DO  
      airuscv2_2d(iip1,j) = airuscv2_2d(1,j)  
      aiuscv2gam_2d(iip1,j) = aiuscv2gam_2d(1,j)  
   END DO  
   
   
   !   calcul des aires aux  poles :  
   !   -----------------------------  
   
   apoln = sum(aire_2d(:iim, 1))  
   apols = sum(aire_2d(:iim, jjp1))  
   unsapolnga1 = 1./(apoln**(-gamdi_gdiv))  
   unsapolsga1 = 1./(apols**(-gamdi_gdiv))  
   unsapolnga2 = 1./(apoln**(-gamdi_h))  
   unsapolsga2 = 1./(apols**(-gamdi_h))  
   
   !----------------------------------------------------------------  
   !     gtitre='Coriolis version ancienne'  
   !     gfichier='fext1'  
   !     CALL writestd(fext_2d,iip1*jjm)  
   
   !   changement F. Hourdin calcul conservatif pour fext_2d  
   !   constang_2d contient le produit a * cos ( latitude ) * omega  
   
   DO i = 1, iim  
      constang_2d(i,1) = 0.  
   END DO  
   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))  
      END DO  
   END DO  
   DO i = 1, iim  
      constang_2d(i,jjp1) = 0.  
   END DO  
   
   !   periodicite en longitude  
   
   DO j = 1, jjm  
      fext_2d(iip1,j) = fext_2d(1,j)  
   END DO  
   DO j = 1, jjp1  
      constang_2d(iip1,j) = constang_2d(1,j)  
   END DO  
   
   ! fin du changement  
   
   
   !----------------------------------------------------------------  
   
   WRITE (6,*) '   ***  Coordonnees de la grille  *** '  
   WRITE (6,995)  
   
   WRITE (6,*) '   LONGITUDES  aux pts.   V  ( degres )  '  
   WRITE (6,995)  
   DO i = 1, iip1  
      rlonvv(i) = rlonv(i)*180./pi  
   END DO  
   WRITE (6,400) rlonvv  
   
   WRITE (6,995)  
   WRITE (6,*) '   LATITUDES   aux pts.   V  ( degres )  '  
   WRITE (6,995)  
   DO i = 1, jjm  
      rlatuu(i) = rlatv(i)*180./pi  
   END DO  
   WRITE (6,400) (rlatuu(i),i=1,jjm)  
   
   DO i = 1, iip1  
      rlonvv(i) = rlonu(i)*180./pi  
   END DO  
   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  
   END DO  
   WRITE (6,400) (rlatuu(i),i=1,jjp1)  
   WRITE (6,995)  
   
 400 FORMAT (1X,8F8.2)  
612  990 FORMAT (//)  990 FORMAT (//)
613  995 FORMAT (/)  995 FORMAT (/)
614    
615  END SUBROUTINE inigeom    END SUBROUTINE inigeom
616    
617    end module inigeom_m

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

  ViewVC Help
Powered by ViewVC 1.1.21