/[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 37 by guez, Fri Mar 5 16:43:45 2010 UTC revision 38 by guez, Thu Jan 6 17:52:19 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      !       on en tire :  u(covariant) = cu_2d * cu_2d * u(contravariant)  
26      !                     v(covariant) = cv_2d * cv_2d * v(contravariant)      ! on en tire :
27        ! u(covariant) = cu_2d * cu_2d * u(contravariant)
28      !     on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2      ! v(covariant) = cv_2d * cv_2d * v(contravariant)
29      !                                                          =     =  
30      !                                           et   - jm/2    <  Y  < jm/2      ! 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       .      ! x est la longitude du point en radians.
34      !      .  y  est la  latitude du point  en radians       .      ! y est la latitude du point en radians.
35      !      .                                                 .      !
36      !      .  on a :  cu_2d(i, j) = rad * COS(y) * dx/dX         .      ! on a : cu_2d(i, j) = rad * cos(y) * dx/dX
37      !      .          cv( j ) = rad          * dy/dY         .      ! cv(j) = rad * dy/dY
38      !      .        aire_2d(i, j) =  cu_2d(i, j) * cv(j)             .      ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
39      !      .                                                 .      !
40      !      . y, dx/dX, dy/dY calcules aux points concernes   .      ! y, dx/dX, dy/dY calcules aux points concernes
41      !                                                           ,      ! cv, bien que dependant de j uniquement, sera ici indice aussi en i
42      !    cv , bien que dependant de j uniquement, sera ici indice aussi en i      ! pour un adressage plus facile en ij.
43      !          pour un adressage plus facile en  ij  .  
44        ! aux points u et v,
45      !  **************  aux points  u  et  v ,           *****************      ! xprimu et xprimv sont respectivement les valeurs de dx/dX
46      !      xprimu et xprimv sont respectivement les valeurs de  dx/dX      ! yprimu et yprimv sont respectivement les valeurs de dy/dY
47      !      yprimu et yprimv    .  .  .  .  .  .  .  .  .  .  .  dy/dY      ! rlatu et rlatv sont respectivement les valeurs de la latitude
48      !      rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude      ! cvu et cv_2d sont respectivement les valeurs de cv_2d
49      !      cvu    et   cv_2d      .  .  .  .  .  .  .  .  .  .  .    cv_2d  
50        ! aux points u, v, scalaires, et z
51      !  **************  aux points u, v, scalaires, et z  ****************      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
52      !      cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de    cu_2d      ! Cf. "inigeom.txt".
   
     !         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)  
53    
54      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
55      USE paramet_m, ONLY : iip1, jjp1      USE paramet_m, ONLY : iip1, jjp1
# Line 139  contains Line 92  contains
92      real aireij2_2d(iim + 1, jjm + 1)      real aireij2_2d(iim + 1, jjm + 1)
93      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)
94      real airuscv2_2d(iim + 1, jjm)      real airuscv2_2d(iim + 1, jjm)
95      real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)        real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
96      real aivscu2gam_2d(iim + 1, jjm)      real aivscu2gam_2d(iim + 1, jjm)
97    
98      !------------------------------------------------------------------      !------------------------------------------------------------------
# Line 168  contains Line 121  contains
121    
122      WRITE (6, 990)      WRITE (6, 990)
123    
124      IF ( .NOT. fxyhypb) THEN      IF (.NOT. fxyhypb) THEN
125         IF (ysinus) THEN         IF (ysinus) THEN
126            print *, ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '            print *, ' Inigeom, Y = Sinus (Latitude) '
127              ! utilisation de f(x, y) avec y = sinus de la latitude
           ! utilisation de f(x, y )  avec  y  =  sinus de la latitude  ...  
   
128            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &            CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
129                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &                 rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
130                 xprimm025, rlonp025, xprimp025)                 xprimm025, rlonp025, xprimp025)
131         ELSE         ELSE
132            print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'            print *, 'Inigeom, Y = Latitude, der. sinusoid .'
133            ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit            ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
134    
135            pxo = clon*pi/180.            pxo = clon*pi/180.
136            pyo = 2.*clat*pi/180.            pyo = 2.*clat*pi/180.
137    
138            ! determination de  transx ( pour le zoom ) par Newton-Raphson .            ! determination de transx (pour le zoom) par Newton-Raphson
139    
140            itmax = 10            itmax = 10
141            eps = .1E-7            eps = .1E-7
# Line 223  contains Line 174  contains
174                 rlonp025, xprimp025)                 rlonp025, xprimp025)
175         END IF         END IF
176      ELSE      ELSE
177         ! ....  Utilisation  de fxyhyper , f(x, y) a derivee tangente hyperbol.         ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
178         print *, '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'         print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
179         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &         CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
180              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &              taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
181              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &              yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
# Line 234  contains Line 185  contains
185      rlatu(1) = asin(1.)      rlatu(1) = asin(1.)
186      rlatu(jjp1) = -rlatu(1)      rlatu(jjp1) = -rlatu(1)
187    
188      !   ....  calcul  aux  poles  ....      ! calcul aux poles
189    
190      yprimu(1) = 0.      yprimu(1) = 0.
191      yprimu(jjp1) = 0.      yprimu(jjp1) = 0.
192    
193      un4rad2 = 0.25*rad*rad      un4rad2 = 0.25*rad*rad
194    
195      ! 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)
196      !   -      et de   fext_2d ,  force de coriolis  extensive  .      ! - et de fext_2d, force de coriolis extensive
197    
198      !   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
199      !   affectees 4 aires entourant P , calculees respectivement aux points      ! affectees 4 aires entourant P, calculees respectivement aux points
200      !            ( i + 1/4, j - 1/4 )    :    aireij1_2d (i, j)      ! (i + 1/4, j - 1/4) : aireij1_2d (i, j)
201      !            ( i + 1/4, j + 1/4 )    :    aireij2_2d (i, j)      ! (i + 1/4, j + 1/4) : aireij2_2d (i, j)
202      !            ( i - 1/4, j + 1/4 )    :    aireij3_2d (i, j)      ! (i - 1/4, j + 1/4) : aireij3_2d (i, j)
203      !            ( i - 1/4, j - 1/4 )    :    aireij4_2d (i, j)      ! (i - 1/4, j - 1/4) : aireij4_2d (i, j)
204    
205      !           ,      !,
206      ! 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).
207      !   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
208      ! 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
209      ! affectees au      ! affectees au
210      !   point (i, j) .      ! point (i, j).
211      !   On definit en outre les coefficients  alpha comme etant egaux a      ! On definit en outre les coefficients alpha comme etant egaux a
212      ! (aireij / aire_2d), c.a.d par exp.      ! (aireij / aire_2d), c.a.d par exp.
213      ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)      ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)
214    
215      !   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
216      ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant      ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
217      ! le point U.      ! le point U.
218      !    Idem pour  airev_2d, airez .      ! Idem pour airev_2d, airez.
219    
220      !       On a , pour chaque maille :    dX = dY = 1      ! On a, pour chaque maille : dX = dY = 1
221    
222      !                             . V      ! V
223    
224      !                 aireij4_2d .        . aireij1_2d      ! aireij4_2d . . aireij1_2d
225    
226      !                   U .       . P      . U      ! U . . P . U
227    
228      !                 aireij3_2d .        . aireij2_2d      ! aireij3_2d . . aireij2_2d
229    
230      !                             . V      ! V
231    
232      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,      ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
233      ! aireij3_2d, aireij4_2d      ! aireij3_2d, aireij4_2d
234      ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations      ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations
235      ! elementaires      ! elementaires
236      ! cuij et les 4 elongat. cvij qui sont calculees aux memes      ! cuij et les 4 elongat. cvij qui sont calculees aux memes
237      !     endroits  que les aireij   .      ! endroits que les aireij.
238    
239      !     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....      ! do 35 : boucle sur les jjm + 1 latitudes
240    
241      DO j = 1, jjp1      DO j = 1, jjp1
242    
# Line 382  contains Line 333  contains
333    
334         END IF         END IF
335    
336         !    ........       periodicite   ............         ! periodicite
337    
338         cvij1(iip1, j) = cvij1(1, j)         cvij1(iip1, j) = cvij1(1, j)
339         cvij2(iip1, j) = cvij2(1, j)         cvij2(iip1, j) = cvij2(1, j)
# Line 460  contains Line 411  contains
411    
412      END DO      END DO
413    
414      !    .....      Calcul  des elongations cu_2d, cv_2d, cvu     .........      ! Calcul des elongations cu_2d, cv_2d, cvu
415    
416      DO j = 1, jjm      DO j = 1, jjm
417         DO i = 1, iim         DO i = 1, iim
# Line 508  contains Line 459  contains
459         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)         cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
460      END DO      END DO
461    
462      !   ....  calcul aux  poles  ....      ! calcul aux poles
463    
464      DO i = 1, iip1      DO i = 1, iip1
465         cu_2d(i, 1) = 0.         cu_2d(i, 1) = 0.
# Line 538  contains Line 489  contains
489         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)         aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
490      END DO      END DO
491    
492      !   calcul des aires aux  poles :      ! calcul des aires aux poles :
493    
494      apoln = sum(aire_2d(:iim, 1))      apoln = sum(aire_2d(:iim, 1))
495      apols = sum(aire_2d(:iim, jjp1))      apols = sum(aire_2d(:iim, jjp1))
# Line 547  contains Line 498  contains
498      unsapolnga2 = 1./(apoln**(-gamdi_h))      unsapolnga2 = 1./(apoln**(-gamdi_h))
499      unsapolsga2 = 1./(apols**(-gamdi_h))      unsapolsga2 = 1./(apols**(-gamdi_h))
500    
501      !   changement F. Hourdin calcul conservatif pour fext_2d      ! changement F. Hourdin calcul conservatif pour fext_2d
502      !   constang_2d contient le produit a * cos ( latitude ) * omega      ! constang_2d contient le produit a * cos (latitude) * omega
503    
504      DO i = 1, iim      DO i = 1, iim
505         constang_2d(i, 1) = 0.         constang_2d(i, 1) = 0.
# Line 562  contains Line 513  contains
513         constang_2d(i, jjp1) = 0.         constang_2d(i, jjp1) = 0.
514      END DO      END DO
515    
516      !   periodicite en longitude      ! periodicite en longitude
517    
518      DO j = 1, jjm      DO j = 1, jjm
519         fext_2d(iip1, j) = fext_2d(1, j)         fext_2d(iip1, j) = fext_2d(1, j)
# Line 573  contains Line 524  contains
524    
525      ! fin du changement      ! fin du changement
526    
527      print *, '   ***  Coordonnees de la grille  *** '      print *, ' Coordonnees de la grille '
528      WRITE (6, 995)      WRITE (6, 995)
529    
530      print *, '   LONGITUDES  aux pts.   V  ( degres )  '      print *, ' LONGITUDES aux pts. V (degres) '
531      WRITE (6, 995)      WRITE (6, 995)
532      DO i = 1, iip1      DO i = 1, iip1
533         rlonvv(i) = rlonv(i)*180./pi         rlonvv(i) = rlonv(i)*180./pi
# Line 584  contains Line 535  contains
535      WRITE (6, 400) rlonvv      WRITE (6, 400) rlonvv
536    
537      WRITE (6, 995)      WRITE (6, 995)
538      print *, '   LATITUDES   aux pts.   V  ( degres )  '      print *, ' LATITUDES aux pts. V (degres) '
539      WRITE (6, 995)      WRITE (6, 995)
540      DO i = 1, jjm      DO i = 1, jjm
541         rlatuu(i) = rlatv(i)*180./pi         rlatuu(i) = rlatv(i)*180./pi
# Line 595  contains Line 546  contains
546         rlonvv(i) = rlonu(i)*180./pi         rlonvv(i) = rlonu(i)*180./pi
547      END DO      END DO
548      WRITE (6, 995)      WRITE (6, 995)
549      print *, '   LONGITUDES  aux pts.   U  ( degres )  '      print *, ' LONGITUDES aux pts. U (degres) '
550      WRITE (6, 995)      WRITE (6, 995)
551      WRITE (6, 400) rlonvv      WRITE (6, 400) rlonvv
552      WRITE (6, 995)      WRITE (6, 995)
553    
554      print *, '   LATITUDES   aux pts.   U  ( degres )  '      print *, ' LATITUDES aux pts. U (degres) '
555      WRITE (6, 995)      WRITE (6, 995)
556      DO i = 1, jjp1      DO i = 1, jjp1
557         rlatuu(i) = rlatu(i)*180./pi         rlatuu(i) = rlatu(i)*180./pi
Legend:
Removed from v.37  
changed lines
  Added in v.38
  ViewVC Help
Powered by ViewVC 1.1.21