/[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 70 by guez, Mon Jun 24 15:39:52 2013 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 fxyhyper_m, only: fxyhyper
66        use jumble, only: new_unit
67        use nr_util, only: pi
68        USE paramet_m, ONLY : iip1, jjp1
69      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
70           grossismy, pxo, pyo, taux, tauy, transx, transy           grossismy, pxo, pyo, taux, tauy, transx, transy
71    
72      ! Variables locales      ! Variables locales
73    
74      INTEGER i, j, itmax, itmay, iter      INTEGER i, j, itmax, itmay, iter, unit
75      REAL cvu(iip1, jjp1), cuv(iip1, jjm)      REAL cvu(iip1, jjp1), cuv(iip1, jjm)
76      REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm      REAL ai14, ai23, airez, un4rad2
77      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
78      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
79      REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
80           cuij4(iip1, jjp1)      REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
81      REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &      REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
82           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)  
83      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
   
84      REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)      REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
85      SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
86      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)  
87      real airuscv2_2d(iim + 1, jjm)      real airuscv2_2d(iim + 1, jjm)
88      real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
89      real aivscu2gam_2d(iim + 1, jjm)      real aivscu2gam_2d(iim + 1, jjm)
90    
91      !------------------------------------------------------------------      !------------------------------------------------------------------
# Line 147  contains Line 93  contains
93      PRINT *, 'Call sequence information: inigeom'      PRINT *, 'Call sequence information: inigeom'
94    
95      IF (nitergdiv/=2) THEN      IF (nitergdiv/=2) THEN
96         gamdi_gdiv = coefdis/(real(nitergdiv)-2.)         gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
97      ELSE      ELSE
98         gamdi_gdiv = 0.         gamdi_gdiv = 0.
99      END IF      END IF
100      IF (nitergrot/=2) THEN      IF (nitergrot/=2) THEN
101         gamdi_grot = coefdis/(real(nitergrot)-2.)         gamdi_grot = coefdis / (real(nitergrot)-2.)
102      ELSE      ELSE
103         gamdi_grot = 0.         gamdi_grot = 0.
104      END IF      END IF
105      IF (niterh/=2) THEN      IF (niterh/=2) THEN
106         gamdi_h = coefdis/(real(niterh)-2.)         gamdi_h = coefdis / (real(niterh)-2.)
107      ELSE      ELSE
108         gamdi_h = 0.         gamdi_h = 0.
109      END IF      END IF
# Line 166  contains Line 112  contains
112      print *, "gamdi_grot = ", gamdi_grot      print *, "gamdi_grot = ", gamdi_grot
113      print *, "gamdi_h = ", gamdi_h      print *, "gamdi_h = ", gamdi_h
114    
115      WRITE (6, 990)      IF (.NOT. fxyhypb) THEN
   
     IF ( .NOT. fxyhypb) THEN  
116         IF (ysinus) THEN         IF (ysinus) THEN
117            print *, ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '            print *, ' Inigeom, Y = Sinus (Latitude) '
118              ! utilisation de f(x, y) avec y = sinus de la latitude
           ! utilisation de f(x, y )  avec  y  =  sinus de la latitude  ...  
   
119            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
120                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
121                 xprimm025, rlonp025, xprimp025)                 xprimm025, rlonp025, xprimp025)
122         ELSE         ELSE
123            print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'            print *, 'Inigeom, Y = Latitude, der. sinusoid .'
124            ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit            ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
125    
126            pxo = clon*pi/180.            pxo = clon * pi / 180.
127            pyo = 2.*clat*pi/180.            pyo = 2. * clat * pi / 180.
128    
129            ! determination de  transx ( pour le zoom ) par Newton-Raphson .            ! determination de transx (pour le zoom) par Newton-Raphson
130    
131            itmax = 10            itmax = 10
132            eps = .1E-7            eps = .1E-7
# Line 192  contains Line 134  contains
134            xo1 = 0.            xo1 = 0.
135            DO iter = 1, itmax            DO iter = 1, itmax
136               x1 = xo1               x1 = xo1
137               f = x1 + alphax*sin(x1-pxo)               f = x1 + alphax * sin(x1-pxo)
138               df = 1. + alphax*cos(x1-pxo)               df = 1. + alphax * cos(x1-pxo)
139               x1 = x1 - f/df               x1 = x1 - f / df
140               xdm = abs(x1-xo1)               xdm = abs(x1-xo1)
141               IF (xdm<=eps) EXIT               IF (xdm<=eps) EXIT
142               xo1 = x1               xo1 = x1
# Line 208  contains Line 150  contains
150            yo1 = 0.            yo1 = 0.
151            DO iter = 1, itmay            DO iter = 1, itmay
152               y1 = yo1               y1 = yo1
153               f = y1 + alphay*sin(y1-pyo)               f = y1 + alphay * sin(y1-pyo)
154               df = 1. + alphay*cos(y1-pyo)               df = 1. + alphay * cos(y1-pyo)
155               y1 = y1 - f/df               y1 = y1 - f / df
156               ydm = abs(y1-yo1)               ydm = abs(y1-yo1)
157               IF (ydm<=eps) EXIT               IF (ydm<=eps) EXIT
158               yo1 = y1               yo1 = y1
# Line 223  contains Line 165  contains
165                 rlonp025, xprimp025)                 rlonp025, xprimp025)
166         END IF         END IF
167      ELSE      ELSE
168         ! ....  Utilisation  de fxyhyper , f(x, y) a derivee tangente hyperbol.         ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
169         print *, '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'         print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
170         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
171              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
172              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
173              rlonp025, xprimp025)              rlonp025, xprimp025)
174      END IF      END IF
175    
176      rlatu(1) = asin(1.)      rlatu(1) = pi / 2.
177      rlatu(jjp1) = -rlatu(1)      rlatu(jjp1) = -rlatu(1)
178    
179      !   ....  calcul  aux  poles  ....      ! Calcul aux pôles
180    
181      yprimu(1) = 0.      yprimu(1) = 0.
182      yprimu(jjp1) = 0.      yprimu(jjp1) = 0.
183    
184      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 .  
185    
186      !       On a , pour chaque maille :    dX = dY = 1      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
187        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
188      !                             . V      ! chaque aire_2d(i, j), ainsi que les quatre élongations
189        ! élémentaires cuij et les quatre élongations cvij qui sont
190      !                 aireij4_2d .        . aireij1_2d      ! calculées aux mêmes endroits que les aireij.
191    
192      !                   U .       . P      . U      coslatm = cos(rlatu1(1))
193        radclatm = 0.5 * rad * coslatm
194      !                 aireij3_2d .        . aireij2_2d  
195        aireij1_2d(:iim, 1) = 0.
196      !                             . V      aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
197        aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
198      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      aireij4_2d(:iim, 1) = 0.
199      ! aireij3_2d, aireij4_2d  
200      ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations      cuij1(:iim, 1) = 0.
201      ! elementaires      cuij2(:iim, 1) = radclatm * xprimp025(:iim)
202      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      cuij3(:iim, 1) = radclatm * xprimm025(:iim)
203      !     endroits  que les aireij   .      cuij4(:iim, 1) = 0.
204    
205      !     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....      cvij1(:iim, 1) = 0.
206        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
207      DO j = 1, jjp1      cvij3(:iim, 1) = cvij2(:iim, 1)
208        cvij4(:iim, 1) = 0.
209         IF (j==1) THEN  
210        do j = 2, jjm
211            yprm = yprimu1(j)         coslatm = cos(rlatu1(j))
212            rlatm = rlatu1(j)         coslatp = cos(rlatu2(j-1))
213           radclatp = 0.5 * rad * coslatp
214            coslatm = cos(rlatm)         radclatm = 0.5 * rad * coslatm
215            radclatm = 0.5*rad*coslatm         ai14 = un4rad2 * coslatp * yprimu2(j-1)
216           ai23 = un4rad2 * coslatm * yprimu1(j)
217            DO i = 1, iim  
218               xprp = xprimp025(i)         aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
219               xprm = xprimm025(i)         aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
220               aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm         aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
221               aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm         aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
222               cuij2(i, 1) = radclatm*xprp         cuij1(:iim, j) = radclatp * xprimp025(:iim)
223               cuij3(i, 1) = radclatm*xprm         cuij2(:iim, j) = radclatm * xprimp025(:iim)
224               cvij2(i, 1) = 0.5*rad*yprm         cuij3(:iim, j) = radclatm * xprimm025(:iim)
225               cvij3(i, 1) = cvij2(i, 1)         cuij4(:iim, j) = radclatp * xprimm025(:iim)
226            END DO         cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
227           cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
228            DO i = 1, iim         cvij3(:iim, j) = cvij2(:iim, j)
229               aireij1_2d(i, 1) = 0.         cvij4(:iim, j) = cvij1(:iim, j)
230               aireij4_2d(i, 1) = 0.      end do
231               cuij1(i, 1) = 0.  
232               cuij4(i, 1) = 0.      coslatp = cos(rlatu2(jjm))
233               cvij1(i, 1) = 0.      radclatp = 0.5 * rad * coslatp
234               cvij4(i, 1) = 0.  
235            END DO      aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
236        aireij2_2d(:iim, jjp1) = 0.
237         END IF      aireij3_2d(:iim, jjp1) = 0.
238        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
239         IF (j==jjp1) THEN  
240            yprp = yprimu2(j-1)      cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
241            rlatp = rlatu2(j-1)      cuij2(:iim, jjp1) = 0.
242        cuij3(:iim, jjp1) = 0.
243            coslatp = cos(rlatp)      cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
244            radclatp = 0.5*rad*coslatp  
245        cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
246            DO i = 1, iim      cvij2(:iim, jjp1) = 0.
247               xprp = xprimp025(i)      cvij3(:iim, jjp1) = 0.
248               xprm = xprimm025(i)      cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
249               aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp  
250               aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp      ! Périodicité :
251               cuij1(i, jjp1) = radclatp*xprp  
252               cuij4(i, jjp1) = radclatp*xprm      cvij1(iip1, :) = cvij1(1, :)
253               cvij1(i, jjp1) = 0.5*rad*yprp      cvij2(iip1, :) = cvij2(1, :)
254               cvij4(i, jjp1) = cvij1(i, jjp1)      cvij3(iip1, :) = cvij3(1, :)
255            END DO      cvij4(iip1, :) = cvij4(1, :)
256    
257            DO i = 1, iim      cuij1(iip1, :) = cuij1(1, :)
258               aireij2_2d(i, jjp1) = 0.      cuij2(iip1, :) = cuij2(1, :)
259               aireij3_2d(i, jjp1) = 0.      cuij3(iip1, :) = cuij3(1, :)
260               cvij2(i, jjp1) = 0.      cuij4(iip1, :) = cuij4(1, :)
261               cvij3(i, jjp1) = 0.  
262               cuij2(i, jjp1) = 0.      aireij1_2d(iip1, :) = aireij1_2d(1, :)
263               cuij3(i, jjp1) = 0.      aireij2_2d(iip1, :) = aireij2_2d(1, :)
264            END DO      aireij3_2d(iip1, :) = aireij3_2d(1, :)
265        aireij4_2d(iip1, :) = aireij4_2d(1, :)
        END IF  
   
        IF (j>1 .AND. j<jjp1) THEN  
   
           rlatp = rlatu2(j-1)  
           yprp = yprimu2(j-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  
266    
267      DO j = 1, jjp1      DO j = 1, jjp1
268         DO i = 1, iim         DO i = 1, iim
269            aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &            aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
270                 + aireij3_2d(i, j) + aireij4_2d(i, j)                 + aireij3_2d(i, j) + aireij4_2d(i, j)
271            alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)            alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
272            alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)            alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
273            alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)            alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
274            alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)            alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
275            alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)            alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
276            alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)            alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
277            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 292  contains
292      DO j = 1, jjp1      DO j = 1, jjp1
293         DO i = 1, iim         DO i = 1, iim
294            aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &            aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
295                 aireij4_2d(i+1, j) + aireij3_2d(i+1, j)                 aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
296            unsaire_2d(i, j) = 1./aire_2d(i, j)            unsaire_2d(i, j) = 1. / aire_2d(i, j)
297            unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)            unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
298            unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)            unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
299            airesurg_2d(i, j) = aire_2d(i, j)/g            airesurg_2d(i, j) = aire_2d(i, j) / g
300         END DO         END DO
301         aireu_2d(iip1, j) = aireu_2d(1, j)         aireu_2d(iip1, j) = aireu_2d(1, j)
302         unsaire_2d(iip1, j) = unsaire_2d(1, j)         unsaire_2d(iip1, j) = unsaire_2d(1, j)
# Line 441  contains Line 306  contains
306      END DO      END DO
307    
308      DO j = 1, jjm      DO j = 1, jjm
   
309         DO i = 1, iim         DO i = 1, iim
310            airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &            airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
311                 aireij1_2d(i, j+1) + aireij4_2d(i, j+1)                 aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
312         END DO         END DO
313         DO i = 1, iim         DO i = 1, iim
314            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) &
315                 + aireij4_2d(i+1, j+1)                 + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
316            unsairez_2d(i, j) = 1./airez            unsairez_2d(i, j) = 1. / airez
317            unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)            unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
318            fext_2d(i, j) = airez*sin(rlatv(j))*2.*omeg            fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
319         END DO         END DO
320         airev_2d(iip1, j) = airev_2d(1, j)         airev_2d(iip1, j) = airev_2d(1, j)
321         unsairez_2d(iip1, j) = unsairez_2d(1, j)         unsairez_2d(iip1, j) = unsairez_2d(1, j)
322         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
323         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
   
324      END DO      END DO
325    
326      !    .....      Calcul  des elongations cu_2d, cv_2d, cvu     .........      ! Calcul des élongations cu_2d, cv_2d, cvu
327    
328      DO j = 1, jjm      DO j = 1, jjm
329         DO i = 1, iim         DO i = 1, iim
330            cv_2d(i, j) = 0.5 * &            cv_2d(i, j) = 0.5 * &
331                 (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))
332            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) &
333            cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))                 + cvij3(i, j))
334            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) &
335                   + cuij4(i, j + 1))
336              unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
337         END DO         END DO
338         DO i = 1, iim         DO i = 1, iim
339            cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)            cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
340            cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)            cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
341            cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)            cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
342            cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)            cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
343            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
# Line 490  contains Line 355  contains
355    
356      DO j = 2, jjm      DO j = 2, jjm
357         DO i = 1, iim         DO i = 1, iim
358            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) &
359                 + cuij3(i+1, j))                 + cuij3(i + 1, j))
360            unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))            unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
361            cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)            cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
362            cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)            cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
363            cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)            cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
364            cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)            cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
365            cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)            cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
# Line 508  contains Line 373  contains
373         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
374      END DO      END DO
375    
376      !   ....  calcul aux  poles  ....      ! Calcul aux pôles
377    
378      DO i = 1, iip1      cu_2d(:, 1) = 0.
379         cu_2d(i, 1) = 0.      unscu2_2d(:, 1) = 0.
380         unscu2_2d(i, 1) = 0.      cvu(:, 1) = 0.
381         cvu(i, 1) = 0.  
382        cu_2d(:, jjp1) = 0.
383         cu_2d(i, jjp1) = 0.      unscu2_2d(:, jjp1) = 0.
384         unscu2_2d(i, jjp1) = 0.      cvu(:, jjp1) = 0.
        cvu(i, jjp1) = 0.  
     END DO  
385    
386      DO j = 1, jjm      DO j = 1, jjm
387         DO i = 1, iim         DO i = 1, iim
388            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))
389            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
390         END DO         END DO
391         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
# Line 531  contains Line 394  contains
394    
395      DO j = 2, jjm      DO j = 2, jjm
396         DO i = 1, iim         DO i = 1, iim
397            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))
398            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
399         END DO         END DO
400         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
401         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
402      END DO      END DO
403    
404      !   calcul des aires aux  poles :      ! Calcul des aires aux pôles :
405    
406      apoln = sum(aire_2d(:iim, 1))      apoln = sum(aire_2d(:iim, 1))
407      apols = sum(aire_2d(:iim, jjp1))      apols = sum(aire_2d(:iim, jjp1))
408      unsapolnga1 = 1./(apoln**(-gamdi_gdiv))      unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
409      unsapolsga1 = 1./(apols**(-gamdi_gdiv))      unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
410      unsapolnga2 = 1./(apoln**(-gamdi_h))      unsapolnga2 = 1. / (apoln**(-gamdi_h))
411      unsapolsga2 = 1./(apols**(-gamdi_h))      unsapolsga2 = 1. / (apols**(-gamdi_h))
412    
413      !   changement F. Hourdin calcul conservatif pour fext_2d      ! Changement F. Hourdin calcul conservatif pour fext_2d
414      !   constang_2d contient le produit a * cos ( latitude ) * omega      ! constang_2d contient le produit a * cos (latitude) * omega
415    
416      DO i = 1, iim      DO i = 1, iim
417         constang_2d(i, 1) = 0.         constang_2d(i, 1) = 0.
418      END DO      END DO
419      DO j = 1, jjm - 1      DO j = 1, jjm - 1
420         DO i = 1, iim         DO i = 1, iim
421            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) &
422                   * cos(rlatu(j + 1))
423         END DO         END DO
424      END DO      END DO
425      DO i = 1, iim      DO i = 1, iim
426         constang_2d(i, jjp1) = 0.         constang_2d(i, jjp1) = 0.
427      END DO      END DO
428    
429      !   periodicite en longitude      ! Périodicité en longitude
430    
431      DO j = 1, jjm      DO j = 1, jjm
432         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
# Line 571  contains Line 435  contains
435         constang_2d(iip1, j) = constang_2d(1, j)         constang_2d(iip1, j) = constang_2d(1, j)
436      END DO      END DO
437    
438      ! fin du changement      call new_unit(unit)
439        open(unit, file="longitude_latitude.txt", status="replace", action="write")
440      print *, '   ***  Coordonnees de la grille  *** '      write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
441      WRITE (6, 995)      write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
442        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
443      print *, '   LONGITUDES  aux pts.   V  ( degres )  '      write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
444      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 (/)  
445    
446    END SUBROUTINE inigeom    END SUBROUTINE inigeom
447    
Legend:
Removed from v.25  
changed lines
  Added in v.70
  ViewVC Help
Powered by ViewVC 1.1.21