/[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 39 by guez, Tue Jan 25 15:11:05 2011 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      ! les coef. (cu_2d, cv_2d) permettent de passer des vitesses naturelles      ! covariantes et contravariantes, ou vice-versa.
19      ! aux vitesses covariantes et contravariantes, ou vice-versa  
20        ! On a :
21      ! on a :      ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
22      ! u (covariant) = cu_2d * u (naturel), u(contrav)= u(nat)/cu_2d      ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
     ! v (covariant) = cv_2d * v (naturel), v(contrav)= v(nat)/cv_2d  
23    
24      ! on en tire :      ! On en tire :
25      ! u(covariant) = cu_2d * cu_2d * u(contravariant)      ! u(covariant) = cu_2d * cu_2d * u(contravariant)
26      ! v(covariant) = cv_2d * cv_2d * v(contravariant)      ! 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      ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
29      ! et - jm/2 <= Y <= jm/2      ! et - jm / 2 <= Y <= jm / 2
30    
31      ! x est la longitude du point en radians.      ! x est la longitude du point en radians.
32      ! y est la latitude du point en radians.      ! y est la latitude du point en radians.
33      !      !
34      ! on a : cu_2d(i, j) = rad * cos(y) * dx/dX      ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
35      ! cv(j) = rad * dy/dY      ! cv(j) = rad * dy / dY
36      ! aire_2d(i, j) = cu_2d(i, j) * cv(j)      ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
37      !      !
38      ! y, dx/dX, dy/dY calcules aux points concernes      ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
39      ! cv, bien que dependant de j uniquement, sera ici indice aussi en i      ! dépendant de j uniquement, sera ici indicé aussi en i pour un
40      ! pour un adressage plus facile en ij.      ! adressage plus facile en ij.
41    
42      ! aux points u et v,      ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
43      ! xprimu et xprimv sont respectivement les valeurs de dx/dX      ! points u et v. yprimu et yprimv sont respectivement les valeurs
44      ! yprimu et yprimv sont respectivement les valeurs de dy/dY      ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
45      ! rlatu et rlatv sont respectivement les valeurs de la latitude      ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
46      ! cvu et cv_2d sont respectivement les valeurs de cv_2d      ! respectivement les valeurs de cv_2d aux points u et v.
47    
     ! aux points u, v, scalaires, et z  
48      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
49      ! Cf. "inigeom.txt".      ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
50    
51      USE comconst, ONLY : g, omeg, rad      USE comconst, ONLY : g, omeg, rad
52      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
# Line 62  contains Line 59  contains
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      USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
62        use conf_gcm_m, ONLY : fxyhypb, ysinus
63      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
64      USE logic, ONLY : fxyhypb, ysinus      use fxy_m, only: fxy
65        use jumble, only: new_unit
66      use nr_util, only: pi      use nr_util, only: pi
67      USE paramet_m, ONLY : iip1, jjp1      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, &
# Line 71  contains Line 70  contains
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)
# Line 101  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 120  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    
     WRITE (6, 990)  
   
114      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) '
# Line 133  contains Line 122  contains
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    
# Line 144  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 160  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 183  contains Line 172  contains
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  
   
     ! U . . P . U  
   
     ! aireij3_2d . . aireij2_2d  
   
     ! V  
184    
185      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
186      ! aireij3_2d, aireij4_2d      ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
187      ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations      ! chaque aire_2d(i, j), ainsi que les quatre élongations
188      ! elementaires      ! élémentaires cuij et les quatre élongations cvij qui sont
189      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      ! calculées aux mêmes endroits que les aireij.
190      ! endroits que les aireij.  
191        coslatm = cos(rlatu1(1))
192      ! do 35 : boucle sur les jjm + 1 latitudes      radclatm = 0.5 * rad * coslatm
193    
194      DO j = 1, jjp1      aireij1_2d(:iim, 1) = 0.
195        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
196         IF (j==1) THEN      aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
197        aireij4_2d(:iim, 1) = 0.
198            yprm = yprimu1(j)  
199            rlatm = rlatu1(j)      cuij1(:iim, 1) = 0.
200        cuij2(:iim, 1) = radclatm * xprimp025(:iim)
201            coslatm = cos(rlatm)      cuij3(:iim, 1) = radclatm * xprimm025(:iim)
202            radclatm = 0.5*rad*coslatm      cuij4(:iim, 1) = 0.
203    
204            DO i = 1, iim      cvij1(:iim, 1) = 0.
205               xprp = xprimp025(i)      cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
206               xprm = xprimm025(i)      cvij3(:iim, 1) = cvij2(:iim, 1)
207               aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm      cvij4(:iim, 1) = 0.
208               aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm  
209               cuij2(i, 1) = radclatm*xprp      do j = 2, jjm
210               cuij3(i, 1) = radclatm*xprm         coslatm = cos(rlatu1(j))
211               cvij2(i, 1) = 0.5*rad*yprm         coslatp = cos(rlatu2(j-1))
212               cvij3(i, 1) = cvij2(i, 1)         radclatp = 0.5 * rad * coslatp
213            END DO         radclatm = 0.5 * rad * coslatm
214           ai14 = un4rad2 * coslatp * yprimu2(j-1)
215            DO i = 1, iim         ai23 = un4rad2 * coslatm * yprimu1(j)
216               aireij1_2d(i, 1) = 0.  
217               aireij4_2d(i, 1) = 0.         aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
218               cuij1(i, 1) = 0.         aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
219               cuij4(i, 1) = 0.         aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
220               cvij1(i, 1) = 0.         aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
221               cvij4(i, 1) = 0.         cuij1(:iim, j) = radclatp * xprimp025(:iim)
222            END DO         cuij2(:iim, j) = radclatm * xprimp025(:iim)
223           cuij3(:iim, j) = radclatm * xprimm025(:iim)
224         END IF         cuij4(:iim, j) = radclatp * xprimm025(:iim)
225           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
226         IF (j==jjp1) THEN         cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
227            yprp = yprimu2(j-1)         cvij3(:iim, j) = cvij2(:iim, j)
228            rlatp = rlatu2(j-1)         cvij4(:iim, j) = cvij1(:iim, j)
229        end do
230            coslatp = cos(rlatp)  
231            radclatp = 0.5*rad*coslatp      coslatp = cos(rlatu2(jjm))
232        radclatp = 0.5 * rad * coslatp
233            DO i = 1, iim  
234               xprp = xprimp025(i)      aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
235               xprm = xprimm025(i)      aireij2_2d(:iim, jjp1) = 0.
236               aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp      aireij3_2d(:iim, jjp1) = 0.
237               aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp      aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
238               cuij1(i, jjp1) = radclatp*xprp  
239               cuij4(i, jjp1) = radclatp*xprm      cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
240               cvij1(i, jjp1) = 0.5*rad*yprp      cuij2(:iim, jjp1) = 0.
241               cvij4(i, jjp1) = cvij1(i, jjp1)      cuij3(:iim, jjp1) = 0.
242            END DO      cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
243    
244            DO i = 1, iim      cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
245               aireij2_2d(i, jjp1) = 0.      cvij2(:iim, jjp1) = 0.
246               aireij3_2d(i, jjp1) = 0.      cvij3(:iim, jjp1) = 0.
247               cvij2(i, jjp1) = 0.      cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
248               cvij3(i, jjp1) = 0.  
249               cuij2(i, jjp1) = 0.      ! Périodicité :
250               cuij3(i, jjp1) = 0.  
251            END DO      cvij1(iip1, :) = cvij1(1, :)
252        cvij2(iip1, :) = cvij2(1, :)
253         END IF      cvij3(iip1, :) = cvij3(1, :)
254        cvij4(iip1, :) = cvij4(1, :)
255         IF (j>1 .AND. j<jjp1) THEN  
256        cuij1(iip1, :) = cuij1(1, :)
257            rlatp = rlatu2(j-1)      cuij2(iip1, :) = cuij2(1, :)
258            yprp = yprimu2(j-1)      cuij3(iip1, :) = cuij3(1, :)
259            rlatm = rlatu1(j)      cuij4(iip1, :) = cuij4(1, :)
260            yprm = yprimu1(j)  
261        aireij1_2d(iip1, :) = aireij1_2d(1, :)
262            coslatm = cos(rlatm)      aireij2_2d(iip1, :) = aireij2_2d(1, :)
263            coslatp = cos(rlatp)      aireij3_2d(iip1, :) = aireij3_2d(1, :)
264            radclatp = 0.5*rad*coslatp      aireij4_2d(iip1, :) = aireij4_2d(1, :)
           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 379  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 393  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 442  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 460  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 474  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 483  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
# Line 507  contains Line 419  contains
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 523  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.39  
changed lines
  Added in v.57
  ViewVC Help
Powered by ViewVC 1.1.21