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

Diff of /trunk/Sources/dyn3d/comgeom.f

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

trunk/libf/dyn3d/comgeom.f90 revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/dyn3d/comgeom.f revision 120 by guez, Tue Jan 13 14:56:15 2015 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)    real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
10    real cu(ip1jmp1), cv(ip1jm)    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)    real unscu2_2d(iim + 1, jjm + 1) ! in m-2
14    real unscu2(ip1jmp1)    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)    real unscv2_2d(iim + 1, jjm) ! in m-2
18    real unscv2(ip1jm)    real unscv2((iim + 1) * jjm) ! in m-2
19    equivalence (unscv2, unscv2_2d)    equivalence (unscv2, unscv2_2d)
20    
21    real aire_2d(iim + 1,jjm + 1), airesurg_2d(iim + 1,jjm + 1)    real aire((iim + 1) * (jjm + 1)), aire_2d(iim + 1, jjm + 1) ! in m2
22    real aire(ip1jmp1), 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)              real aireu_2d(iim + 1, jjm + 1) ! in m2
26    real aireu(ip1jmp1)    real aireu((iim + 1) * (jjm + 1)) ! in m2
27    equivalence (aireu, aireu_2d)    equivalence (aireu, aireu_2d)
28    
29    real airev_2d(iim + 1,jjm),unsaire_2d(iim + 1,jjm + 1)    real airev((iim + 1) * jjm), airev_2d(iim + 1, jjm) ! in m2
30    real airev(ip1jm),unsaire(ip1jmp1)    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    real apoln, apols ! in m2
34    
35    real unsairez_2d(iim + 1,jjm),airuscv2_2d(iim + 1,jjm)    real unsairez_2d(iim + 1, jjm)
36    real unsairez(ip1jm),airuscv2(ip1jm)    real unsairez((iim + 1) * jjm)
37    equivalence (unsairez, unsairez_2d), (airuscv2, airuscv2_2d)    equivalence (unsairez, unsairez_2d)
38    
39    real airvscu2_2d(iim + 1,jjm)          real alpha1_2d(iim + 1, jjm + 1)
40    real airvscu2(ip1jm)    real alpha1((iim + 1) * (jjm + 1))
41    equivalence (airvscu2, airvscu2_2d)    equivalence (alpha1, alpha1_2d)
   
   real aireij1_2d(iim + 1,jjm + 1),aireij2_2d(iim + 1,jjm + 1)  
   real aireij1(ip1jmp1),aireij2(ip1jmp1)  
   equivalence (aireij1, aireij1_2d), (aireij2, aireij2_2d)  
   
   real aireij3(ip1jmp1)  
   real aireij3_2d(iim + 1,jjm + 1)        
   equivalence (aireij3, aireij3_2d)  
   
   real aireij4_2d(iim + 1,jjm + 1), alpha1_2d(iim + 1,jjm + 1)  
   real aireij4(ip1jmp1), alpha1(ip1jmp1)  
   equivalence (aireij4, aireij4_2d), (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 88  module comgeom Line 75  module comgeom
75    real rlonv(iim + 1)    real rlonv(iim + 1)
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)      real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
79    real cuvsurcv(ip1jm),cvsurcuv(ip1jm)    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    real cvusurcu(ip1jmp1),cusurcvu(ip1jmp1)    ! no dimension
84      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 aivscu2gam_2d(iim + 1,jjm),aiuscv2gam_2d(iim + 1,jjm)    real xprimu(iim + 1), xprimv(iim + 1)
   real aivscu2gam(ip1jm),aiuscv2gam(ip1jm)  
   equivalence (aivscu2gam, aivscu2gam_2d), (aiuscv2gam, aiuscv2gam_2d)  
   
   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 fxhyp_m, only: fxhyp
167        use fyhyp_m, only: fyhyp
168        use jumble, only: new_unit
169        use nr_util, only: pi
170        USE paramet_m, ONLY : iip1, jjp1
171    
172        ! Local:
173        INTEGER i, j, unit
174        REAL cvu(iip1, jjp1), cuv(iip1, jjm)
175        REAL ai14, ai23, airez, un4rad2
176        REAL coslatm, coslatp, radclatm, radclatp
177        REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
178        REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
179        REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
180        real yprimu(jjp1)
181        REAL gamdi_gdiv, gamdi_grot, gamdi_h
182        REAL xprimm025(iip1), xprimp025(iip1)
183        real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
184             aireij4_2d ! in m2
185        real airuscv2_2d(iim + 1, jjm)
186        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
187        real aivscu2gam_2d(iim + 1, jjm)
188    
189        !------------------------------------------------------------------
190    
191        PRINT *, 'Call sequence information: inigeom'
192    
193        IF (nitergdiv/=2) THEN
194           gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
195        ELSE
196           gamdi_gdiv = 0.
197        END IF
198        IF (nitergrot/=2) THEN
199           gamdi_grot = coefdis / (real(nitergrot)-2.)
200        ELSE
201           gamdi_grot = 0.
202        END IF
203        IF (niterh/=2) THEN
204           gamdi_h = coefdis / (real(niterh)-2.)
205        ELSE
206           gamdi_h = 0.
207        END IF
208    
209        print *, 'gamdi_gdiv = ', gamdi_gdiv
210        print *, "gamdi_grot = ", gamdi_grot
211        print *, "gamdi_h = ", gamdi_h
212    
213        CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
214        CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
215    
216        rlatu(1) = pi / 2.
217        rlatu(jjp1) = -rlatu(1)
218    
219        ! Calcul aux pôles
220    
221        yprimu(1) = 0.
222        yprimu(jjp1) = 0.
223    
224        un4rad2 = 0.25 * rad * rad
225    
226        ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
227        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
228        ! chaque aire_2d(i, j), ainsi que les quatre élongations
229        ! élémentaires cuij et les quatre élongations cvij qui sont
230        ! calculées aux mêmes endroits que les aireij.
231    
232        coslatm = cos(rlatu1(1))
233        radclatm = 0.5 * rad * coslatm
234    
235        aireij1_2d(:iim, 1) = 0.
236        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
237        aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
238        aireij4_2d(:iim, 1) = 0.
239    
240        cuij1(:iim, 1) = 0.
241        cuij2(:iim, 1) = radclatm * xprimp025(:iim)
242        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
243        cuij4(:iim, 1) = 0.
244    
245        cvij1(:iim, 1) = 0.
246        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
247        cvij3(:iim, 1) = cvij2(:iim, 1)
248        cvij4(:iim, 1) = 0.
249    
250        do j = 2, jjm
251           coslatm = cos(rlatu1(j))
252           coslatp = cos(rlatu2(j-1))
253           radclatp = 0.5 * rad * coslatp
254           radclatm = 0.5 * rad * coslatm
255           ai14 = un4rad2 * coslatp * yprimu2(j-1)
256           ai23 = un4rad2 * coslatm * yprimu1(j)
257    
258           aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
259           aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
260           aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
261           aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
262           cuij1(:iim, j) = radclatp * xprimp025(:iim)
263           cuij2(:iim, j) = radclatm * xprimp025(:iim)
264           cuij3(:iim, j) = radclatm * xprimm025(:iim)
265           cuij4(:iim, j) = radclatp * xprimm025(:iim)
266           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
267           cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
268           cvij3(:iim, j) = cvij2(:iim, j)
269           cvij4(:iim, j) = cvij1(:iim, j)
270        end do
271    
272        coslatp = cos(rlatu2(jjm))
273        radclatp = 0.5 * rad * coslatp
274    
275        aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
276        aireij2_2d(:iim, jjp1) = 0.
277        aireij3_2d(:iim, jjp1) = 0.
278        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
279    
280        cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
281        cuij2(:iim, jjp1) = 0.
282        cuij3(:iim, jjp1) = 0.
283        cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
284    
285        cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
286        cvij2(:iim, jjp1) = 0.
287        cvij3(:iim, jjp1) = 0.
288        cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
289    
290        ! Périodicité :
291    
292        cvij1(iip1, :) = cvij1(1, :)
293        cvij2(iip1, :) = cvij2(1, :)
294        cvij3(iip1, :) = cvij3(1, :)
295        cvij4(iip1, :) = cvij4(1, :)
296    
297        cuij1(iip1, :) = cuij1(1, :)
298        cuij2(iip1, :) = cuij2(1, :)
299        cuij3(iip1, :) = cuij3(1, :)
300        cuij4(iip1, :) = cuij4(1, :)
301    
302        aireij1_2d(iip1, :) = aireij1_2d(1, :)
303        aireij2_2d(iip1, :) = aireij2_2d(1, :)
304        aireij3_2d(iip1, :) = aireij3_2d(1, :)
305        aireij4_2d(iip1, :) = aireij4_2d(1, :)
306    
307        DO j = 1, jjp1
308           DO i = 1, iim
309              aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
310                   + aireij3_2d(i, j) + aireij4_2d(i, j)
311              alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
312              alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
313              alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
314              alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
315              alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
316              alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
317              alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
318              alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
319           END DO
320    
321           aire_2d(iip1, j) = aire_2d(1, j)
322           alpha1_2d(iip1, j) = alpha1_2d(1, j)
323           alpha2_2d(iip1, j) = alpha2_2d(1, j)
324           alpha3_2d(iip1, j) = alpha3_2d(1, j)
325           alpha4_2d(iip1, j) = alpha4_2d(1, j)
326           alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
327           alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
328           alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
329           alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
330        END DO
331    
332        DO j = 1, jjp1
333           DO i = 1, iim
334              aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
335                   aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
336              unsaire_2d(i, j) = 1. / aire_2d(i, j)
337              unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
338              unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
339              airesurg_2d(i, j) = aire_2d(i, j) / g
340           END DO
341           aireu_2d(iip1, j) = aireu_2d(1, j)
342           unsaire_2d(iip1, j) = unsaire_2d(1, j)
343           unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
344           unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
345           airesurg_2d(iip1, j) = airesurg_2d(1, j)
346        END DO
347    
348        DO j = 1, jjm
349           DO i = 1, iim
350              airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
351                   aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
352           END DO
353           DO i = 1, iim
354              airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
355                   + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
356              unsairez_2d(i, j) = 1. / airez
357              unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
358              fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
359           END DO
360           airev_2d(iip1, j) = airev_2d(1, j)
361           unsairez_2d(iip1, j) = unsairez_2d(1, j)
362           fext_2d(iip1, j) = fext_2d(1, j)
363           unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
364        END DO
365    
366        ! Calcul des élongations cu_2d, cv_2d, cvu
367    
368        DO j = 1, jjm
369           DO i = 1, iim
370              cv_2d(i, j) = 0.5 * &
371                   (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
372              cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
373                   + cvij3(i, j))
374              cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
375                   + cuij4(i, j + 1))
376              unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
377           END DO
378           DO i = 1, iim
379              cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
380              cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
381              cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
382              cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
383              cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
384           END DO
385           cv_2d(iip1, j) = cv_2d(1, j)
386           cvu(iip1, j) = cvu(1, j)
387           unscv2_2d(iip1, j) = unscv2_2d(1, j)
388           cuv(iip1, j) = cuv(1, j)
389           cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
390           cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
391           cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
392           cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
393           cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
394        END DO
395    
396        DO j = 2, jjm
397           DO i = 1, iim
398              cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
399                   + cuij3(i + 1, j))
400              unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
401              cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
402              cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
403              cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
404              cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
405              cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
406           END DO
407           cu_2d(iip1, j) = cu_2d(1, j)
408           unscu2_2d(iip1, j) = unscu2_2d(1, j)
409           cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
410           cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
411           cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
412           cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
413           cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
414        END DO
415    
416        ! Calcul aux pôles
417    
418        cu_2d(:, 1) = 0.
419        unscu2_2d(:, 1) = 0.
420        cvu(:, 1) = 0.
421    
422        cu_2d(:, jjp1) = 0.
423        unscu2_2d(:, jjp1) = 0.
424        cvu(:, jjp1) = 0.
425    
426        DO j = 1, jjm
427           DO i = 1, iim
428              airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
429              aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
430           END DO
431           airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
432           aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
433        END DO
434    
435        DO j = 2, jjm
436           DO i = 1, iim
437              airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
438              aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
439           END DO
440           airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
441           aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
442        END DO
443    
444        ! Calcul des aires aux pôles :
445    
446        apoln = sum(aire_2d(:iim, 1))
447        apols = sum(aire_2d(:iim, jjp1))
448        unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
449        unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
450        unsapolnga2 = 1. / (apoln**(-gamdi_h))
451        unsapolsga2 = 1. / (apols**(-gamdi_h))
452    
453        ! Changement F. Hourdin calcul conservatif pour fext_2d
454        ! constang_2d contient le produit a * cos (latitude) * omega
455    
456        DO i = 1, iim
457           constang_2d(i, 1) = 0.
458        END DO
459        DO j = 1, jjm - 1
460           DO i = 1, iim
461              constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
462                   * cos(rlatu(j + 1))
463           END DO
464        END DO
465        DO i = 1, iim
466           constang_2d(i, jjp1) = 0.
467        END DO
468    
469        ! Périodicité en longitude
470        DO j = 1, jjp1
471           constang_2d(iip1, j) = constang_2d(1, j)
472        END DO
473    
474        call new_unit(unit)
475        open(unit, file="longitude_latitude.txt", status="replace", action="write")
476        write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
477        write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
478        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
479        write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
480        close(unit)
481    
482      END SUBROUTINE inigeom
483    
484  end module comgeom  end module comgeom
Legend:
Removed from v.3  
changed lines
  Added in v.120
  ViewVC Help
Powered by ViewVC 1.1.21