Changeset 5159 for LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90

Timestamp:

Aug 2, 2024, 9:58:25 PM (7 weeks ago)

Author:

abarral

Message:

Put dimensions.h and paramet.h into modules

File:

• LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90 (modified) (36 diffs)

LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90

@@ -2 +2 @@
 ! $Id$
 
-!
-!
+
+
 SUBROUTINE inigeom
-  !
+
   ! Auteur :  P. Le Van
-  !
+
   !   ............      Version  du 01/04/2001     ........................
-  !
+
   !  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-
   ! endroits que les aires aireij1,..aireij4 .
 
   !  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
-  !
-  !
+
+
   USE fxhyp_m, ONLY: fxhyp
   USE fyhyp_m, ONLY: fyhyp
@@ -27 +27 @@
   USE lmdz_comgeom2
 
+USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
+  USE lmdz_paramet
   IMPLICIT NONE
   !
-  INCLUDE "dimensions.h"
-  INCLUDE "paramet.h"
+
+
 
   !-----------------------------------------------------------------------
   !   ....  Variables  locales   ....
-  !
+
   INTEGER :: i,j,itmax,itmay,iter
   REAL :: cvu(iip1,jjp1),cuv(iip1,jjm)
@@ -53 +55 @@
   SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
   SAVE rlonm025,xprimm025,rlonp025,xprimp025
-  !
-  !
+
+
   !   ------------------------------------------------------------------
   !   -                                                                -
@@ -60 +62 @@
   !   -                                                                -
   !   ------------------------------------------------------------------
-  !
+
   !  les coef. ( cu, cv ) permettent de passer des vitesses naturelles
   !  aux vitesses covariantes et contravariantes , ou vice-versa ...
-  !
-  !
+
+
   ! on a :  u (covariant) = cu * u (naturel)   , u(contrav)= u(nat)/cu
   !         v (covariant) = cv * v (naturel)   , v(contrav)= v(nat)/cv
-  !
+
   !   on en tire :  u(covariant) = cu * cu * u(contravariant)
   !                 v(covariant) = cv * cv * v(contravariant)
-  !
-  !
+
+
   ! on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2
   !                                                      =     =
   !                                       et   - jm/2    <  Y  < jm/2
   !                                                      =     =
-  !
+
   !  ...................................................
   !  ...................................................
@@ -90 +92 @@
   !  ...................................................
   !  ...................................................
-  !
-  !
-  !
+
+
+
   !                                                       ,
   !    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
@@ -104 +106 @@
   !  rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude
   !  cvu    et   cv      .  .  .  .  .  .  .  .  .  .  .    cv
-  !
+
   !  **************  aux points u, v, scalaires, et z  ****************
   !  cu, cuv, cuscal, cuz sont respectiv. les valeurs de    cu
-  !
-  !
-  !
+
+
+
   !     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)
-  !
-  !
-  !
+
+
+
   WRITE(6,3)
  3   FORMAT( // 10x,' ....  INIGEOM  date du 01/06/98   ..... ', &
@@ -169 +171 @@
          5  x,' elong. cuij1, .. 4  , cvij1,.. 4  qui les entourent , aux&
            & '/ 5x,' memes endroits que les aires aireij1,...j4   . ' / )
-  !
-  !
+
+
   IF( nitergdiv/=2 ) THEN
     gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
@@ -189 +191 @@
   WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis, &
         nitergdiv,nitergrot,niterh
-  !
+
   pi    = 2.* ASIN(1.)
-  !
+
   WRITE(6,990)
 
   ! ----------------------------------------------------------------
-  !
+
   IF( .NOT.fxyhypb )   THEN
-  !
-  !
+
+
    IF( ysinus )  THEN
-  !
+
     WRITE(6,*) ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
-  !
+
   !   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
 
@@ -210 +212 @@
 
    ELSE
-  !
+
     WRITE(6,*) '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
 
@@ -218 +220 @@
     pxo   = clon *pi /180.
     pyo   = 2.* clat* pi /180.
-  !
+
   !  ....  determination de  transx ( pour le zoom ) par Newton-Raphson ...
-  !
+
     itmax = 10
     eps   = .1e-7
-  !
+
     xo1 = 0.
     DO iter = 1, itmax
@@ -235 +237 @@
   END DO
  11   CONTINUE
-  !
+
     transx = xo1
 
     itmay = 10
     eps   = .1e-7
-  !
+
     yo1  = 0.
     DO iter = 1,itmay
@@ -251 +253 @@
     yo1  = y1
   END DO
-  !
+
  17   CONTINUE
     transy = yo1
@@ -260 +262 @@
 
    ENDIF
-  !
+
   ELSE
-  !
+
   !   ....  Utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
   !   .....................................................................
@@ -272 +274 @@
 
   ENDIF
-  !
+
   !  -------------------------------------------------------------------
 
-  !
+
   rlatu(1)    =     ASIN(1.)
   rlatu(jjp1) =  - rlatu(1)
-  !
-  !
+
+
   !   ....  calcul  aux  poles  ....
-  !
+
   yprimu(1)      = 0.
   yprimu(jjp1)   = 0.
-  !
-  !
+
+
   un4rad2 = 0.25 * rad * rad
-  !
+
   !   --------------------------------------------------------------------
   !   --------------------------------------------------------------------
@@ -296 +298 @@
   !   --------------------------------------------------------------------
   !   --------------------------------------------------------------------
-  !
-  !
-  !
+
+
+
   !   A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
   !   affectees 4 aires entourant P , calculees respectivement aux points
@@ -305 +307 @@
   !        ( i - 1/4, j + 1/4 )    :    aireij3 (i,j)
   !        ( i - 1/4, j - 1/4 )    :    aireij4 (i,j)
-  !
+
   !       ,
   !   Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
@@ -313 +315 @@
   !   On definit en outre les coefficients  alpha comme etant egaux a
   !    (aireij / aire), c.a.d par exp.  alpha1(i,j)=aireij1(i,j)/aire(i,j)
-  !
+
   !   De meme, toute aire centree en 1 point U est egale a la somme des
   !   4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
   !    Idem pour  airev, airez .
-  !
+
   !   On a ,pour chaque maille :    dX = dY = 1
-  !
-  !
+
+
   !                         . V
-  !
+
   !             aireij4 .        . aireij1
-  !
+
   !               U .       . P      . U
-  !
+
   !             aireij3 .        . aireij2
-  !
+
   !                         . V
-  !
-  !
-  !
-  !
-  !
+
+
+
+
+
   !   ....................................................................
-  !
+
   !    Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
   !   qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
   !   taires cuij et les 4 elongat. cvij qui sont calculees aux memes
   ! endroits  que les aireij   .
-  !
+
   !   ....................................................................
-  !
+
   ! .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....
-  !
-  !
+
+
   DO j = 1, jjp1
-  !
+
   IF ( j == 1 )  THEN
-  !
+
   yprm           = yprimu1(j)
   rlatm          = rlatu1(j)
-  !
+
   coslatm        = COS( rlatm )
   radclatm       = 0.5* rad * coslatm
-  !
+
   DO i = 1, iim
   xprp           = xprimp025( i )
@@ -367 +369 @@
   cvij3  ( i,1 ) = cvij2(i,1)
   END DO
-  !
+
   DO  i = 1, iim
   aireij1( i,1 ) = 0.
@@ -376 +378 @@
   cvij4  ( i,1 ) = 0.
   ENDDO
-  !
+
   END IF
-  !
+
   IF ( j == jjp1 )  THEN
    yprp               = yprimu2(j-1)
@@ -384 +386 @@
   !cc       yprp             = fyprim( REAL(j) - 0.25 )
   !cc       rlatp            = fy    ( REAL(j) - 0.25 )
-  !
+
   coslatp             = COS( rlatp )
   radclatp            = 0.5* rad * coslatp
-  !
+
   DO i = 1,iim
     xprp              = xprimp025( i )
@@ -398 +400 @@
     cvij4(i,jjp1)     = cvij1(i,jjp1)
   END DO
-  !
+
    DO   i    = 1, iim
     aireij2( i,jjp1 ) = 0.
@@ -407 +409 @@
     cuij3  ( i,jjp1 ) = 0.
    ENDDO
-  !
+
   END IF
   !
 
   IF ( j > 1 .AND. j < jjp1 )  THEN
-  !
+
     rlatp    = rlatu2 ( j-1 )
     yprp     = yprimu2( j-1 )
@@ -426 +428 @@
      radclatp = 0.5* rad * coslatp
      radclatm = 0.5* rad * coslatm
-  !
+
      ai14            = un4rad2 * coslatp * yprp
      ai23            = un4rad2 * coslatm * yprm
@@ -446 +448 @@
      cvij4   ( i,j ) = cvij1(i,j)
   END DO
-  !
+
   END IF
-  !
+
   !    ........       periodicite   ............
-  !
+
      cvij1   (iip1,j) = cvij1   (1,j)
      cvij2   (iip1,j) = cvij2   (1,j)
@@ -465 +467 @@
 
   END DO
-  !
+
   !    ..............................................................
-  !
+
   DO j = 1, jjp1
   DO i = 1, iim
@@ -481 +483 @@
   alpha3p4( i,j )  = alpha3 (i,j) + alpha4 (i,j)
   END DO
-  !
-  !
+
+
   aire    (iip1,j) = aire    (1,j)
   alpha1  (iip1,j) = alpha1  (1,j)
@@ -510 +512 @@
   airesurg   (iip1,j) = airesurg(1,j)
   END DO
-  !
-  !
+
+
   DO j = 1,jjm
-  !
+
     DO i=1,iim
      airev     (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) + &
@@ -529 +531 @@
     fext      (iip1,j)  = fext(1,j)
     unsairz_gam(iip1,j) = unsairz_gam(1,j)
-  !
-  END DO
-  !
-  !
+
+  END DO
+
+
   !    .....      Calcul  des elongations cu,cv, cvu     .........
-  !
+
   DO    j   = 1, jjm
    DO   i  = 1, iim
@@ -579 +581 @@
   ENDDO
 
-  !
+
   !   ....  calcul aux  poles  ....
-  !
+
   DO    i      =  1, iip1
     cu    ( i, 1 )  =   0.
     unscu2( i, 1 )  =   0.
     cvu   ( i, 1 )  =   0.
-  !
+
     cu    (i, jjp1) =   0.
     unscu2(i, jjp1) =   0.
     cvu   (i, jjp1) =   0.
   ENDDO
-  !
+
   !    ..............................................................
-  !
+
   DO j = 1, jjm
     DO i= 1, iim
@@ -612 +614 @@
   ENDDO
 
-  !
+
   !   calcul des aires aux  poles :
   !   -----------------------------
-  !
+
   apoln       = SSUM(iim,aire(1,1),1)
   apols       = SSUM(iim,aire(1,jjp1),1)
@@ -622 +624 @@
   unsapolnga2 = 1./ ( apoln ** ( - gamdi_h    ) )
   unsapolsga2 = 1./ ( apols ** ( - gamdi_h    ) )
-  !
+
   !-----------------------------------------------------------------------
   ! gtitre='Coriolis version ancienne'
   ! gfichier='fext1'
   ! CALL writestd(fext,iip1*jjm)
-  !
+
   !   changement F. Hourdin calcul conservatif pour fext
   !   constang contient le produit a * cos ( latitude ) * omega
-  !
+
   DO i=1,iim
      constang(i,1) = 0.
@@ -642 +644 @@
      constang(i,jjp1) = 0.
   ENDDO
-  !
+
   !   periodicite en longitude
-  !
+
   DO j=1,jjm
     fext(iip1,j)     = fext(1,j)
@@ -654 +656 @@
   ! fin du changement
 
-  !
+
   !-----------------------------------------------------------------------
-  !
+
    WRITE(6,*) '   ***  Coordonnees de la grille  *** '
    WRITE(6,995)
-  !
+
    WRITE(6,*) '   LONGITUDES  aux pts.   V  ( degres )  '
    WRITE(6,995)
@@ -666 +668 @@
     ENDDO
    WRITE(6,400) rlonvv
-  !
+
    WRITE(6,995)
    WRITE(6,*) '   LATITUDES   aux pts.   V  ( degres )  '
@@ -674 +676 @@
     ENDDO
    WRITE(6,400) (rlatuu(i),i=1,jjm)
-  !
+
     DO i=1,iip1
       rlonvv(i)=rlonu(i)*180./pi
@@ -691 +693 @@
    WRITE(6,400) (rlatuu(i),i=1,jjp1)
    WRITE(6,995)
-  !
+
 444   format(f10.3,f6.0)
 400   FORMAT(1x,8f8.2)
 990   FORMAT(//)
 995   FORMAT(/)
-  !
+
   RETURN
 END SUBROUTINE inigeom