/[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 140 by guez, Fri Jun 5 18:58:06 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 cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
68    ! (latitudes of points of the "scalar" and "u" grid, in rad)    real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
   
   real rlatv(jjm)  
   ! (latitudes of points of the "v" grid, in rad, in decreasing order)  
   
   real rlonu(iim + 1) ! longitudes of points of the "u" grid, in rad  
   
   real rlonv(iim + 1)  
   ! (longitudes of points of the "scalar" and "v" grid, in rad)  
   
   real cuvsurcv_2d(iim + 1,jjm),cvsurcuv_2d(iim + 1,jjm)    
   real cuvsurcv(ip1jm),cvsurcuv(ip1jm)  
69    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
70    
71    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)
72    real cvusurcu(ip1jmp1),cusurcvu(ip1jmp1)    ! no dimension
73      real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
74      ! no dimension
75    equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)    equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
76    
77    real cuvscvgam1_2d(iim + 1,jjm)    real cuvscvgam1_2d(iim + 1, jjm)
78    real cuvscvgam1(ip1jm)    real cuvscvgam1((iim + 1) * jjm)
79    equivalence (cuvscvgam1, cuvscvgam1_2d)    equivalence (cuvscvgam1, cuvscvgam1_2d)
80    
81    real cuvscvgam2_2d(iim + 1,jjm),cvuscugam1_2d(iim + 1,jjm + 1)    real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
82    real cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1)    real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
83    equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)    equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
84    
85    real cvuscugam2_2d(iim + 1,jjm + 1),cvscuvgam_2d(iim + 1,jjm)    real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
86    real cvuscugam2(ip1jmp1),cvscuvgam(ip1jm)    real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
87    equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)    equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
88    
89    real cuscvugam(ip1jmp1)    real cuscvugam((iim + 1) * (jjm + 1))
90    real cuscvugam_2d(iim + 1,jjm + 1)    real cuscvugam_2d(iim + 1, jjm + 1)
91    equivalence (cuscvugam, cuscvugam_2d)    equivalence (cuscvugam, cuscvugam_2d)
92    
93    real unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2                    real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                
94    
95    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)
96    real unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1)    real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
97    equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)    equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
98    
99    real unsairz_gam_2d(iim + 1,jjm)    real unsairz_gam_2d(iim + 1, jjm)
100    real unsairz_gam(ip1jm)    real unsairz_gam((iim + 1) * jjm)
101    equivalence (unsairz_gam, unsairz_gam_2d)    equivalence (unsairz_gam, unsairz_gam_2d)
102    
103    real aivscu2gam_2d(iim + 1,jjm),aiuscv2gam_2d(iim + 1,jjm)    save
104    real aivscu2gam(ip1jm),aiuscv2gam(ip1jm)  
105    equivalence (aivscu2gam, aivscu2gam_2d), (aiuscv2gam, aiuscv2gam_2d)  contains
106    
107    real xprimu(iim + 1),xprimv(iim + 1)    SUBROUTINE inigeom
108    
109    save      ! Auteur : P. Le Van
110    
111        ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
112        ! endroits que les aires aireij1_2d, ..., aireij4_2d.
113    
114        ! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des
115        ! coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les
116        ! coefficients cu_2d et cv_2d permettent de passer des vitesses
117        ! naturelles aux vitesses covariantes et contravariantes, ou
118        ! vice-versa.
119    
120        ! On a :
121        ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
122        ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
123    
124        ! On en tire :
125        ! u(covariant) = cu_2d * cu_2d * u(contravariant)
126        ! v(covariant) = cv_2d * cv_2d * v(contravariant)
127    
128        ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
129        ! et - jm / 2 <= Y <= jm / 2
130    
131        ! x est la longitude du point en radians.
132        ! y est la latitude du point en radians.
133        !
134        ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
135        ! cv(j) = rad * dy / dY
136        ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
137        !
138        ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
139        ! dépendant de j uniquement, sera ici indicé aussi en i pour un
140        ! adressage plus facile en ij.
141    
142        ! cv_2d est aux points v. cu_2d est aux points
143        ! u. Cf. "inigeom.txt".
144    
145        USE comconst, ONLY : g, omeg, rad
146        USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
147        use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
148             yprimu1, yprimu2, rlonu, rlonv
149        use jumble, only: new_unit
150        use nr_util, only: pi
151        USE paramet_m, ONLY : iip1, jjp1
152    
153        ! Local:
154        INTEGER i, j, unit
155        REAL ai14, ai23, airez, un4rad2
156        REAL coslatm, coslatp, radclatm, radclatp
157        REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
158        REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
159        REAL gamdi_gdiv, gamdi_grot, gamdi_h
160        real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
161             aireij4_2d ! in m2
162    
163        !------------------------------------------------------------------
164    
165        PRINT *, 'Call sequence information: inigeom'
166    
167        IF (nitergdiv /= 2) THEN
168           gamdi_gdiv = coefdis / (nitergdiv - 2)
169        ELSE
170           gamdi_gdiv = 0.
171        END IF
172    
173        IF (nitergrot /= 2) THEN
174           gamdi_grot = coefdis / (nitergrot - 2)
175        ELSE
176           gamdi_grot = 0.
177        END IF
178    
179        IF (niterh /= 2) THEN
180           gamdi_h = coefdis / (niterh - 2)
181        ELSE
182           gamdi_h = 0.
183        END IF
184    
185        print *, 'gamdi_gdiv = ', gamdi_gdiv
186        print *, "gamdi_grot = ", gamdi_grot
187        print *, "gamdi_h = ", gamdi_h
188    
189        un4rad2 = 0.25 * rad * rad
190    
191        ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
192        ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
193        ! chaque aire_2d(i, j), ainsi que les quatre élongations
194        ! élémentaires cuij et les quatre élongations cvij qui sont
195        ! calculées aux mêmes endroits que les aireij.
196    
197        coslatm = cos(rlatu1(1))
198        radclatm = 0.5 * rad * coslatm
199    
200        aireij1_2d(:iim, 1) = 0.
201        aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
202        aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
203        aireij4_2d(:iim, 1) = 0.
204    
205        cuij1(:iim, 1) = 0.
206        cuij2(:iim, 1) = radclatm * xprimp025(:iim)
207        cuij3(:iim, 1) = radclatm * xprimm025(:iim)
208        cuij4(:iim, 1) = 0.
209    
210        cvij1(:iim, 1) = 0.
211        cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
212        cvij3(:iim, 1) = cvij2(:iim, 1)
213        cvij4(:iim, 1) = 0.
214    
215        do j = 2, jjm
216           coslatm = cos(rlatu1(j))
217           coslatp = cos(rlatu2(j-1))
218           radclatp = 0.5 * rad * coslatp
219           radclatm = 0.5 * rad * coslatm
220           ai14 = un4rad2 * coslatp * yprimu2(j-1)
221           ai23 = un4rad2 * coslatm * yprimu1(j)
222    
223           aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
224           aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
225           aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
226           aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
227           cuij1(:iim, j) = radclatp * xprimp025(:iim)
228           cuij2(:iim, j) = radclatm * xprimp025(:iim)
229           cuij3(:iim, j) = radclatm * xprimm025(:iim)
230           cuij4(:iim, j) = radclatp * xprimm025(:iim)
231           cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
232           cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
233           cvij3(:iim, j) = cvij2(:iim, j)
234           cvij4(:iim, j) = cvij1(:iim, j)
235        end do
236    
237        coslatp = cos(rlatu2(jjm))
238        radclatp = 0.5 * rad * coslatp
239    
240        aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
241        aireij2_2d(:iim, jjp1) = 0.
242        aireij3_2d(:iim, jjp1) = 0.
243        aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
244    
245        cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
246        cuij2(:iim, jjp1) = 0.
247        cuij3(:iim, jjp1) = 0.
248        cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
249    
250        cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
251        cvij2(:iim, jjp1) = 0.
252        cvij3(:iim, jjp1) = 0.
253        cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
254    
255        ! Périodicité :
256    
257        cvij1(iip1, :) = cvij1(1, :)
258        cvij2(iip1, :) = cvij2(1, :)
259        cvij3(iip1, :) = cvij3(1, :)
260        cvij4(iip1, :) = cvij4(1, :)
261    
262        cuij1(iip1, :) = cuij1(1, :)
263        cuij2(iip1, :) = cuij2(1, :)
264        cuij3(iip1, :) = cuij3(1, :)
265        cuij4(iip1, :) = cuij4(1, :)
266    
267        aireij1_2d(iip1, :) = aireij1_2d(1, :)
268        aireij2_2d(iip1, :) = aireij2_2d(1, :)
269        aireij3_2d(iip1, :) = aireij3_2d(1, :)
270        aireij4_2d(iip1, :) = aireij4_2d(1, :)
271    
272        DO j = 1, jjp1
273           DO i = 1, iim
274              aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
275                   + aireij3_2d(i, j) + aireij4_2d(i, j)
276              alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
277              alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
278              alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
279              alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
280              alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
281              alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
282              alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
283              alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
284           END DO
285    
286           aire_2d(iip1, j) = aire_2d(1, j)
287           alpha1_2d(iip1, j) = alpha1_2d(1, j)
288           alpha2_2d(iip1, j) = alpha2_2d(1, j)
289           alpha3_2d(iip1, j) = alpha3_2d(1, j)
290           alpha4_2d(iip1, j) = alpha4_2d(1, j)
291           alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
292           alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
293           alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
294           alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
295        END DO
296    
297        DO j = 1, jjp1
298           DO i = 1, iim
299              aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
300                   aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
301              unsaire_2d(i, j) = 1. / aire_2d(i, j)
302              unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
303              unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
304              airesurg_2d(i, j) = aire_2d(i, j) / g
305           END DO
306           aireu_2d(iip1, j) = aireu_2d(1, j)
307           unsaire_2d(iip1, j) = unsaire_2d(1, j)
308           unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
309           unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
310           airesurg_2d(iip1, j) = airesurg_2d(1, j)
311        END DO
312    
313        DO j = 1, jjm
314           DO i = 1, iim
315              airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
316                   aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
317           END DO
318           DO i = 1, iim
319              airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
320                   + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
321              unsairez_2d(i, j) = 1. / airez
322              unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
323              fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
324           END DO
325           airev_2d(iip1, j) = airev_2d(1, j)
326           unsairez_2d(iip1, j) = unsairez_2d(1, j)
327           fext_2d(iip1, j) = fext_2d(1, j)
328           unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
329        END DO
330    
331        ! Calcul des élongations cu_2d, cv_2d
332    
333        DO j = 1, jjm
334           DO i = 1, iim
335              cv_2d(i, j) = 0.5 * &
336                   (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
337              unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
338           END DO
339           DO i = 1, iim
340              cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
341              cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
342              cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
343              cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
344              cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
345           END DO
346           cv_2d(iip1, j) = cv_2d(1, j)
347           unscv2_2d(iip1, j) = unscv2_2d(1, j)
348           cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
349           cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
350           cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
351           cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
352           cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
353        END DO
354    
355        DO j = 2, jjm
356           DO i = 1, iim
357              cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
358                   + cuij3(i + 1, j))
359              unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
360              cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
361              cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
362              cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
363              cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
364              cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
365           END DO
366           cu_2d(iip1, j) = cu_2d(1, j)
367           unscu2_2d(iip1, j) = unscu2_2d(1, j)
368           cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
369           cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
370           cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
371           cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
372           cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
373        END DO
374    
375        ! Calcul aux pôles
376    
377        cu_2d(:, 1) = 0.
378        unscu2_2d(:, 1) = 0.
379    
380        cu_2d(:, jjp1) = 0.
381        unscu2_2d(:, jjp1) = 0.
382    
383        ! Calcul des aires aux pôles :
384    
385        apoln = sum(aire_2d(:iim, 1))
386        apols = sum(aire_2d(:iim, jjp1))
387        unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
388        unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
389        unsapolnga2 = 1. / (apoln**(-gamdi_h))
390        unsapolsga2 = 1. / (apols**(-gamdi_h))
391    
392        ! Changement F. Hourdin calcul conservatif pour fext_2d
393        ! constang_2d contient le produit a * cos (latitude) * omega
394    
395        DO i = 1, iim
396           constang_2d(i, 1) = 0.
397        END DO
398        DO j = 1, jjm - 1
399           DO i = 1, iim
400              constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
401                   * cos(rlatu(j + 1))
402           END DO
403        END DO
404        DO i = 1, iim
405           constang_2d(i, jjp1) = 0.
406        END DO
407    
408        ! Périodicité en longitude
409        DO j = 1, jjp1
410           constang_2d(iip1, j) = constang_2d(1, j)
411        END DO
412    
413        call new_unit(unit)
414        open(unit, file="longitude_latitude.txt", status="replace", action="write")
415        write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
416        write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
417        write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
418        write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
419        close(unit)
420    
421      END SUBROUTINE inigeom
422    
423  end module comgeom  end module comgeom
Legend:
Removed from v.3  
changed lines
  Added in v.140
  ViewVC Help
Powered by ViewVC 1.1.21