/[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 61 by guez, Fri Apr 20 14:58:43 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 ! in m
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  
184    
185      ! V      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
186        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
187      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      ! chaque aire_2d(i, j), ainsi que les quatre élongations
188      ! aireij3_2d, aireij4_2d      ! élémentaires cuij et les quatre élongations cvij qui sont
189      ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations      ! calculées aux mêmes endroits que les aireij.
190      ! elementaires  
191      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      coslatm = cos(rlatu1(1))
192      ! endroits que les aireij.      radclatm = 0.5 * rad * coslatm
193    
194      ! do 35 : boucle sur les jjm + 1 latitudes      aireij1_2d(:iim, 1) = 0.
195        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
196      DO j = 1, jjp1      aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
197        aireij4_2d(:iim, 1) = 0.
198         IF (j==1) THEN  
199        cuij1(:iim, 1) = 0.
200            yprm = yprimu1(j)      cuij2(:iim, 1) = radclatm * xprimp025(:iim)
201            rlatm = rlatu1(j)      cuij3(:iim, 1) = radclatm * xprimm025(:iim)
202        cuij4(:iim, 1) = 0.
203            coslatm = cos(rlatm)  
204            radclatm = 0.5*rad*coslatm      cvij1(:iim, 1) = 0.
205        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
206            DO i = 1, iim      cvij3(:iim, 1) = cvij2(:iim, 1)
207               xprp = xprimp025(i)      cvij4(:iim, 1) = 0.
208               xprm = xprimm025(i)  
209               aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm      do j = 2, jjm
210               aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm         coslatm = cos(rlatu1(j))
211               cuij2(i, 1) = radclatm*xprp         coslatp = cos(rlatu2(j-1))
212               cuij3(i, 1) = radclatm*xprm         radclatp = 0.5 * rad * coslatp
213               cvij2(i, 1) = 0.5*rad*yprm         radclatm = 0.5 * rad * coslatm
214               cvij3(i, 1) = cvij2(i, 1)         ai14 = un4rad2 * coslatp * yprimu2(j-1)
215            END DO         ai23 = un4rad2 * coslatm * yprimu1(j)
216    
217            DO i = 1, iim         aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
218               aireij1_2d(i, 1) = 0.         aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
219               aireij4_2d(i, 1) = 0.         aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
220               cuij1(i, 1) = 0.         aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
221               cuij4(i, 1) = 0.         cuij1(:iim, j) = radclatp * xprimp025(:iim)
222               cvij1(i, 1) = 0.         cuij2(:iim, j) = radclatm * xprimp025(:iim)
223               cvij4(i, 1) = 0.         cuij3(:iim, j) = radclatm * xprimm025(:iim)
224            END DO         cuij4(:iim, j) = radclatp * xprimm025(:iim)
225           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
226         END IF         cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
227           cvij3(:iim, j) = cvij2(:iim, j)
228         IF (j==jjp1) THEN         cvij4(:iim, j) = cvij1(:iim, j)
229            yprp = yprimu2(j-1)      end do
230            rlatp = rlatu2(j-1)  
231        coslatp = cos(rlatu2(jjm))
232            coslatp = cos(rlatp)      radclatp = 0.5 * rad * coslatp
233            radclatp = 0.5*rad*coslatp  
234        aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
235            DO i = 1, iim      aireij2_2d(:iim, jjp1) = 0.
236               xprp = xprimp025(i)      aireij3_2d(:iim, jjp1) = 0.
237               xprm = xprimm025(i)      aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
238               aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp  
239               aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp      cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
240               cuij1(i, jjp1) = radclatp*xprp      cuij2(:iim, jjp1) = 0.
241               cuij4(i, jjp1) = radclatp*xprm      cuij3(:iim, jjp1) = 0.
242               cvij1(i, jjp1) = 0.5*rad*yprp      cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
243               cvij4(i, jjp1) = cvij1(i, jjp1)  
244            END DO      cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
245        cvij2(:iim, jjp1) = 0.
246            DO i = 1, iim      cvij3(:iim, jjp1) = 0.
247               aireij2_2d(i, jjp1) = 0.      cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
248               aireij3_2d(i, jjp1) = 0.  
249               cvij2(i, jjp1) = 0.      ! Périodicité :
250               cvij3(i, jjp1) = 0.  
251               cuij2(i, jjp1) = 0.      cvij1(iip1, :) = cvij1(1, :)
252               cuij3(i, jjp1) = 0.      cvij2(iip1, :) = cvij2(1, :)
253            END DO      cvij3(iip1, :) = cvij3(1, :)
254        cvij4(iip1, :) = cvij4(1, :)
255         END IF  
256        cuij1(iip1, :) = cuij1(1, :)
257         IF (j>1 .AND. j<jjp1) THEN      cuij2(iip1, :) = cuij2(1, :)
258        cuij3(iip1, :) = cuij3(1, :)
259            rlatp = rlatu2(j-1)      cuij4(iip1, :) = cuij4(1, :)
260            yprp = yprimu2(j-1)  
261            rlatm = rlatu1(j)      aireij1_2d(iip1, :) = aireij1_2d(1, :)
262            yprm = yprimu1(j)      aireij2_2d(iip1, :) = aireij2_2d(1, :)
263        aireij3_2d(iip1, :) = aireij3_2d(1, :)
264            coslatm = cos(rlatm)      aireij4_2d(iip1, :) = aireij4_2d(1, :)
           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 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)**2
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)**2
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      cu_2d(:, 1) = 0.
378         cu_2d(i, 1) = 0.      unscu2_2d(:, 1) = 0.
379         unscu2_2d(i, 1) = 0.      cvu(:, 1) = 0.
380         cvu(i, 1) = 0.  
381        cu_2d(:, jjp1) = 0.
382         cu_2d(i, jjp1) = 0.      unscu2_2d(:, jjp1) = 0.
383         unscu2_2d(i, jjp1) = 0.      cvu(:, jjp1) = 0.
        cvu(i, jjp1) = 0.  
     END DO  
384    
385      DO j = 1, jjm      DO j = 1, jjm
386         DO i = 1, iim         DO i = 1, iim
387            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))
388            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)            aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
389         END DO         END DO
390         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)         airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
# Line 483  contains Line 393  contains
393    
394      DO j = 2, jjm      DO j = 2, jjm
395         DO i = 1, iim         DO i = 1, iim
396            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))
397            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)            aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
398         END DO         END DO
399         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)         airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
400         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
401      END DO      END DO
402    
403      ! calcul des aires aux poles :      ! Calcul des aires aux pôles :
404    
405      apoln = sum(aire_2d(:iim, 1))      apoln = sum(aire_2d(:iim, 1))
406      apols = sum(aire_2d(:iim, jjp1))      apols = sum(aire_2d(:iim, jjp1))
407      unsapolnga1 = 1./(apoln**(-gamdi_gdiv))      unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
408      unsapolsga1 = 1./(apols**(-gamdi_gdiv))      unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
409      unsapolnga2 = 1./(apoln**(-gamdi_h))      unsapolnga2 = 1. / (apoln**(-gamdi_h))
410      unsapolsga2 = 1./(apols**(-gamdi_h))      unsapolsga2 = 1. / (apols**(-gamdi_h))
411    
412      ! changement F. Hourdin calcul conservatif pour fext_2d      ! Changement F. Hourdin calcul conservatif pour fext_2d
413      ! constang_2d contient le produit a * cos (latitude) * omega      ! constang_2d contient le produit a * cos (latitude) * omega
414    
415      DO i = 1, iim      DO i = 1, iim
# Line 507  contains Line 417  contains
417      END DO      END DO
418      DO j = 1, jjm - 1      DO j = 1, jjm - 1
419         DO i = 1, iim         DO i = 1, iim
420            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) &
421                   * cos(rlatu(j + 1))
422         END DO         END DO
423      END DO      END DO
424      DO i = 1, iim      DO i = 1, iim
425         constang_2d(i, jjp1) = 0.         constang_2d(i, jjp1) = 0.
426      END DO      END DO
427    
428      ! periodicite en longitude      ! Périodicité en longitude
429    
430      DO j = 1, jjm      DO j = 1, jjm
431         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
# Line 523  contains Line 434  contains
434         constang_2d(iip1, j) = constang_2d(1, j)         constang_2d(iip1, j) = constang_2d(1, j)
435      END DO      END DO
436    
437      ! fin du changement      call new_unit(unit)
438        open(unit, file="longitude_latitude.txt", status="replace", action="write")
439      print *, ' Coordonnees de la grille '      write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
440      WRITE (6, 995)      write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
441        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
442      print *, ' LONGITUDES aux pts. V (degres) '      write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
443      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 (/)  
444    
445    END SUBROUTINE inigeom    END SUBROUTINE inigeom
446    
Legend:
Removed from v.39  
changed lines
  Added in v.61
  ViewVC Help
Powered by ViewVC 1.1.21