/[lmdze]/trunk/libf/dyn3d/inigeom.f90
ViewVC logotype

Diff of /trunk/libf/dyn3d/inigeom.f90

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

revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC revision 39 by guez, Tue Jan 25 15:11:05 2011 UTC
# Line 16  contains Line 16  contains
16      ! tangente hyperbolique      ! tangente hyperbolique
17      ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2)      ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2)
18    
19      ! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles      ! les coef. (cu_2d, cv_2d) permettent de passer des vitesses naturelles
20      !      aux vitesses covariantes et contravariantes , ou vice-versa ...      ! aux vitesses covariantes et contravariantes, ou vice-versa
21    
22      ! on a :  u (covariant) = cu_2d * u (naturel)   , u(contrav)= u(nat)/cu_2d      ! on a :
23      !         v (covariant) = cv_2d * v (naturel)   , v(contrav)= v(nat)/cv_2d      ! u (covariant) = cu_2d * u (naturel), u(contrav)= u(nat)/cu_2d
24        ! v (covariant) = cv_2d * v (naturel), v(contrav)= v(nat)/cv_2d
25    
26        ! on en tire :
27        ! u(covariant) = cu_2d * cu_2d * u(contravariant)
28        ! v(covariant) = cv_2d * cv_2d * v(contravariant)
29    
30        ! on a l'application (x(X), y(Y)) avec - im/2 +1 <= X <= im/2
31        ! et - jm/2 <= Y <= jm/2
32    
33        ! x est la longitude du point en radians.
34        ! y est la latitude du point en radians.
35        !
36        ! on a : cu_2d(i, j) = rad * cos(y) * dx/dX
37        ! cv(j) = rad * dy/dY
38        ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
39        !
40        ! y, dx/dX, dy/dY calcules aux points concernes
41        ! cv, bien que dependant de j uniquement, sera ici indice aussi en i
42        ! pour un adressage plus facile en ij.
43    
44        ! aux points u et v,
45        ! xprimu et xprimv sont respectivement les valeurs de dx/dX
46        ! yprimu et yprimv sont respectivement les valeurs de dy/dY
47        ! rlatu et rlatv sont respectivement les valeurs de la latitude
48        ! cvu et cv_2d sont respectivement les valeurs de cv_2d
49    
50        ! aux points u, v, scalaires, et z
51        ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
52        ! Cf. "inigeom.txt".
53    
54      !       on en tire :  u(covariant) = cu_2d * cu_2d * u(contravariant)      USE comconst, ONLY : g, omeg, rad
     !                     v(covariant) = cv_2d * cv_2d * v(contravariant)  
   
     !     on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2  
     !                                                          =     =  
     !                                           et   - jm/2    <  Y  < jm/2  
     !                                                          =     =  
   
     !      .  x  est la longitude du point  en radians       .  
     !      .  y  est la  latitude du point  en radians       .  
     !      .                                                 .  
     !      .  on a :  cu_2d(i, j) = rad * COS(y) * dx/dX         .  
     !      .          cv( j ) = rad          * dy/dY         .  
     !      .        aire_2d(i, j) =  cu_2d(i, j) * cv(j)             .  
     !      .                                                 .  
     !      . y, dx/dX, dy/dY calcules aux points concernes   .  
     !                                                           ,  
     !    cv , bien que dependant de j uniquement, sera ici indice aussi en i  
     !          pour un adressage plus facile en  ij  .  
   
     !  **************  aux points  u  et  v ,           *****************  
     !      xprimu et xprimv sont respectivement les valeurs de  dx/dX  
     !      yprimu et yprimv    .  .  .  .  .  .  .  .  .  .  .  dy/dY  
     !      rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude  
     !      cvu    et   cv_2d      .  .  .  .  .  .  .  .  .  .  .    cv_2d  
   
     !  **************  aux points u, v, scalaires, et z  ****************  
     !      cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de    cu_2d  
   
     !         Exemple de distribution de variables sur la grille dans le  
     !             domaine de travail ( X, Y ) .  
     !                  DX=DY= 1  
   
     !        +     represente  un  point scalaire ( p.exp  la pression )  
     !        >     represente  la composante zonale du  vent  
     !        V     represente  la composante meridienne du vent  
     !        o     represente  la  vorticite  
   
     !     ----  , car aux poles , les comp.zonales covariantes sont nulles  
   
     !         i ->  
   
     !         1      2      3      4      5      6      7      8  
     !  j  
     !  v  1   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --  
   
     !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
     !     2   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
     !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
     !     3   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
     !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
     !     4   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >  
   
     !         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o  
   
     !     5   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --  
   
     !      Ci-dessus,  on voit que le nombre de pts.en longitude est egal  
     !                 a   IM = 8  
     !      De meme ,   le nombre d'intervalles entre les 2 poles est egal  
     !                 a   JM = 4  
   
     !      Les points scalaires ( + ) correspondent donc a des valeurs  
     !       entieres  de  i ( 1 a IM )   et  de  j ( 1 a  JM +1 )   .  
   
     !      Les vents    U       ( > ) correspondent a des valeurs  semi-  
     !       entieres  de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)  
   
     !      Les vents    V       ( V ) correspondent a des valeurs entieres  
     !     de     i ( 1 a  IM ) et semi-entieres de  j ( 1 +0.5  a JM +0.5)  
   
     USE dimens_m, ONLY : iim, jjm  
     USE paramet_m, ONLY : iip1, jjp1  
     USE comconst, ONLY : g, omeg, pi, rad  
     USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh  
     USE logic, ONLY : fxyhypb, ysinus  
55      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &      USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
56           alpha1p2_2d, alpha1p4_2d, alpha1_2d, &           alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
57           alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &           alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
# Line 112  contains Line 61  contains
61           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &           rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
62           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &           unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
63           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv           unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
64        USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
65        USE dimens_m, ONLY : iim, jjm
66        USE logic, ONLY : fxyhypb, ysinus
67        use nr_util, only: pi
68        USE paramet_m, ONLY : iip1, jjp1
69      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &      USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
70           grossismy, pxo, pyo, taux, tauy, transx, transy           grossismy, pxo, pyo, taux, tauy, transx, transy
71    
# Line 139  contains Line 93  contains
93      real aireij2_2d(iim + 1, jjm + 1)      real aireij2_2d(iim + 1, jjm + 1)
94      real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1)      real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1)
95      real airuscv2_2d(iim + 1, jjm)      real airuscv2_2d(iim + 1, jjm)
96      real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
97      real aivscu2gam_2d(iim + 1, jjm)      real aivscu2gam_2d(iim + 1, jjm)
98    
99      !------------------------------------------------------------------      !------------------------------------------------------------------
# Line 168  contains Line 122  contains
122    
123      WRITE (6, 990)      WRITE (6, 990)
124    
125      IF ( .NOT. fxyhypb) THEN      IF (.NOT. fxyhypb) THEN
126         IF (ysinus) THEN         IF (ysinus) THEN
127            print *, ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '            print *, ' Inigeom, Y = Sinus (Latitude) '
128              ! utilisation de f(x, y) avec y = sinus de la latitude
           ! utilisation de f(x, y )  avec  y  =  sinus de la latitude  ...  
   
129            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
130                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
131                 xprimm025, rlonp025, xprimp025)                 xprimm025, rlonp025, xprimp025)
132         ELSE         ELSE
133            print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'            print *, 'Inigeom, Y = Latitude, der. sinusoid .'
134            ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit            ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
135    
136            pxo = clon*pi/180.            pxo = clon*pi/180.
137            pyo = 2.*clat*pi/180.            pyo = 2.*clat*pi/180.
138    
139            ! determination de  transx ( pour le zoom ) par Newton-Raphson .            ! determination de transx (pour le zoom) par Newton-Raphson
140    
141            itmax = 10            itmax = 10
142            eps = .1E-7            eps = .1E-7
# Line 223  contains Line 175  contains
175                 rlonp025, xprimp025)                 rlonp025, xprimp025)
176         END IF         END IF
177      ELSE      ELSE
178         ! ....  Utilisation  de fxyhyper , f(x, y) a derivee tangente hyperbol.         ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
179         print *, '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'         print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
180         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
181              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
182              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
# Line 234  contains Line 186  contains
186      rlatu(1) = asin(1.)      rlatu(1) = asin(1.)
187      rlatu(jjp1) = -rlatu(1)      rlatu(jjp1) = -rlatu(1)
188    
189      !   ....  calcul  aux  poles  ....      ! calcul aux poles
190    
191      yprimu(1) = 0.      yprimu(1) = 0.
192      yprimu(jjp1) = 0.      yprimu(jjp1) = 0.
193    
194      un4rad2 = 0.25*rad*rad      un4rad2 = 0.25*rad*rad
195    
196      ! calcul  des aires ( aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez )      ! calcul des aires (aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez)
197      !   -      et de   fext_2d ,  force de coriolis  extensive  .      ! - et de fext_2d, force de coriolis extensive
198    
199      !   A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y) , sont      ! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y), sont
200      !   affectees 4 aires entourant P , calculees respectivement aux points      ! affectees 4 aires entourant P, calculees respectivement aux points
201      !            ( i + 1/4, j - 1/4 )    :    aireij1_2d (i, j)      ! (i + 1/4, j - 1/4) : aireij1_2d (i, j)
202      !            ( i + 1/4, j + 1/4 )    :    aireij2_2d (i, j)      ! (i + 1/4, j + 1/4) : aireij2_2d (i, j)
203      !            ( i - 1/4, j + 1/4 )    :    aireij3_2d (i, j)      ! (i - 1/4, j + 1/4) : aireij3_2d (i, j)
204      !            ( i - 1/4, j - 1/4 )    :    aireij4_2d (i, j)      ! (i - 1/4, j - 1/4) : aireij4_2d (i, j)
205    
206      !           ,      !,
207      ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).      ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).
208      !   Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme      ! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme
209      ! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont      ! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont
210      ! affectees au      ! affectees au
211      !   point (i, j) .      ! point (i, j).
212      !   On definit en outre les coefficients  alpha comme etant egaux a      ! On definit en outre les coefficients alpha comme etant egaux a
213      ! (aireij / aire_2d), c.a.d par exp.      ! (aireij / aire_2d), c.a.d par exp.
214      ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)      ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)
215    
216      !   De meme, toute aire centree en 1 point U est egale a la somme des      ! De meme, toute aire centree en 1 point U est egale a la somme des
217      ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant      ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
218      ! le point U.      ! le point U.
219      !    Idem pour  airev_2d, airez .      ! Idem pour airev_2d, airez.
220    
221      !       On a , pour chaque maille :    dX = dY = 1      ! On a, pour chaque maille : dX = dY = 1
222    
223      !                             . V      ! V
224    
225      !                 aireij4_2d .        . aireij1_2d      ! aireij4_2d . . aireij1_2d
226    
227      !                   U .       . P      . U      ! U . . P . U
228    
229      !                 aireij3_2d .        . aireij2_2d      ! aireij3_2d . . aireij2_2d
230    
231      !                             . V      ! V
232    
233      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
234      ! aireij3_2d, aireij4_2d      ! aireij3_2d, aireij4_2d
235      ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations      ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations
236      ! elementaires      ! elementaires
237      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      ! cuij et les 4 elongat. cvij qui sont calculees aux memes
238      !     endroits  que les aireij   .      ! endroits que les aireij.
239    
240      !     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....      ! do 35 : boucle sur les jjm + 1 latitudes
241    
242      DO j = 1, jjp1      DO j = 1, jjp1
243    
# Line 382  contains Line 334  contains
334    
335         END IF         END IF
336    
337         !    ........       periodicite   ............         ! periodicite
338    
339         cvij1(iip1, j) = cvij1(1, j)         cvij1(iip1, j) = cvij1(1, j)
340         cvij2(iip1, j) = cvij2(1, j)         cvij2(iip1, j) = cvij2(1, j)
# Line 460  contains Line 412  contains
412    
413      END DO      END DO
414    
415      !    .....      Calcul  des elongations cu_2d, cv_2d, cvu     .........      ! Calcul des elongations cu_2d, cv_2d, cvu
416    
417      DO j = 1, jjm      DO j = 1, jjm
418         DO i = 1, iim         DO i = 1, iim
# Line 508  contains Line 460  contains
460         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
461      END DO      END DO
462    
463      !   ....  calcul aux  poles  ....      ! calcul aux poles
464    
465      DO i = 1, iip1      DO i = 1, iip1
466         cu_2d(i, 1) = 0.         cu_2d(i, 1) = 0.
# Line 538  contains Line 490  contains
490         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
491      END DO      END DO
492    
493      !   calcul des aires aux  poles :      ! calcul des aires aux poles :
494    
495      apoln = sum(aire_2d(:iim, 1))      apoln = sum(aire_2d(:iim, 1))
496      apols = sum(aire_2d(:iim, jjp1))      apols = sum(aire_2d(:iim, jjp1))
# Line 547  contains Line 499  contains
499      unsapolnga2 = 1./(apoln**(-gamdi_h))      unsapolnga2 = 1./(apoln**(-gamdi_h))
500      unsapolsga2 = 1./(apols**(-gamdi_h))      unsapolsga2 = 1./(apols**(-gamdi_h))
501    
502      !   changement F. Hourdin calcul conservatif pour fext_2d      ! changement F. Hourdin calcul conservatif pour fext_2d
503      !   constang_2d contient le produit a * cos ( latitude ) * omega      ! constang_2d contient le produit a * cos (latitude) * omega
504    
505      DO i = 1, iim      DO i = 1, iim
506         constang_2d(i, 1) = 0.         constang_2d(i, 1) = 0.
# Line 562  contains Line 514  contains
514         constang_2d(i, jjp1) = 0.         constang_2d(i, jjp1) = 0.
515      END DO      END DO
516    
517      !   periodicite en longitude      ! periodicite en longitude
518    
519      DO j = 1, jjm      DO j = 1, jjm
520         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
# Line 573  contains Line 525  contains
525    
526      ! fin du changement      ! fin du changement
527    
528      print *, '   ***  Coordonnees de la grille  *** '      print *, ' Coordonnees de la grille '
529      WRITE (6, 995)      WRITE (6, 995)
530    
531      print *, '   LONGITUDES  aux pts.   V  ( degres )  '      print *, ' LONGITUDES aux pts. V (degres) '
532      WRITE (6, 995)      WRITE (6, 995)
533      DO i = 1, iip1      DO i = 1, iip1
534         rlonvv(i) = rlonv(i)*180./pi         rlonvv(i) = rlonv(i)*180./pi
# Line 584  contains Line 536  contains
536      WRITE (6, 400) rlonvv      WRITE (6, 400) rlonvv
537    
538      WRITE (6, 995)      WRITE (6, 995)
539      print *, '   LATITUDES   aux pts.   V  ( degres )  '      print *, ' LATITUDES aux pts. V (degres) '
540      WRITE (6, 995)      WRITE (6, 995)
541      DO i = 1, jjm      DO i = 1, jjm
542         rlatuu(i) = rlatv(i)*180./pi         rlatuu(i) = rlatv(i)*180./pi
# Line 595  contains Line 547  contains
547         rlonvv(i) = rlonu(i)*180./pi         rlonvv(i) = rlonu(i)*180./pi
548      END DO      END DO
549      WRITE (6, 995)      WRITE (6, 995)
550      print *, '   LONGITUDES  aux pts.   U  ( degres )  '      print *, ' LONGITUDES aux pts. U (degres) '
551      WRITE (6, 995)      WRITE (6, 995)
552      WRITE (6, 400) rlonvv      WRITE (6, 400) rlonvv
553      WRITE (6, 995)      WRITE (6, 995)
554    
555      print *, '   LATITUDES   aux pts.   U  ( degres )  '      print *, ' LATITUDES aux pts. U (degres) '
556      WRITE (6, 995)      WRITE (6, 995)
557      DO i = 1, jjp1      DO i = 1, jjp1
558         rlatuu(i) = rlatu(i)*180./pi         rlatuu(i) = rlatu(i)*180./pi
Legend:
Removed from v.25  
changed lines
  Added in v.39
  ViewVC Help
Powered by ViewVC 1.1.21