/[lmdze]/trunk/dyn3d/comgeom.f
ViewVC logotype

Diff of /trunk/dyn3d/comgeom.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/comgeom.f90 revision 60 by guez, Mon Jan 30 14:37:26 2012 UTC trunk/dyn3d/comgeom.f revision 112 by guez, Thu Sep 18 13:36:51 2014 UTC
# Line 1  Line 1 
1  module comgeom  module comgeom
2    
3    use dimens_m, only: iim, jjm    use dimens_m, only: iim, jjm
   use paramet_m, only: ip1jmp1, ip1jm  
4    
5    implicit none    implicit none
6    
7    private iim, jjm, ip1jmp1, ip1jm    private iim, jjm
8    
9    real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m    real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
10    real cu(ip1jmp1), cv(ip1jm) ! in m    real cu((iim + 1) * (jjm + 1)), cv((iim + 1) * jjm) ! in m
11    equivalence (cu, cu_2d), (cv, cv_2d)    equivalence (cu, cu_2d), (cv, cv_2d)
12    
13    real unscu2_2d(iim + 1, jjm + 1) ! in m-2    real unscu2_2d(iim + 1, jjm + 1) ! in m-2
14    real unscu2(ip1jmp1) ! in m-2    real unscu2((iim + 1) * (jjm + 1)) ! in m-2
15    equivalence (unscu2, unscu2_2d)    equivalence (unscu2, unscu2_2d)
16    
17    real unscv2_2d(iim + 1, jjm) ! in m-2    real unscv2_2d(iim + 1, jjm) ! in m-2
18    real unscv2(ip1jm) ! in m-2    real unscv2((iim + 1) * jjm) ! in m-2
19    equivalence (unscv2, unscv2_2d)    equivalence (unscv2, unscv2_2d)
20    
21    real aire(ip1jmp1), aire_2d(iim + 1, jjm + 1) ! in m2    real aire((iim + 1) * (jjm + 1)), aire_2d(iim + 1, jjm + 1) ! in m2
22    real airesurg_2d(iim + 1, jjm + 1), airesurg(ip1jmp1)    real airesurg_2d(iim + 1, jjm + 1), airesurg((iim + 1) * (jjm + 1))
23    equivalence (aire, aire_2d), (airesurg, airesurg_2d)    equivalence (aire, aire_2d), (airesurg, airesurg_2d)
24    
25    real aireu_2d(iim + 1, jjm + 1) ! in m2    real aireu_2d(iim + 1, jjm + 1) ! in m2
26    real aireu(ip1jmp1) ! in m2    real aireu((iim + 1) * (jjm + 1)) ! in m2
27    equivalence (aireu, aireu_2d)    equivalence (aireu, aireu_2d)
28    
29    real airev(ip1jm), airev_2d(iim + 1, jjm) ! in m2    real airev((iim + 1) * jjm), airev_2d(iim + 1, jjm) ! in m2
30    real unsaire(ip1jmp1), unsaire_2d(iim + 1, jjm + 1) ! in m-2    real unsaire((iim + 1) * (jjm + 1)), unsaire_2d(iim + 1, jjm + 1) ! in m-2
31    equivalence (airev, airev_2d), (unsaire, unsaire_2d)    equivalence (airev, airev_2d), (unsaire, unsaire_2d)
32    
33    real apoln, apols ! in m2    real apoln, apols ! in m2
34    
35    real unsairez_2d(iim + 1, jjm)    real unsairez_2d(iim + 1, jjm)
36    real unsairez(ip1jm)    real unsairez((iim + 1) * jjm)
37    equivalence (unsairez, unsairez_2d)    equivalence (unsairez, unsairez_2d)
38    
39    real alpha1_2d(iim + 1, jjm + 1)    real alpha1_2d(iim + 1, jjm + 1)
40    real alpha1(ip1jmp1)    real alpha1((iim + 1) * (jjm + 1))
41    equivalence (alpha1, alpha1_2d)    equivalence (alpha1, alpha1_2d)
42    
43    real alpha2_2d(iim + 1, jjm + 1)            real alpha2_2d(iim + 1, jjm + 1)        
44    real alpha2(ip1jmp1)    real alpha2((iim + 1) * (jjm + 1))
45    equivalence (alpha2, alpha2_2d)    equivalence (alpha2, alpha2_2d)
46    
47    real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)    real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
48    real alpha3(ip1jmp1), alpha4(ip1jmp1)    real alpha3((iim + 1) * (jjm + 1)), alpha4((iim + 1) * (jjm + 1))
49    equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)    equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)
50    
51    real alpha1p2_2d(iim + 1, jjm + 1)            real alpha1p2_2d(iim + 1, jjm + 1)        
52    real alpha1p2(ip1jmp1)    real alpha1p2((iim + 1) * (jjm + 1))
53    equivalence (alpha1p2, alpha1p2_2d)    equivalence (alpha1p2, alpha1p2_2d)
54    
55    real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)    real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
56    real alpha1p4(ip1jmp1), alpha2p3(ip1jmp1)    real alpha1p4((iim + 1) * (jjm + 1)), alpha2p3((iim + 1) * (jjm + 1))
57    equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)    equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)
58    
59    real alpha3p4(ip1jmp1)    real alpha3p4((iim + 1) * (jjm + 1))
60    real alpha3p4_2d(iim + 1, jjm + 1)        real alpha3p4_2d(iim + 1, jjm + 1)    
61    equivalence (alpha3p4, alpha3p4_2d)    equivalence (alpha3p4, alpha3p4_2d)
62    
63    real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)    real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
64    real fext(ip1jm), constang(ip1jmp1)    real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
65    equivalence (fext, fext_2d), (constang, constang_2d)    equivalence (fext, fext_2d), (constang, constang_2d)
66    
67    real rlatu(jjm + 1)    real rlatu(jjm + 1)
# Line 77  module comgeom Line 76  module comgeom
76    ! (longitudes of points of the "scalar" and "v" grid, in rad)    ! (longitudes of points of the "scalar" and "v" grid, in rad)
77    
78    real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension    real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
79    real cuvsurcv(ip1jm), cvsurcuv(ip1jm) ! no dimension    real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
80    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
81    
82    real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)    real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
83    ! no dimension    ! no dimension
84    real cvusurcu(ip1jmp1), cusurcvu(ip1jmp1) ! no dimension    real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
85      ! no dimension
86    equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)    equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
87    
88    real cuvscvgam1_2d(iim + 1, jjm)    real cuvscvgam1_2d(iim + 1, jjm)
89    real cuvscvgam1(ip1jm)    real cuvscvgam1((iim + 1) * jjm)
90    equivalence (cuvscvgam1, cuvscvgam1_2d)    equivalence (cuvscvgam1, cuvscvgam1_2d)
91    
92    real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)    real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
93    real cuvscvgam2(ip1jm), cvuscugam1(ip1jmp1)    real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
94    equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)    equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
95    
96    real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)    real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
97    real cvuscugam2(ip1jmp1), cvscuvgam(ip1jm)    real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
98    equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)    equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
99    
100    real cuscvugam(ip1jmp1)    real cuscvugam((iim + 1) * (jjm + 1))
101    real cuscvugam_2d(iim + 1, jjm + 1)    real cuscvugam_2d(iim + 1, jjm + 1)
102    equivalence (cuscvugam, cuscvugam_2d)    equivalence (cuscvugam, cuscvugam_2d)
103    
104    real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                    real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                
105    
106    real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)    real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
107    real unsair_gam1(ip1jmp1), unsair_gam2(ip1jmp1)    real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
108    equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)    equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
109    
110    real unsairz_gam_2d(iim + 1, jjm)    real unsairz_gam_2d(iim + 1, jjm)
111    real unsairz_gam(ip1jm)    real unsairz_gam((iim + 1) * jjm)
112    equivalence (unsairz_gam, unsairz_gam_2d)    equivalence (unsairz_gam, unsairz_gam_2d)
113    
114    real xprimu(iim + 1), xprimv(iim + 1)    real xprimu(iim + 1), xprimv(iim + 1)
115    
116    save    save
117    
118    contains
119    
120      SUBROUTINE inigeom
121    
122        ! Auteur : P. Le Van
123    
124        ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
125        ! endroits que les aires aireij1_2d, ..., aireij4_2d.
126    
127        ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à
128        ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,
129        ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d
130        ! permettent de passer des vitesses naturelles aux vitesses
131        ! covariantes et contravariantes, ou vice-versa.
132    
133        ! On a :
134        ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
135        ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
136    
137        ! On en tire :
138        ! u(covariant) = cu_2d * cu_2d * u(contravariant)
139        ! v(covariant) = cv_2d * cv_2d * v(contravariant)
140    
141        ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
142        ! et - jm / 2 <= Y <= jm / 2
143    
144        ! x est la longitude du point en radians.
145        ! y est la latitude du point en radians.
146        !
147        ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
148        ! cv(j) = rad * dy / dY
149        ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
150        !
151        ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
152        ! dépendant de j uniquement, sera ici indicé aussi en i pour un
153        ! adressage plus facile en ij.
154    
155        ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
156        ! points u et v. yprimu et yprimv sont respectivement les valeurs
157        ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
158        ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
159        ! respectivement les valeurs de cv_2d aux points u et v.
160    
161        ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
162        ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
163    
164        USE comconst, ONLY : g, omeg, rad
165        USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
166        use conf_gcm_m, ONLY : fxyhypb, ysinus
167        use fxy_m, only: fxy
168        use fxyhyper_m, only: fxyhyper
169        use fxysinus_m, only: fxysinus
170        use jumble, only: new_unit
171        use nr_util, only: pi
172        USE paramet_m, ONLY : iip1, jjp1
173        USE serre, ONLY : alphax, alphay, clat, clon, pxo, pyo, transx, transy
174        ! Modifiés pxo, pyo, transx, transy
175    
176        ! Local:
177        INTEGER i, j, itmax, itmay, iter, unit
178        REAL cvu(iip1, jjp1), cuv(iip1, jjm)
179        REAL ai14, ai23, airez, un4rad2
180        REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
181        REAL coslatm, coslatp, radclatm, radclatp
182        REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
183        REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
184        REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
185        real yprimv(jjm), yprimu(jjp1)
186        REAL gamdi_gdiv, gamdi_grot, gamdi_h
187        REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
188        real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
189             aireij4_2d ! in m2
190        real airuscv2_2d(iim + 1, jjm)
191        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
192        real aivscu2gam_2d(iim + 1, jjm)
193    
194        !------------------------------------------------------------------
195    
196        PRINT *, 'Call sequence information: inigeom'
197    
198        IF (nitergdiv/=2) THEN
199           gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
200        ELSE
201           gamdi_gdiv = 0.
202        END IF
203        IF (nitergrot/=2) THEN
204           gamdi_grot = coefdis / (real(nitergrot)-2.)
205        ELSE
206           gamdi_grot = 0.
207        END IF
208        IF (niterh/=2) THEN
209           gamdi_h = coefdis / (real(niterh)-2.)
210        ELSE
211           gamdi_h = 0.
212        END IF
213    
214        print *, 'gamdi_gdiv = ', gamdi_gdiv
215        print *, "gamdi_grot = ", gamdi_grot
216        print *, "gamdi_h = ", gamdi_h
217    
218        IF (fxyhypb) THEN
219           print *, 'inigeom: Y = latitude, dérivée tangente hyperbolique'
220           CALL fxyhyper(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
221                yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
222                rlonp025, xprimp025)
223        ELSE
224           IF (ysinus) THEN
225              print *, 'inigeom: Y = sin(latitude)'
226              ! Utilisation de f(x, y) avec y = sinus de la latitude
227              CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
228                   rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
229                   xprimm025, rlonp025, xprimp025)
230           ELSE
231              print *, 'Inigeom, Y = Latitude, der. sinusoid .'
232              ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
233    
234              pxo = clon * pi / 180.
235              pyo = 2. * clat * pi / 180.
236    
237              ! determination de transx (pour le zoom) par Newton-Raphson
238    
239              itmax = 10
240              eps = .1E-7
241    
242              xo1 = 0.
243              DO iter = 1, itmax
244                 x1 = xo1
245                 f = x1 + alphax * sin(x1-pxo)
246                 df = 1. + alphax * cos(x1-pxo)
247                 x1 = x1 - f / df
248                 xdm = abs(x1-xo1)
249                 IF (xdm<=eps) EXIT
250                 xo1 = x1
251              END DO
252    
253              transx = xo1
254    
255              itmay = 10
256              eps = .1E-7
257    
258              yo1 = 0.
259              DO iter = 1, itmay
260                 y1 = yo1
261                 f = y1 + alphay * sin(y1-pyo)
262                 df = 1. + alphay * cos(y1-pyo)
263                 y1 = y1 - f / df
264                 ydm = abs(y1-yo1)
265                 IF (ydm<=eps) EXIT
266                 yo1 = y1
267              END DO
268    
269              transy = yo1
270    
271              CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
272                   yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
273                   rlonp025, xprimp025)
274           END IF
275        END IF
276    
277        rlatu(1) = pi / 2.
278        rlatu(jjp1) = -rlatu(1)
279    
280        ! Calcul aux pôles
281    
282        yprimu(1) = 0.
283        yprimu(jjp1) = 0.
284    
285        un4rad2 = 0.25 * rad * rad
286    
287        ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
288        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
289        ! chaque aire_2d(i, j), ainsi que les quatre élongations
290        ! élémentaires cuij et les quatre élongations cvij qui sont
291        ! calculées aux mêmes endroits que les aireij.
292    
293        coslatm = cos(rlatu1(1))
294        radclatm = 0.5 * rad * coslatm
295    
296        aireij1_2d(:iim, 1) = 0.
297        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
298        aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
299        aireij4_2d(:iim, 1) = 0.
300    
301        cuij1(:iim, 1) = 0.
302        cuij2(:iim, 1) = radclatm * xprimp025(:iim)
303        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
304        cuij4(:iim, 1) = 0.
305    
306        cvij1(:iim, 1) = 0.
307        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
308        cvij3(:iim, 1) = cvij2(:iim, 1)
309        cvij4(:iim, 1) = 0.
310    
311        do j = 2, jjm
312           coslatm = cos(rlatu1(j))
313           coslatp = cos(rlatu2(j-1))
314           radclatp = 0.5 * rad * coslatp
315           radclatm = 0.5 * rad * coslatm
316           ai14 = un4rad2 * coslatp * yprimu2(j-1)
317           ai23 = un4rad2 * coslatm * yprimu1(j)
318    
319           aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
320           aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
321           aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
322           aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
323           cuij1(:iim, j) = radclatp * xprimp025(:iim)
324           cuij2(:iim, j) = radclatm * xprimp025(:iim)
325           cuij3(:iim, j) = radclatm * xprimm025(:iim)
326           cuij4(:iim, j) = radclatp * xprimm025(:iim)
327           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
328           cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
329           cvij3(:iim, j) = cvij2(:iim, j)
330           cvij4(:iim, j) = cvij1(:iim, j)
331        end do
332    
333        coslatp = cos(rlatu2(jjm))
334        radclatp = 0.5 * rad * coslatp
335    
336        aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
337        aireij2_2d(:iim, jjp1) = 0.
338        aireij3_2d(:iim, jjp1) = 0.
339        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
340    
341        cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
342        cuij2(:iim, jjp1) = 0.
343        cuij3(:iim, jjp1) = 0.
344        cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
345    
346        cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
347        cvij2(:iim, jjp1) = 0.
348        cvij3(:iim, jjp1) = 0.
349        cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
350    
351        ! Périodicité :
352    
353        cvij1(iip1, :) = cvij1(1, :)
354        cvij2(iip1, :) = cvij2(1, :)
355        cvij3(iip1, :) = cvij3(1, :)
356        cvij4(iip1, :) = cvij4(1, :)
357    
358        cuij1(iip1, :) = cuij1(1, :)
359        cuij2(iip1, :) = cuij2(1, :)
360        cuij3(iip1, :) = cuij3(1, :)
361        cuij4(iip1, :) = cuij4(1, :)
362    
363        aireij1_2d(iip1, :) = aireij1_2d(1, :)
364        aireij2_2d(iip1, :) = aireij2_2d(1, :)
365        aireij3_2d(iip1, :) = aireij3_2d(1, :)
366        aireij4_2d(iip1, :) = aireij4_2d(1, :)
367    
368        DO j = 1, jjp1
369           DO i = 1, iim
370              aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
371                   + aireij3_2d(i, j) + aireij4_2d(i, j)
372              alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
373              alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
374              alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
375              alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
376              alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
377              alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
378              alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
379              alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
380           END DO
381    
382           aire_2d(iip1, j) = aire_2d(1, j)
383           alpha1_2d(iip1, j) = alpha1_2d(1, j)
384           alpha2_2d(iip1, j) = alpha2_2d(1, j)
385           alpha3_2d(iip1, j) = alpha3_2d(1, j)
386           alpha4_2d(iip1, j) = alpha4_2d(1, j)
387           alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
388           alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
389           alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
390           alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
391        END DO
392    
393        DO j = 1, jjp1
394           DO i = 1, iim
395              aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
396                   aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
397              unsaire_2d(i, j) = 1. / aire_2d(i, j)
398              unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
399              unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
400              airesurg_2d(i, j) = aire_2d(i, j) / g
401           END DO
402           aireu_2d(iip1, j) = aireu_2d(1, j)
403           unsaire_2d(iip1, j) = unsaire_2d(1, j)
404           unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
405           unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
406           airesurg_2d(iip1, j) = airesurg_2d(1, j)
407        END DO
408    
409        DO j = 1, jjm
410           DO i = 1, iim
411              airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
412                   aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
413           END DO
414           DO i = 1, iim
415              airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
416                   + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
417              unsairez_2d(i, j) = 1. / airez
418              unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
419              fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
420           END DO
421           airev_2d(iip1, j) = airev_2d(1, j)
422           unsairez_2d(iip1, j) = unsairez_2d(1, j)
423           fext_2d(iip1, j) = fext_2d(1, j)
424           unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
425        END DO
426    
427        ! Calcul des élongations cu_2d, cv_2d, cvu
428    
429        DO j = 1, jjm
430           DO i = 1, iim
431              cv_2d(i, j) = 0.5 * &
432                   (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
433              cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
434                   + cvij3(i, j))
435              cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
436                   + cuij4(i, j + 1))
437              unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
438           END DO
439           DO i = 1, iim
440              cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
441              cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
442              cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
443              cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
444              cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
445           END DO
446           cv_2d(iip1, j) = cv_2d(1, j)
447           cvu(iip1, j) = cvu(1, j)
448           unscv2_2d(iip1, j) = unscv2_2d(1, j)
449           cuv(iip1, j) = cuv(1, j)
450           cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
451           cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
452           cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
453           cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
454           cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
455        END DO
456    
457        DO j = 2, jjm
458           DO i = 1, iim
459              cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
460                   + cuij3(i + 1, j))
461              unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
462              cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
463              cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
464              cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
465              cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
466              cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
467           END DO
468           cu_2d(iip1, j) = cu_2d(1, j)
469           unscu2_2d(iip1, j) = unscu2_2d(1, j)
470           cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
471           cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
472           cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
473           cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
474           cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
475        END DO
476    
477        ! Calcul aux pôles
478    
479        cu_2d(:, 1) = 0.
480        unscu2_2d(:, 1) = 0.
481        cvu(:, 1) = 0.
482    
483        cu_2d(:, jjp1) = 0.
484        unscu2_2d(:, jjp1) = 0.
485        cvu(:, jjp1) = 0.
486    
487        DO j = 1, jjm
488           DO i = 1, iim
489              airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
490              aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
491           END DO
492           airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
493           aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
494        END DO
495    
496        DO j = 2, jjm
497           DO i = 1, iim
498              airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
499              aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
500           END DO
501           airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
502           aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
503        END DO
504    
505        ! Calcul des aires aux pôles :
506    
507        apoln = sum(aire_2d(:iim, 1))
508        apols = sum(aire_2d(:iim, jjp1))
509        unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
510        unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
511        unsapolnga2 = 1. / (apoln**(-gamdi_h))
512        unsapolsga2 = 1. / (apols**(-gamdi_h))
513    
514        ! Changement F. Hourdin calcul conservatif pour fext_2d
515        ! constang_2d contient le produit a * cos (latitude) * omega
516    
517        DO i = 1, iim
518           constang_2d(i, 1) = 0.
519        END DO
520        DO j = 1, jjm - 1
521           DO i = 1, iim
522              constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
523                   * cos(rlatu(j + 1))
524           END DO
525        END DO
526        DO i = 1, iim
527           constang_2d(i, jjp1) = 0.
528        END DO
529    
530        ! Périodicité en longitude
531        DO j = 1, jjp1
532           constang_2d(iip1, j) = constang_2d(1, j)
533        END DO
534    
535        call new_unit(unit)
536        open(unit, file="longitude_latitude.txt", status="replace", action="write")
537        write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
538        write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
539        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
540        write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
541        close(unit)
542    
543      END SUBROUTINE inigeom
544    
545  end module comgeom  end module comgeom
Legend:
Removed from v.60  
changed lines
  Added in v.112
  ViewVC Help
Powered by ViewVC 1.1.21