/[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 38 by guez, Thu Jan 6 17:52:19 2011 UTC revision 60 by guez, Mon Jan 30 14:37:26 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 dimens_m, ONLY : iim, jjm      USE comconst, ONLY : g, omeg, rad
     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 65  contains Line 58  contains
58           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
59           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
60           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
61        USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
62        use conf_gcm_m, ONLY : fxyhypb, ysinus
63        USE dimens_m, ONLY : iim, jjm
64        use fxy_m, only: fxy
65        use jumble, only: new_unit
66        use nr_util, only: pi
67        USE paramet_m, ONLY : iip1, jjp1
68      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
69           grossismy, pxo, pyo, taux, tauy, transx, transy           grossismy, pxo, pyo, taux, tauy, transx, transy
70    
71      ! Variables locales      ! Variables locales
72    
73      INTEGER i, j, itmax, itmay, iter      INTEGER i, j, itmax, itmay, iter, unit
74      REAL cvu(iip1, jjp1), cuv(iip1, jjm)      REAL cvu(iip1, jjp1), cuv(iip1, jjm)
75      REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm      REAL ai14, ai23, airez, un4rad2
76      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm      REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
77      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
78      REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
79           cuij4(iip1, jjp1)      REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! 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 100  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 119  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 132  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 143  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 159  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 182  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  
184    
185      ! aireij4_2d . . aireij1_2d      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
186        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
187      ! U . . P . U      ! chaque aire_2d(i, j), ainsi que les quatre élongations
188        ! élémentaires cuij et les quatre élongations cvij qui sont
189      ! aireij3_2d . . aireij2_2d      ! calculées aux mêmes endroits que les aireij.
190    
191      ! V      coslatm = cos(rlatu1(1))
192        radclatm = 0.5 * rad * coslatm
193      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,  
194      ! aireij3_2d, aireij4_2d      aireij1_2d(:iim, 1) = 0.
195      ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations      aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
196      ! elementaires      aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
197      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      aireij4_2d(:iim, 1) = 0.
198      ! endroits que les aireij.  
199        cuij1(:iim, 1) = 0.
200      ! do 35 : boucle sur les jjm + 1 latitudes      cuij2(:iim, 1) = radclatm * xprimp025(:iim)
201        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
202      DO j = 1, jjp1      cuij4(:iim, 1) = 0.
203    
204         IF (j==1) THEN      cvij1(:iim, 1) = 0.
205        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
206            yprm = yprimu1(j)      cvij3(:iim, 1) = cvij2(:iim, 1)
207            rlatm = rlatu1(j)      cvij4(:iim, 1) = 0.
208    
209            coslatm = cos(rlatm)      do j = 2, jjm
210            radclatm = 0.5*rad*coslatm         coslatm = cos(rlatu1(j))
211           coslatp = cos(rlatu2(j-1))
212            DO i = 1, iim         radclatp = 0.5 * rad * coslatp
213               xprp = xprimp025(i)         radclatm = 0.5 * rad * coslatm
214               xprm = xprimm025(i)         ai14 = un4rad2 * coslatp * yprimu2(j-1)
215               aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm         ai23 = un4rad2 * coslatm * yprimu1(j)
216               aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm  
217               cuij2(i, 1) = radclatm*xprp         aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
218               cuij3(i, 1) = radclatm*xprm         aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
219               cvij2(i, 1) = 0.5*rad*yprm         aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
220               cvij3(i, 1) = cvij2(i, 1)         aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
221            END DO         cuij1(:iim, j) = radclatp * xprimp025(:iim)
222           cuij2(:iim, j) = radclatm * xprimp025(:iim)
223            DO i = 1, iim         cuij3(:iim, j) = radclatm * xprimm025(:iim)
224               aireij1_2d(i, 1) = 0.         cuij4(:iim, j) = radclatp * xprimm025(:iim)
225               aireij4_2d(i, 1) = 0.         cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
226               cuij1(i, 1) = 0.         cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
227               cuij4(i, 1) = 0.         cvij3(:iim, j) = cvij2(:iim, j)
228               cvij1(i, 1) = 0.         cvij4(:iim, j) = cvij1(:iim, j)
229               cvij4(i, 1) = 0.      end do
230            END DO  
231        coslatp = cos(rlatu2(jjm))
232         END IF      radclatp = 0.5 * rad * coslatp
233    
234         IF (j==jjp1) THEN      aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
235            yprp = yprimu2(j-1)      aireij2_2d(:iim, jjp1) = 0.
236            rlatp = rlatu2(j-1)      aireij3_2d(:iim, jjp1) = 0.
237        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
238            coslatp = cos(rlatp)  
239            radclatp = 0.5*rad*coslatp      cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
240        cuij2(:iim, jjp1) = 0.
241            DO i = 1, iim      cuij3(:iim, jjp1) = 0.
242               xprp = xprimp025(i)      cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
243               xprm = xprimm025(i)  
244               aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp      cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
245               aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp      cvij2(:iim, jjp1) = 0.
246               cuij1(i, jjp1) = radclatp*xprp      cvij3(:iim, jjp1) = 0.
247               cuij4(i, jjp1) = radclatp*xprm      cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
248               cvij1(i, jjp1) = 0.5*rad*yprp  
249               cvij4(i, jjp1) = cvij1(i, jjp1)      ! Périodicité :
250            END DO  
251        cvij1(iip1, :) = cvij1(1, :)
252            DO i = 1, iim      cvij2(iip1, :) = cvij2(1, :)
253               aireij2_2d(i, jjp1) = 0.      cvij3(iip1, :) = cvij3(1, :)
254               aireij3_2d(i, jjp1) = 0.      cvij4(iip1, :) = cvij4(1, :)
255               cvij2(i, jjp1) = 0.  
256               cvij3(i, jjp1) = 0.      cuij1(iip1, :) = cuij1(1, :)
257               cuij2(i, jjp1) = 0.      cuij2(iip1, :) = cuij2(1, :)
258               cuij3(i, jjp1) = 0.      cuij3(iip1, :) = cuij3(1, :)
259            END DO      cuij4(iip1, :) = cuij4(1, :)
260    
261         END IF      aireij1_2d(iip1, :) = aireij1_2d(1, :)
262        aireij2_2d(iip1, :) = aireij2_2d(1, :)
263         IF (j>1 .AND. j<jjp1) THEN      aireij3_2d(iip1, :) = aireij3_2d(1, :)
264        aireij4_2d(iip1, :) = aireij4_2d(1, :)
           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  
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 378  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 392  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 441  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 459  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 473  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 482  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 506  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 522  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.38  
changed lines
  Added in v.60
  ViewVC Help
Powered by ViewVC 1.1.21