/[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 25 by guez, Fri Mar 5 16:43:45 2010 UTC revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC
# Line 7  contains Line 7  contains
7    SUBROUTINE inigeom    SUBROUTINE inigeom
8    
9      ! Auteur : P. Le Van      ! Auteur : P. Le Van
     ! Version du 01/04/2001  
10    
11      ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes      ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
12      ! endroits que les aires aireij1_2d, ..., aireij4_2d.      ! endroits que les aires aireij1_2d, ..., aireij4_2d.
13    
14      ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée      ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à
15      ! tangente hyperbolique      ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,
16      ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2)      ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d
17        ! permettent de passer des vitesses naturelles aux vitesses
18        ! covariantes et contravariantes, ou vice-versa.
19    
20        ! On a :
21        ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
22        ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
23    
24        ! On en tire :
25        ! 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        ! et - jm / 2 <= Y <= jm / 2
30    
31        ! x est la longitude du point en radians.
32        ! y est la latitude du point en radians.
33        !
34        ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
35        ! cv(j) = rad * dy / dY
36        ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
37        !
38        ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
39        ! dépendant de j uniquement, sera ici indicé aussi en i pour un
40        ! adressage plus facile en ij.
41    
42        ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
43        ! points u et v. yprimu et yprimv sont respectivement les valeurs
44        ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
45        ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
46        ! respectivement les valeurs de cv_2d aux points u et v.
47    
48      ! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
49      !      aux vitesses covariantes et contravariantes , ou vice-versa ...      ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
50    
51      ! on a :  u (covariant) = cu_2d * u (naturel)   , u(contrav)= u(nat)/cu_2d      USE comconst, ONLY : g, omeg, rad
     !         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)  
   
     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  
52      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
53           alpha1p2_2d, alpha1p4_2d, alpha1_2d, &           alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
54           alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &           alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
# Line 112  contains Line 58  contains
58           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
59           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
60           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
61        USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
62        use conf_gcm_m, ONLY : fxyhypb, ysinus
63        USE dimens_m, ONLY : iim, jjm
64        use fxy_m, only: fxy
65        use jumble, only: new_unit
66        use nr_util, only: pi
67        USE paramet_m, ONLY : iip1, jjp1
68      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
69           grossismy, pxo, pyo, taux, tauy, transx, transy           grossismy, pxo, pyo, taux, tauy, transx, transy
70    
71      ! Variables locales      ! Variables locales
72    
73      INTEGER i, j, itmax, itmay, iter      INTEGER i, j, itmax, itmay, iter, unit
74      REAL cvu(iip1, jjp1), cuv(iip1, jjm)      REAL cvu(iip1, jjp1), cuv(iip1, jjm)
75      REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm      REAL ai14, ai23, airez, un4rad2
76      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
77      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
78      REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
79           cuij4(iip1, jjp1)      REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4
80      REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &      REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
81           cvij4(iip1, jjp1)      real yprimv(jjm), yprimu(jjp1)
     REAL rlonvv(iip1), rlatuu(jjp1)  
     REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), &  
          yprimu(jjp1)  
82      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
   
83      REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)      REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
84      SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
85      SAVE rlonm025, xprimm025, rlonp025, xprimp025           aireij4_2d ! in m2
   
     real aireij1_2d(iim + 1, jjm + 1)  
     real aireij2_2d(iim + 1, jjm + 1)  
     real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1)  
86      real airuscv2_2d(iim + 1, jjm)      real airuscv2_2d(iim + 1, jjm)
87      real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
88      real aivscu2gam_2d(iim + 1, jjm)      real aivscu2gam_2d(iim + 1, jjm)
89    
90      !------------------------------------------------------------------      !------------------------------------------------------------------
# Line 147  contains Line 92  contains
92      PRINT *, 'Call sequence information: inigeom'      PRINT *, 'Call sequence information: inigeom'
93    
94      IF (nitergdiv/=2) THEN      IF (nitergdiv/=2) THEN
95         gamdi_gdiv = coefdis/(real(nitergdiv)-2.)         gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
96      ELSE      ELSE
97         gamdi_gdiv = 0.         gamdi_gdiv = 0.
98      END IF      END IF
99      IF (nitergrot/=2) THEN      IF (nitergrot/=2) THEN
100         gamdi_grot = coefdis/(real(nitergrot)-2.)         gamdi_grot = coefdis / (real(nitergrot)-2.)
101      ELSE      ELSE
102         gamdi_grot = 0.         gamdi_grot = 0.
103      END IF      END IF
104      IF (niterh/=2) THEN      IF (niterh/=2) THEN
105         gamdi_h = coefdis/(real(niterh)-2.)         gamdi_h = coefdis / (real(niterh)-2.)
106      ELSE      ELSE
107         gamdi_h = 0.         gamdi_h = 0.
108      END IF      END IF
# Line 166  contains Line 111  contains
111      print *, "gamdi_grot = ", gamdi_grot      print *, "gamdi_grot = ", gamdi_grot
112      print *, "gamdi_h = ", gamdi_h      print *, "gamdi_h = ", gamdi_h
113    
114      WRITE (6, 990)      IF (.NOT. fxyhypb) THEN
   
     IF ( .NOT. fxyhypb) THEN  
115         IF (ysinus) THEN         IF (ysinus) THEN
116            print *, ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '            print *, ' Inigeom, Y = Sinus (Latitude) '
117              ! utilisation de f(x, y) avec y = sinus de la latitude
           ! utilisation de f(x, y )  avec  y  =  sinus de la latitude  ...  
   
118            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
119                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
120                 xprimm025, rlonp025, xprimp025)                 xprimm025, rlonp025, xprimp025)
121         ELSE         ELSE
122            print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'            print *, 'Inigeom, Y = Latitude, der. sinusoid .'
123            ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit            ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
124    
125            pxo = clon*pi/180.            pxo = clon * pi / 180.
126            pyo = 2.*clat*pi/180.            pyo = 2. * clat * pi / 180.
127    
128            ! determination de  transx ( pour le zoom ) par Newton-Raphson .            ! determination de transx (pour le zoom) par Newton-Raphson
129    
130            itmax = 10            itmax = 10
131            eps = .1E-7            eps = .1E-7
# Line 192  contains Line 133  contains
133            xo1 = 0.            xo1 = 0.
134            DO iter = 1, itmax            DO iter = 1, itmax
135               x1 = xo1               x1 = xo1
136               f = x1 + alphax*sin(x1-pxo)               f = x1 + alphax * sin(x1-pxo)
137               df = 1. + alphax*cos(x1-pxo)               df = 1. + alphax * cos(x1-pxo)
138               x1 = x1 - f/df               x1 = x1 - f / df
139               xdm = abs(x1-xo1)               xdm = abs(x1-xo1)
140               IF (xdm<=eps) EXIT               IF (xdm<=eps) EXIT
141               xo1 = x1               xo1 = x1
# Line 208  contains Line 149  contains
149            yo1 = 0.            yo1 = 0.
150            DO iter = 1, itmay            DO iter = 1, itmay
151               y1 = yo1               y1 = yo1
152               f = y1 + alphay*sin(y1-pyo)               f = y1 + alphay * sin(y1-pyo)
153               df = 1. + alphay*cos(y1-pyo)               df = 1. + alphay * cos(y1-pyo)
154               y1 = y1 - f/df               y1 = y1 - f / df
155               ydm = abs(y1-yo1)               ydm = abs(y1-yo1)
156               IF (ydm<=eps) EXIT               IF (ydm<=eps) EXIT
157               yo1 = y1               yo1 = y1
# Line 223  contains Line 164  contains
164                 rlonp025, xprimp025)                 rlonp025, xprimp025)
165         END IF         END IF
166      ELSE      ELSE
167         ! ....  Utilisation  de fxyhyper , f(x, y) a derivee tangente hyperbol.         ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
168         print *, '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'         print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
169         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
170              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
171              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
172              rlonp025, xprimp025)              rlonp025, xprimp025)
173      END IF      END IF
174    
175      rlatu(1) = asin(1.)      rlatu(1) = pi / 2.
176      rlatu(jjp1) = -rlatu(1)      rlatu(jjp1) = -rlatu(1)
177    
178      !   ....  calcul  aux  poles  ....      ! Calcul aux pôles
179    
180      yprimu(1) = 0.      yprimu(1) = 0.
181      yprimu(jjp1) = 0.      yprimu(jjp1) = 0.
182    
183      un4rad2 = 0.25*rad*rad      un4rad2 = 0.25 * rad * rad
   
     ! 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  
184    
185      !                   U .       . P      . U      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
186        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
187      !                 aireij3_2d .        . aireij2_2d      ! chaque aire_2d(i, j), ainsi que les quatre élongations
188        ! élémentaires cuij et les quatre élongations cvij qui sont
189      !                             . V      ! calculées aux mêmes endroits que les aireij.
190    
191      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      coslatm = cos(rlatu1(1))
192      ! aireij3_2d, aireij4_2d      radclatm = 0.5 * rad * coslatm
193      ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations  
194      ! elementaires      aireij1_2d(:iim, 1) = 0.
195      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
196      !     endroits  que les aireij   .      aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
197        aireij4_2d(:iim, 1) = 0.
198      !     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....  
199        cuij1(:iim, 1) = 0.
200      DO j = 1, jjp1      cuij2(:iim, 1) = radclatm * xprimp025(:iim)
201        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
202         IF (j==1) THEN      cuij4(:iim, 1) = 0.
203    
204            yprm = yprimu1(j)      cvij1(:iim, 1) = 0.
205            rlatm = rlatu1(j)      cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
206        cvij3(:iim, 1) = cvij2(:iim, 1)
207            coslatm = cos(rlatm)      cvij4(:iim, 1) = 0.
208            radclatm = 0.5*rad*coslatm  
209        do j = 2, jjm
210            DO i = 1, iim         coslatm = cos(rlatu1(j))
211               xprp = xprimp025(i)         coslatp = cos(rlatu2(j-1))
212               xprm = xprimm025(i)         radclatp = 0.5 * rad * coslatp
213               aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm         radclatm = 0.5 * rad * coslatm
214               aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm         ai14 = un4rad2 * coslatp * yprimu2(j-1)
215               cuij2(i, 1) = radclatm*xprp         ai23 = un4rad2 * coslatm * yprimu1(j)
216               cuij3(i, 1) = radclatm*xprm  
217               cvij2(i, 1) = 0.5*rad*yprm         aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
218               cvij3(i, 1) = cvij2(i, 1)         aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
219            END DO         aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
220           aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
221            DO i = 1, iim         cuij1(:iim, j) = radclatp * xprimp025(:iim)
222               aireij1_2d(i, 1) = 0.         cuij2(:iim, j) = radclatm * xprimp025(:iim)
223               aireij4_2d(i, 1) = 0.         cuij3(:iim, j) = radclatm * xprimm025(:iim)
224               cuij1(i, 1) = 0.         cuij4(:iim, j) = radclatp * xprimm025(:iim)
225               cuij4(i, 1) = 0.         cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
226               cvij1(i, 1) = 0.         cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
227               cvij4(i, 1) = 0.         cvij3(:iim, j) = cvij2(:iim, j)
228            END DO         cvij4(:iim, j) = cvij1(:iim, j)
229        end do
230         END IF  
231        coslatp = cos(rlatu2(jjm))
232         IF (j==jjp1) THEN      radclatp = 0.5 * rad * coslatp
233            yprp = yprimu2(j-1)  
234            rlatp = rlatu2(j-1)      aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
235        aireij2_2d(:iim, jjp1) = 0.
236            coslatp = cos(rlatp)      aireij3_2d(:iim, jjp1) = 0.
237            radclatp = 0.5*rad*coslatp      aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
238    
239            DO i = 1, iim      cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
240               xprp = xprimp025(i)      cuij2(:iim, jjp1) = 0.
241               xprm = xprimm025(i)      cuij3(:iim, jjp1) = 0.
242               aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp      cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
243               aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp  
244               cuij1(i, jjp1) = radclatp*xprp      cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
245               cuij4(i, jjp1) = radclatp*xprm      cvij2(:iim, jjp1) = 0.
246               cvij1(i, jjp1) = 0.5*rad*yprp      cvij3(:iim, jjp1) = 0.
247               cvij4(i, jjp1) = cvij1(i, jjp1)      cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
248            END DO  
249        ! Périodicité :
250            DO i = 1, iim  
251               aireij2_2d(i, jjp1) = 0.      cvij1(iip1, :) = cvij1(1, :)
252               aireij3_2d(i, jjp1) = 0.      cvij2(iip1, :) = cvij2(1, :)
253               cvij2(i, jjp1) = 0.      cvij3(iip1, :) = cvij3(1, :)
254               cvij3(i, jjp1) = 0.      cvij4(iip1, :) = cvij4(1, :)
255               cuij2(i, jjp1) = 0.  
256               cuij3(i, jjp1) = 0.      cuij1(iip1, :) = cuij1(1, :)
257            END DO      cuij2(iip1, :) = cuij2(1, :)
258        cuij3(iip1, :) = cuij3(1, :)
259         END IF      cuij4(iip1, :) = cuij4(1, :)
260    
261         IF (j>1 .AND. j<jjp1) THEN      aireij1_2d(iip1, :) = aireij1_2d(1, :)
262        aireij2_2d(iip1, :) = aireij2_2d(1, :)
263            rlatp = rlatu2(j-1)      aireij3_2d(iip1, :) = aireij3_2d(1, :)
264            yprp = yprimu2(j-1)      aireij4_2d(iip1, :) = aireij4_2d(1, :)
           rlatm = rlatu1(j)  
           yprm = yprimu1(j)  
   
           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  
265    
266      DO j = 1, jjp1      DO j = 1, jjp1
267         DO i = 1, iim         DO i = 1, iim
268            aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &            aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
269                 + aireij3_2d(i, j) + aireij4_2d(i, j)                 + aireij3_2d(i, j) + aireij4_2d(i, j)
270            alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)            alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
271            alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)            alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
272            alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)            alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
273            alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)            alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
274            alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)            alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
275            alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)            alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
276            alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)            alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
# Line 427  contains Line 291  contains
291      DO j = 1, jjp1      DO j = 1, jjp1
292         DO i = 1, iim         DO i = 1, iim
293            aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &            aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
294                 aireij4_2d(i+1, j) + aireij3_2d(i+1, j)                 aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
295            unsaire_2d(i, j) = 1./aire_2d(i, j)            unsaire_2d(i, j) = 1. / aire_2d(i, j)
296            unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)            unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
297            unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)            unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
298            airesurg_2d(i, j) = aire_2d(i, j)/g            airesurg_2d(i, j) = aire_2d(i, j) / g
299         END DO         END DO
300         aireu_2d(iip1, j) = aireu_2d(1, j)         aireu_2d(iip1, j) = aireu_2d(1, j)
301         unsaire_2d(iip1, j) = unsaire_2d(1, j)         unsaire_2d(iip1, j) = unsaire_2d(1, j)
# Line 441  contains Line 305  contains
305      END DO      END DO
306    
307      DO j = 1, jjm      DO j = 1, jjm
   
308         DO i = 1, iim         DO i = 1, iim
309            airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &            airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
310                 aireij1_2d(i, j+1) + aireij4_2d(i, j+1)                 aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
311         END DO         END DO
312         DO i = 1, iim         DO i = 1, iim
313            airez = aireij2_2d(i, j) + aireij1_2d(i, j+1) + aireij3_2d(i+1, j) &            airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
314                 + aireij4_2d(i+1, j+1)                 + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
315            unsairez_2d(i, j) = 1./airez            unsairez_2d(i, j) = 1. / airez
316            unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)            unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
317            fext_2d(i, j) = airez*sin(rlatv(j))*2.*omeg            fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
318         END DO         END DO
319         airev_2d(iip1, j) = airev_2d(1, j)         airev_2d(iip1, j) = airev_2d(1, j)
320         unsairez_2d(iip1, j) = unsairez_2d(1, j)         unsairez_2d(iip1, j) = unsairez_2d(1, j)
321         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
322         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
   
323      END DO      END DO
324    
325      !    .....      Calcul  des elongations cu_2d, cv_2d, cvu     .........      ! Calcul des élongations cu_2d, cv_2d, cvu
326    
327      DO j = 1, jjm      DO j = 1, jjm
328         DO i = 1, iim         DO i = 1, iim
329            cv_2d(i, j) = 0.5 * &            cv_2d(i, j) = 0.5 * &
330                 (cvij2(i, j) + cvij3(i, j) + cvij1(i, j+1) + cvij4(i, j+1))                 (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
331            cvu(i, j) = 0.5*(cvij1(i, j)+cvij4(i, j)+cvij2(i, j)+cvij3(i, j))            cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
332            cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))                 + cvij3(i, j))
333            unscv2_2d(i, j) = 1./(cv_2d(i, j)*cv_2d(i, j))            cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
334                   + cuij4(i, j + 1))
335              unscv2_2d(i, j) = 1. / (cv_2d(i, j) * cv_2d(i, j))
336         END DO         END DO
337         DO i = 1, iim         DO i = 1, iim
338            cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)            cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
339            cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)            cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
340            cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)            cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
341            cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)            cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
342            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
# Line 490  contains Line 354  contains
354    
355      DO j = 2, jjm      DO j = 2, jjm
356         DO i = 1, iim         DO i = 1, iim
357            cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i+1, j) + cuij2(i, j) &            cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
358                 + cuij3(i+1, j))                 + cuij3(i + 1, j))
359            unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))            unscu2_2d(i, j) = 1. / (cu_2d(i, j) * cu_2d(i, j))
360            cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)            cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
361            cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)            cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
362            cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)            cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
363            cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)            cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
364            cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)            cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
# Line 508  contains Line 372  contains
372         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
373      END DO      END DO
374    
375      !   ....  calcul aux  poles  ....      ! Calcul aux pôles
376    
377      DO i = 1, iip1      DO i = 1, iip1
378         cu_2d(i, 1) = 0.         cu_2d(i, 1) = 0.
# Line 522  contains Line 386  contains
386    
387      DO j = 1, jjm      DO j = 1, jjm
388         DO i = 1, iim         DO i = 1, iim
389            airvscu2_2d(i, j) = airev_2d(i, j)/(cuv(i, j)*cuv(i, j))            airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
390            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
391         END DO         END DO
392         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
# Line 531  contains Line 395  contains
395    
396      DO j = 2, jjm      DO j = 2, jjm
397         DO i = 1, iim         DO i = 1, iim
398            airuscv2_2d(i, j) = aireu_2d(i, j)/(cvu(i, j)*cvu(i, j))            airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
399            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
400         END DO         END DO
401         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
402         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
403      END DO      END DO
404    
405      !   calcul des aires aux  poles :      ! Calcul des aires aux pôles :
406    
407      apoln = sum(aire_2d(:iim, 1))      apoln = sum(aire_2d(:iim, 1))
408      apols = sum(aire_2d(:iim, jjp1))      apols = sum(aire_2d(:iim, jjp1))
409      unsapolnga1 = 1./(apoln**(-gamdi_gdiv))      unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
410      unsapolsga1 = 1./(apols**(-gamdi_gdiv))      unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
411      unsapolnga2 = 1./(apoln**(-gamdi_h))      unsapolnga2 = 1. / (apoln**(-gamdi_h))
412      unsapolsga2 = 1./(apols**(-gamdi_h))      unsapolsga2 = 1. / (apols**(-gamdi_h))
413    
414      !   changement F. Hourdin calcul conservatif pour fext_2d      ! Changement F. Hourdin calcul conservatif pour fext_2d
415      !   constang_2d contient le produit a * cos ( latitude ) * omega      ! constang_2d contient le produit a * cos (latitude) * omega
416    
417      DO i = 1, iim      DO i = 1, iim
418         constang_2d(i, 1) = 0.         constang_2d(i, 1) = 0.
419      END DO      END DO
420      DO j = 1, jjm - 1      DO j = 1, jjm - 1
421         DO i = 1, iim         DO i = 1, iim
422            constang_2d(i, j+1) = rad*omeg*cu_2d(i, j+1)*cos(rlatu(j+1))            constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
423                   * cos(rlatu(j + 1))
424         END DO         END DO
425      END DO      END DO
426      DO i = 1, iim      DO i = 1, iim
427         constang_2d(i, jjp1) = 0.         constang_2d(i, jjp1) = 0.
428      END DO      END DO
429    
430      !   periodicite en longitude      ! Périodicité en longitude
431    
432      DO j = 1, jjm      DO j = 1, jjm
433         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
# Line 571  contains Line 436  contains
436         constang_2d(iip1, j) = constang_2d(1, j)         constang_2d(iip1, j) = constang_2d(1, j)
437      END DO      END DO
438    
439      ! fin du changement      call new_unit(unit)
440        open(unit, file="longitude_latitude.txt", status="replace", action="write")
441      print *, '   ***  Coordonnees de la grille  *** '      write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
442      WRITE (6, 995)      write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
443        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
444      print *, '   LONGITUDES  aux pts.   V  ( degres )  '      write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
445      WRITE (6, 995)      close(unit)
     DO i = 1, iip1  
        rlonvv(i) = rlonv(i)*180./pi  
     END DO  
     WRITE (6, 400) rlonvv  
   
     WRITE (6, 995)  
     print *, '   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)  
     print *, '   LONGITUDES  aux pts.   U  ( degres )  '  
     WRITE (6, 995)  
     WRITE (6, 400) rlonvv  
     WRITE (6, 995)  
   
     print *, '   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)  
 990 FORMAT (//)  
 995 FORMAT (/)  
446    
447    END SUBROUTINE inigeom    END SUBROUTINE inigeom
448    
Legend:
Removed from v.25  
changed lines
  Added in v.57
  ViewVC Help
Powered by ViewVC 1.1.21