/[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/Sources/dyn3d/comgeom.f revision 134 by guez, Wed Apr 29 15:47:56 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        ! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des
128        ! coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les
129        ! coefficients cu_2d et cv_2d permettent de passer des vitesses
130        ! naturelles aux vitesses covariantes et contravariantes, ou
131        ! 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 / (nitergdiv - 2)
195        ELSE
196           gamdi_gdiv = 0.
197        END IF
198    
199        IF (nitergrot /= 2) THEN
200           gamdi_grot = coefdis / (nitergrot - 2)
201        ELSE
202           gamdi_grot = 0.
203        END IF
204    
205        IF (niterh /= 2) THEN
206           gamdi_h = coefdis / (niterh - 2)
207        ELSE
208           gamdi_h = 0.
209        END IF
210    
211        print *, 'gamdi_gdiv = ', gamdi_gdiv
212        print *, "gamdi_grot = ", gamdi_grot
213        print *, "gamdi_h = ", gamdi_h
214    
215        CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
216        CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
217    
218        rlatu(1) = pi / 2.
219        rlatu(jjp1) = -rlatu(1)
220    
221        ! Calcul aux pôles
222    
223        yprimu(1) = 0.
224        yprimu(jjp1) = 0.
225    
226        un4rad2 = 0.25 * rad * rad
227    
228        ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
229        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
230        ! chaque aire_2d(i, j), ainsi que les quatre élongations
231        ! élémentaires cuij et les quatre élongations cvij qui sont
232        ! calculées aux mêmes endroits que les aireij.
233    
234        coslatm = cos(rlatu1(1))
235        radclatm = 0.5 * rad * coslatm
236    
237        aireij1_2d(:iim, 1) = 0.
238        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
239        aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
240        aireij4_2d(:iim, 1) = 0.
241    
242        cuij1(:iim, 1) = 0.
243        cuij2(:iim, 1) = radclatm * xprimp025(:iim)
244        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
245        cuij4(:iim, 1) = 0.
246    
247        cvij1(:iim, 1) = 0.
248        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
249        cvij3(:iim, 1) = cvij2(:iim, 1)
250        cvij4(:iim, 1) = 0.
251    
252        do j = 2, jjm
253           coslatm = cos(rlatu1(j))
254           coslatp = cos(rlatu2(j-1))
255           radclatp = 0.5 * rad * coslatp
256           radclatm = 0.5 * rad * coslatm
257           ai14 = un4rad2 * coslatp * yprimu2(j-1)
258           ai23 = un4rad2 * coslatm * yprimu1(j)
259    
260           aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
261           aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
262           aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
263           aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
264           cuij1(:iim, j) = radclatp * xprimp025(:iim)
265           cuij2(:iim, j) = radclatm * xprimp025(:iim)
266           cuij3(:iim, j) = radclatm * xprimm025(:iim)
267           cuij4(:iim, j) = radclatp * xprimm025(:iim)
268           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
269           cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
270           cvij3(:iim, j) = cvij2(:iim, j)
271           cvij4(:iim, j) = cvij1(:iim, j)
272        end do
273    
274        coslatp = cos(rlatu2(jjm))
275        radclatp = 0.5 * rad * coslatp
276    
277        aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
278        aireij2_2d(:iim, jjp1) = 0.
279        aireij3_2d(:iim, jjp1) = 0.
280        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
281    
282        cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
283        cuij2(:iim, jjp1) = 0.
284        cuij3(:iim, jjp1) = 0.
285        cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
286    
287        cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
288        cvij2(:iim, jjp1) = 0.
289        cvij3(:iim, jjp1) = 0.
290        cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
291    
292        ! Périodicité :
293    
294        cvij1(iip1, :) = cvij1(1, :)
295        cvij2(iip1, :) = cvij2(1, :)
296        cvij3(iip1, :) = cvij3(1, :)
297        cvij4(iip1, :) = cvij4(1, :)
298    
299        cuij1(iip1, :) = cuij1(1, :)
300        cuij2(iip1, :) = cuij2(1, :)
301        cuij3(iip1, :) = cuij3(1, :)
302        cuij4(iip1, :) = cuij4(1, :)
303    
304        aireij1_2d(iip1, :) = aireij1_2d(1, :)
305        aireij2_2d(iip1, :) = aireij2_2d(1, :)
306        aireij3_2d(iip1, :) = aireij3_2d(1, :)
307        aireij4_2d(iip1, :) = aireij4_2d(1, :)
308    
309        DO j = 1, jjp1
310           DO i = 1, iim
311              aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
312                   + aireij3_2d(i, j) + aireij4_2d(i, j)
313              alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
314              alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
315              alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
316              alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
317              alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
318              alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
319              alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
320              alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
321           END DO
322    
323           aire_2d(iip1, j) = aire_2d(1, j)
324           alpha1_2d(iip1, j) = alpha1_2d(1, j)
325           alpha2_2d(iip1, j) = alpha2_2d(1, j)
326           alpha3_2d(iip1, j) = alpha3_2d(1, j)
327           alpha4_2d(iip1, j) = alpha4_2d(1, j)
328           alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
329           alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
330           alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
331           alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
332        END DO
333    
334        DO j = 1, jjp1
335           DO i = 1, iim
336              aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
337                   aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
338              unsaire_2d(i, j) = 1. / aire_2d(i, j)
339              unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
340              unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
341              airesurg_2d(i, j) = aire_2d(i, j) / g
342           END DO
343           aireu_2d(iip1, j) = aireu_2d(1, j)
344           unsaire_2d(iip1, j) = unsaire_2d(1, j)
345           unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
346           unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
347           airesurg_2d(iip1, j) = airesurg_2d(1, j)
348        END DO
349    
350        DO j = 1, jjm
351           DO i = 1, iim
352              airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
353                   aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
354           END DO
355           DO i = 1, iim
356              airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
357                   + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
358              unsairez_2d(i, j) = 1. / airez
359              unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
360              fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
361           END DO
362           airev_2d(iip1, j) = airev_2d(1, j)
363           unsairez_2d(iip1, j) = unsairez_2d(1, j)
364           fext_2d(iip1, j) = fext_2d(1, j)
365           unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
366        END DO
367    
368        ! Calcul des élongations cu_2d, cv_2d, cvu
369    
370        DO j = 1, jjm
371           DO i = 1, iim
372              cv_2d(i, j) = 0.5 * &
373                   (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
374              cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
375                   + cvij3(i, j))
376              cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
377                   + cuij4(i, j + 1))
378              unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
379           END DO
380           DO i = 1, iim
381              cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
382              cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
383              cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
384              cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
385              cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
386           END DO
387           cv_2d(iip1, j) = cv_2d(1, j)
388           cvu(iip1, j) = cvu(1, j)
389           unscv2_2d(iip1, j) = unscv2_2d(1, j)
390           cuv(iip1, j) = cuv(1, j)
391           cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
392           cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
393           cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
394           cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
395           cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
396        END DO
397    
398        DO j = 2, jjm
399           DO i = 1, iim
400              cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
401                   + cuij3(i + 1, j))
402              unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
403              cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
404              cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
405              cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
406              cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
407              cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
408           END DO
409           cu_2d(iip1, j) = cu_2d(1, j)
410           unscu2_2d(iip1, j) = unscu2_2d(1, j)
411           cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
412           cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
413           cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
414           cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
415           cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
416        END DO
417    
418        ! Calcul aux pôles
419    
420        cu_2d(:, 1) = 0.
421        unscu2_2d(:, 1) = 0.
422        cvu(:, 1) = 0.
423    
424        cu_2d(:, jjp1) = 0.
425        unscu2_2d(:, jjp1) = 0.
426        cvu(:, jjp1) = 0.
427    
428        DO j = 1, jjm
429           DO i = 1, iim
430              airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
431              aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
432           END DO
433           airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
434           aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
435        END DO
436    
437        DO j = 2, jjm
438           DO i = 1, iim
439              airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
440              aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
441           END DO
442           airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
443           aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
444        END DO
445    
446        ! Calcul des aires aux pôles :
447    
448        apoln = sum(aire_2d(:iim, 1))
449        apols = sum(aire_2d(:iim, jjp1))
450        unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
451        unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
452        unsapolnga2 = 1. / (apoln**(-gamdi_h))
453        unsapolsga2 = 1. / (apols**(-gamdi_h))
454    
455        ! Changement F. Hourdin calcul conservatif pour fext_2d
456        ! constang_2d contient le produit a * cos (latitude) * omega
457    
458        DO i = 1, iim
459           constang_2d(i, 1) = 0.
460        END DO
461        DO j = 1, jjm - 1
462           DO i = 1, iim
463              constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
464                   * cos(rlatu(j + 1))
465           END DO
466        END DO
467        DO i = 1, iim
468           constang_2d(i, jjp1) = 0.
469        END DO
470    
471        ! Périodicité en longitude
472        DO j = 1, jjp1
473           constang_2d(iip1, j) = constang_2d(1, j)
474        END DO
475    
476        call new_unit(unit)
477        open(unit, file="longitude_latitude.txt", status="replace", action="write")
478        write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
479        write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
480        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
481        write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
482        close(unit)
483    
484      END SUBROUTINE inigeom
485    
486  end module comgeom  end module comgeom
Legend:
Removed from v.3  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21