libf/dyn3d/inigeom.f90

module inigeom_m

  IMPLICIT NONE

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van
    ! Version du 01/04/2001

    ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
    ! endroits que les aires aireij1_2d, ..., aireij4_2d.

    ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée
    ! tangente hyperbolique
    ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2)

    ! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles
    !      aux vitesses covariantes et contravariantes , ou vice-versa ...

    ! on a :  u (covariant) = cu_2d * u (naturel)   , u(contrav)= u(nat)/cu_2d
    !         v (covariant) = cv_2d * v (naturel)   , v(contrav)= v(nat)/cv_2d

    !       on en tire :  u(covariant) = cu_2d * cu_2d * u(contravariant)
    !                     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
    USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
         alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
         alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
         apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &
         cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &
         cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &
         rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
         unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
         unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
    USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
         grossismy, pxo, pyo, taux, tauy, transx, transy

    ! Variables locales

    INTEGER i, j, itmax, itmay, iter
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm
    REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
    REAL coslatm, coslatp, radclatm, radclatp
    REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &
         cuij4(iip1, jjp1)
    REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &
         cvij4(iip1, jjp1)
    REAL rlonvv(iip1), rlatuu(jjp1)
    REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), &
         yprimu(jjp1)
    REAL gamdi_gdiv, gamdi_grot, gamdi_h

    REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
    SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu
    SAVE rlonm025, xprimm025, rlonp025, xprimp025

    real aireij1_2d(iim + 1, jjm + 1)
    real aireij2_2d(iim + 1, jjm + 1)
    real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1) 
    real airuscv2_2d(iim + 1, jjm) 
    real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)   
    real aivscu2gam_2d(iim + 1, jjm)

    !------------------------------------------------------------------

    PRINT *, 'Call sequence information: inigeom'

    IF (nitergdiv/=2) THEN
       gamdi_gdiv = coefdis/(real(nitergdiv)-2.)
    ELSE
       gamdi_gdiv = 0.
    END IF
    IF (nitergrot/=2) THEN
       gamdi_grot = coefdis/(real(nitergrot)-2.)
    ELSE
       gamdi_grot = 0.
    END IF
    IF (niterh/=2) THEN
       gamdi_h = coefdis/(real(niterh)-2.)
    ELSE
       gamdi_h = 0.
    END IF

    print *, 'gamdi_gdiv = ', gamdi_gdiv
    print *, "gamdi_grot = ", gamdi_grot
    print *, "gamdi_h = ", gamdi_h

    WRITE (6, 990)

    IF ( .NOT. fxyhypb) THEN
       IF (ysinus) THEN
          print *, ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '

          ! utilisation de f(x, y )  avec  y  =  sinus de la latitude  ...

          CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
               rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
               xprimm025, rlonp025, xprimp025)
       ELSE
          print *, '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
          ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit

          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
             x1 = xo1
             f = x1 + alphax*sin(x1-pxo)
             df = 1. + alphax*cos(x1-pxo)
             x1 = x1 - f/df
             xdm = abs(x1-xo1)
             IF (xdm<=eps) EXIT
             xo1 = x1
          END DO

          transx = xo1

          itmay = 10
          eps = .1E-7

          yo1 = 0.
          DO iter = 1, itmay
             y1 = yo1
             f = y1 + alphay*sin(y1-pyo)
             df = 1. + alphay*cos(y1-pyo)
             y1 = y1 - f/df
             ydm = abs(y1-yo1)
             IF (ydm<=eps) EXIT
             yo1 = y1
          END DO

          transy = yo1

          CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
               yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
               rlonp025, xprimp025)
       END IF
    ELSE
       ! ....  Utilisation  de fxyhyper , f(x, y) a derivee tangente hyperbol.
       print *, '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
       CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
            taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
            yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
            rlonp025, xprimp025)
    END IF

    rlatu(1) = asin(1.)
    rlatu(jjp1) = -rlatu(1)

    !   ....  calcul  aux  poles  ....

    yprimu(1) = 0.
    yprimu(jjp1) = 0.

    un4rad2 = 0.25*rad*rad

    ! calcul  des aires ( aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez )
    !   -      et de   fext_2d ,  force de coriolis  extensive  .

    !   A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y) , sont
    !   affectees 4 aires entourant P , calculees respectivement aux points
    !            ( i + 1/4, j - 1/4 )    :    aireij1_2d (i, j)
    !            ( i + 1/4, j + 1/4 )    :    aireij2_2d (i, j)
    !            ( i - 1/4, j + 1/4 )    :    aireij3_2d (i, j)
    !            ( i - 1/4, j - 1/4 )    :    aireij4_2d (i, j)

    !           ,
    ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).
    !   Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme
    ! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont
    ! affectees au
    !   point (i, j) .
    !   On definit en outre les coefficients  alpha comme etant egaux a
    ! (aireij / aire_2d), c.a.d par exp.
    ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)

    !   De meme, toute aire centree en 1 point U est egale a la somme des
    ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
    ! le point U.
    !    Idem pour  airev_2d, airez .

    !       On a , pour chaque maille :    dX = dY = 1

    !                             . V

    !                 aireij4_2d .        . aireij1_2d

    !                   U .       . P      . U

    !                 aireij3_2d .        . aireij2_2d

    !                             . V

    ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
    ! aireij3_2d, aireij4_2d
    ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations
    ! elementaires
    ! 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)
             xprm = xprimm025(i)
             aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm
             aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm
             cuij2(i, 1) = radclatm*xprp
             cuij3(i, 1) = radclatm*xprm
             cvij2(i, 1) = 0.5*rad*yprm
             cvij3(i, 1) = cvij2(i, 1)
          END DO

          DO i = 1, iim
             aireij1_2d(i, 1) = 0.
             aireij4_2d(i, 1) = 0.
             cuij1(i, 1) = 0.
             cuij4(i, 1) = 0.
             cvij1(i, 1) = 0.
             cvij4(i, 1) = 0.
          END DO

       END IF

       IF (j==jjp1) THEN
          yprp = yprimu2(j-1)
          rlatp = rlatu2(j-1)

          coslatp = cos(rlatp)
          radclatp = 0.5*rad*coslatp

          DO i = 1, iim
             xprp = xprimp025(i)
             xprm = xprimm025(i)
             aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp
             aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp
             cuij1(i, jjp1) = radclatp*xprp
             cuij4(i, jjp1) = radclatp*xprm
             cvij1(i, jjp1) = 0.5*rad*yprp
             cvij4(i, jjp1) = cvij1(i, jjp1)
          END DO

          DO i = 1, iim
             aireij2_2d(i, jjp1) = 0.
             aireij3_2d(i, jjp1) = 0.
             cvij2(i, jjp1) = 0.
             cvij3(i, jjp1) = 0.
             cuij2(i, jjp1) = 0.
             cuij3(i, jjp1) = 0.
          END DO

       END IF

       IF (j>1 .AND. j<jjp1) THEN

          rlatp = rlatu2(j-1)
          yprp = yprimu2(j-1)
          rlatm = rlatu1(j)
          yprm = yprimu1(j)

          coslatm = cos(rlatm)
          coslatp = cos(rlatp)
          radclatp = 0.5*rad*coslatp
          radclatm = 0.5*rad*coslatm

          DO i = 1, iim
             xprp = xprimp025(i)
             xprm = xprimm025(i)

             ai14 = un4rad2*coslatp*yprp
             ai23 = un4rad2*coslatm*yprm
             aireij1_2d(i, j) = ai14*xprp
             aireij2_2d(i, j) = ai23*xprp
             aireij3_2d(i, j) = ai23*xprm
             aireij4_2d(i, j) = ai14*xprm
             cuij1(i, j) = radclatp*xprp
             cuij2(i, j) = radclatm*xprp
             cuij3(i, j) = radclatm*xprm
             cuij4(i, j) = radclatp*xprm
             cvij1(i, j) = 0.5*rad*yprp
             cvij2(i, j) = 0.5*rad*yprm
             cvij3(i, j) = cvij2(i, j)
             cvij4(i, j) = cvij1(i, j)
          END DO

       END IF

       !    ........       periodicite   ............

       cvij1(iip1, j) = cvij1(1, j)
       cvij2(iip1, j) = cvij2(1, j)
       cvij3(iip1, j) = cvij3(1, j)
       cvij4(iip1, j) = cvij4(1, j)
       cuij1(iip1, j) = cuij1(1, j)
       cuij2(iip1, j) = cuij2(1, j)
       cuij3(iip1, j) = cuij3(1, j)
       cuij4(iip1, j) = cuij4(1, j)
       aireij1_2d(iip1, j) = aireij1_2d(1, j)
       aireij2_2d(iip1, j) = aireij2_2d(1, j)
       aireij3_2d(iip1, j) = aireij3_2d(1, j)
       aireij4_2d(iip1, j) = aireij4_2d(1, j)

    END DO

    DO j = 1, jjp1
       DO i = 1, iim
          aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
               + aireij3_2d(i, j) + aireij4_2d(i, j)
          alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)
          alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)
          alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)
          alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)
          alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
          alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
          alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
          alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
       END DO

       aire_2d(iip1, j) = aire_2d(1, j)
       alpha1_2d(iip1, j) = alpha1_2d(1, j)
       alpha2_2d(iip1, j) = alpha2_2d(1, j)
       alpha3_2d(iip1, j) = alpha3_2d(1, j)
       alpha4_2d(iip1, j) = alpha4_2d(1, j)
       alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
       alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
       alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
       alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
    END DO

    DO j = 1, jjp1
       DO i = 1, iim
          aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
               aireij4_2d(i+1, j) + aireij3_2d(i+1, j)
          unsaire_2d(i, j) = 1./aire_2d(i, j)
          unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
          unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
          airesurg_2d(i, j) = aire_2d(i, j)/g
       END DO
       aireu_2d(iip1, j) = aireu_2d(1, j)
       unsaire_2d(iip1, j) = unsaire_2d(1, j)
       unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
       unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
       airesurg_2d(iip1, j) = airesurg_2d(1, j)
    END DO

    DO j = 1, jjm

       DO i = 1, iim
          airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
               aireij1_2d(i, j+1) + aireij4_2d(i, j+1)
       END DO
       DO i = 1, iim
          airez = aireij2_2d(i, j) + aireij1_2d(i, j+1) + aireij3_2d(i+1, j) &
               + aireij4_2d(i+1, j+1)
          unsairez_2d(i, j) = 1./airez
          unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
          fext_2d(i, j) = airez*sin(rlatv(j))*2.*omeg
       END DO
       airev_2d(iip1, j) = airev_2d(1, j)
       unsairez_2d(iip1, j) = unsairez_2d(1, j)
       fext_2d(iip1, j) = fext_2d(1, j)
       unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)

    END DO

    !    .....      Calcul  des elongations cu_2d, cv_2d, cvu     .........

    DO j = 1, jjm
       DO i = 1, iim
          cv_2d(i, j) = 0.5 * &
               (cvij2(i, j) + cvij3(i, j) + cvij1(i, j+1) + cvij4(i, j+1))
          cvu(i, j) = 0.5*(cvij1(i, j)+cvij4(i, j)+cvij2(i, j)+cvij3(i, j))
          cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))
          unscv2_2d(i, j) = 1./(cv_2d(i, j)*cv_2d(i, j))
       END DO
       DO i = 1, iim
          cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)
          cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)
          cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
          cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
          cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
       END DO
       cv_2d(iip1, j) = cv_2d(1, j)
       cvu(iip1, j) = cvu(1, j)
       unscv2_2d(iip1, j) = unscv2_2d(1, j)
       cuv(iip1, j) = cuv(1, j)
       cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
       cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
       cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
       cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
       cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i+1, j) + cuij2(i, j) &
               + cuij3(i+1, j))
          unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))
          cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)
          cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)
          cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
          cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
          cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
       END DO
       cu_2d(iip1, j) = cu_2d(1, j)
       unscu2_2d(iip1, j) = unscu2_2d(1, j)
       cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
       cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
       cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
       cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
       cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
    END DO

    !   ....  calcul aux  poles  ....

    DO i = 1, iip1
       cu_2d(i, 1) = 0.
       unscu2_2d(i, 1) = 0.
       cvu(i, 1) = 0.

       cu_2d(i, jjp1) = 0.
       unscu2_2d(i, jjp1) = 0.
       cvu(i, jjp1) = 0.
    END DO

    DO j = 1, jjm
       DO i = 1, iim
          airvscu2_2d(i, j) = airev_2d(i, j)/(cuv(i, j)*cuv(i, j))
          aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
       END DO
       airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
       aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          airuscv2_2d(i, j) = aireu_2d(i, j)/(cvu(i, j)*cvu(i, j))
          aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
       END DO
       airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
       aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
    END DO

    !   calcul des aires aux  poles :

    apoln = sum(aire_2d(:iim, 1))
    apols = sum(aire_2d(:iim, jjp1))
    unsapolnga1 = 1./(apoln**(-gamdi_gdiv))
    unsapolsga1 = 1./(apols**(-gamdi_gdiv))
    unsapolnga2 = 1./(apoln**(-gamdi_h))
    unsapolsga2 = 1./(apols**(-gamdi_h))

    !   changement F. Hourdin calcul conservatif pour fext_2d
    !   constang_2d contient le produit a * cos ( latitude ) * omega

    DO i = 1, iim
       constang_2d(i, 1) = 0.
    END DO
    DO j = 1, jjm - 1
       DO i = 1, iim
          constang_2d(i, j+1) = rad*omeg*cu_2d(i, j+1)*cos(rlatu(j+1))
       END DO
    END DO
    DO i = 1, iim
       constang_2d(i, jjp1) = 0.
    END DO

    !   periodicite en longitude

    DO j = 1, jjm
       fext_2d(iip1, j) = fext_2d(1, j)
    END DO
    DO j = 1, jjp1
       constang_2d(iip1, j) = constang_2d(1, j)
    END DO

    ! fin du changement

    print *, '   ***  Coordonnees de la grille  *** '
    WRITE (6, 995)

    print *, '   LONGITUDES  aux pts.   V  ( degres )  '
    WRITE (6, 995)
    DO i = 1, iip1
       rlonvv(i) = rlonv(i)*180./pi
    END DO
    WRITE (6, 400) rlonvv

    WRITE (6, 995)
    print *, '   LATITUDES   aux pts.   V  ( degres )  '
    WRITE (6, 995)
    DO i = 1, jjm
       rlatuu(i) = rlatv(i)*180./pi
    END DO
    WRITE (6, 400) (rlatuu(i), i=1, jjm)

    DO i = 1, iip1
       rlonvv(i) = rlonu(i)*180./pi
    END DO
    WRITE (6, 995)
    print *, '   LONGITUDES  aux pts.   U  ( degres )  '
    WRITE (6, 995)
    WRITE (6, 400) rlonvv
    WRITE (6, 995)

    print *, '   LATITUDES   aux pts.   U  ( degres )  '
    WRITE (6, 995)
    DO i = 1, jjp1
       rlatuu(i) = rlatu(i)*180./pi
    END DO
    WRITE (6, 400) (rlatuu(i), i=1, jjp1)
    WRITE (6, 995)

400 FORMAT (1X, 8F8.2)
990 FORMAT (//)
995 FORMAT (/)

  END SUBROUTINE inigeom

end module inigeom_m
1	guez	25	module inigeom_m
2	guez	3
3	guez	7	IMPLICIT NONE
4	guez	3
5	guez	25	contains
6	guez	7
7	guez	25	SUBROUTINE inigeom
8	guez	7
9	guez	25	! Auteur : P. Le Van
10			! Version du 01/04/2001
11	guez	7
12	guez	25	! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
13			! endroits que les aires aireij1_2d, ..., aireij4_2d.
14	guez	7
15	guez	25	! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée
16			! tangente hyperbolique
17			! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d2, 1./cv_2d2)
18	guez	7
19	guez	25	! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles
20			! aux vitesses covariantes et contravariantes , ou vice-versa ...
21	guez	7
22	guez	25	! on a : u (covariant) = cu_2d * u (naturel) , u(contrav)= u(nat)/cu_2d
23			! v (covariant) = cv_2d * v (naturel) , v(contrav)= v(nat)/cv_2d
24	guez	7
25	guez	25	! on en tire : u(covariant) = cu_2d * cu_2d * u(contravariant)
26			! v(covariant) = cv_2d * cv_2d * v(contravariant)
27	guez	7
28	guez	25	! on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
29			! = =
30			! et - jm/2 < Y < jm/2
31			! = =
32	guez	7
33	guez	25	! . 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			! ,
42			! cv , bien que dependant de j uniquement, sera ici indice aussi en i
43			! pour un adressage plus facile en ij .
44	guez	7
45	guez	25	! ************ aux points u et v , ***************
46			! xprimu et xprimv sont respectivement les valeurs de dx/dX
47			! yprimu et yprimv . . . . . . . . . . . dy/dY
48			! rlatu et rlatv . . . . . . . . . . .la latitude
49			! cvu et cv_2d . . . . . . . . . . . cv_2d
50	guez	7
51	guez	25	! ************ aux points u, v, scalaires, et z **************
52			! cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de cu_2d
53	guez	7
54	guez	25	! Exemple de distribution de variables sur la grille dans le
55			! domaine de travail ( X, Y ) .
56			! DX=DY= 1
57	guez	7
58	guez	25	! + represente un point scalaire ( p.exp la pression )
59			! > represente la composante zonale du vent
60			! V represente la composante meridienne du vent
61			! o represente la vorticite
62	guez	7
63	guez	25	! ---- , car aux poles , les comp.zonales covariantes sont nulles
64	guez	7
65	guez	25	! i ->
66	guez	7
67	guez	25	! 1 2 3 4 5 6 7 8
68			! j
69			! v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
70	guez	7
71	guez	25	! V o V o V o V o V o V o V o V o
72	guez	7
73	guez	25	! 2 + > + > + > + > + > + > + > + >
74	guez	7
75	guez	25	! V o V o V o V o V o V o V o V o
76	guez	7
77	guez	25	! 3 + > + > + > + > + > + > + > + >
78	guez	7
79	guez	25	! V o V o V o V o V o V o V o V o
80	guez	7
81	guez	25	! 4 + > + > + > + > + > + > + > + >
82	guez	7
83	guez	25	! V o V o V o V o V o V o V o V o
84	guez	7
85	guez	25	! 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
86	guez	7
87	guez	25	! Ci-dessus, on voit que le nombre de pts.en longitude est egal
88			! a IM = 8
89			! De meme , le nombre d'intervalles entre les 2 poles est egal
90			! a JM = 4
91	guez	7
92	guez	25	! Les points scalaires ( + ) correspondent donc a des valeurs
93			! entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
94	guez	7
95	guez	25	! Les vents U ( > ) correspondent a des valeurs semi-
96			! entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
97	guez	7
98	guez	25	! Les vents V ( V ) correspondent a des valeurs entieres
99			! de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
100	guez	7
101	guez	25	USE dimens_m, ONLY : iim, jjm
102			USE paramet_m, ONLY : iip1, jjp1
103			USE comconst, ONLY : g, omeg, pi, rad
104			USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
105			USE logic, ONLY : fxyhypb, ysinus
106			USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
107			alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
108			alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
109			apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &
110			cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &
111			cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &
112			rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
113			unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
114			unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
115			USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
116			grossismy, pxo, pyo, taux, tauy, transx, transy
117	guez	7
118	guez	25	! Variables locales
119	guez	7
120	guez	25	INTEGER i, j, itmax, itmay, iter
121			REAL cvu(iip1, jjp1), cuv(iip1, jjm)
122			REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm
123			REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
124			REAL coslatm, coslatp, radclatm, radclatp
125			REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &
126			cuij4(iip1, jjp1)
127			REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &
128			cvij4(iip1, jjp1)
129			REAL rlonvv(iip1), rlatuu(jjp1)
130			REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), &
131			yprimu(jjp1)
132			REAL gamdi_gdiv, gamdi_grot, gamdi_h
133	guez	7
134	guez	25	REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
135			SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu
136			SAVE rlonm025, xprimm025, rlonp025, xprimp025
137	guez	7
138	guez	25	real aireij1_2d(iim + 1, jjm + 1)
139			real aireij2_2d(iim + 1, jjm + 1)
140			real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1)
141			real airuscv2_2d(iim + 1, jjm)
142			real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
143			real aivscu2gam_2d(iim + 1, jjm)
144	guez	7
145	guez	25	!------------------------------------------------------------------
146	guez	7
147	guez	25	PRINT *, 'Call sequence information: inigeom'
148	guez	7
149	guez	25	IF (nitergdiv/=2) THEN
150			gamdi_gdiv = coefdis/(real(nitergdiv)-2.)
151			ELSE
152			gamdi_gdiv = 0.
153			END IF
154			IF (nitergrot/=2) THEN
155			gamdi_grot = coefdis/(real(nitergrot)-2.)
156			ELSE
157			gamdi_grot = 0.
158			END IF
159			IF (niterh/=2) THEN
160			gamdi_h = coefdis/(real(niterh)-2.)
161			ELSE
162			gamdi_h = 0.
163			END IF
164	guez	7
165	guez	25	print *, 'gamdi_gdiv = ', gamdi_gdiv
166			print *, "gamdi_grot = ", gamdi_grot
167			print *, "gamdi_h = ", gamdi_h
168	guez	7
169	guez	25	WRITE (6, 990)
170	guez	7
171	guez	25	IF ( .NOT. fxyhypb) THEN
172			IF (ysinus) THEN
173			print , ' Inigeom , Y = Sinus ( Latitude ) * '
174	guez	7
175	guez	25	! utilisation de f(x, y ) avec y = sinus de la latitude ...
176	guez	7
177	guez	25	CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
178			rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
179			xprimm025, rlonp025, xprimp025)
180			ELSE
181			print , ' Inigeom , Y = Latitude , der. sinusoid . *'
182			! utilisation de f(x, y) a tangente sinusoidale , y etant la latit
183	guez	7
184	guez	25	pxo = clon*pi/180.
185			pyo = 2.clatpi/180.
186	guez	7
187	guez	25	! determination de transx ( pour le zoom ) par Newton-Raphson .
188	guez	7
189	guez	25	itmax = 10
190			eps = .1E-7
191	guez	7
192	guez	25	xo1 = 0.
193			DO iter = 1, itmax
194			x1 = xo1
195			f = x1 + alphax*sin(x1-pxo)
196			df = 1. + alphax*cos(x1-pxo)
197			x1 = x1 - f/df
198			xdm = abs(x1-xo1)
199			IF (xdm<=eps) EXIT
200			xo1 = x1
201			END DO
202	guez	7
203	guez	25	transx = xo1
204	guez	7
205	guez	25	itmay = 10
206			eps = .1E-7
207	guez	7
208	guez	25	yo1 = 0.
209			DO iter = 1, itmay
210			y1 = yo1
211			f = y1 + alphay*sin(y1-pyo)
212			df = 1. + alphay*cos(y1-pyo)
213			y1 = y1 - f/df
214			ydm = abs(y1-yo1)
215			IF (ydm<=eps) EXIT
216			yo1 = y1
217			END DO
218	guez	7
219	guez	25	transy = yo1
220	guez	7
221	guez	25	CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
222			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
223			rlonp025, xprimp025)
224			END IF
225			ELSE
226			! .... Utilisation de fxyhyper , f(x, y) a derivee tangente hyperbol.
227			print , ' Inigeom , Y = Latitude , der.tg. hyperbolique *'
228			CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
229			taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
230			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
231			rlonp025, xprimp025)
232			END IF
233	guez	7
234	guez	25	rlatu(1) = asin(1.)
235			rlatu(jjp1) = -rlatu(1)
236	guez	7
237	guez	25	! .... calcul aux poles ....
238	guez	7
239	guez	25	yprimu(1) = 0.
240			yprimu(jjp1) = 0.
241	guez	7
242	guez	25	un4rad2 = 0.25radrad
243	guez	7
244	guez	25	! calcul des aires ( aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez )
245			! - et de fext_2d , force de coriolis extensive .
246	guez	7
247	guez	25	! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y) , sont
248			! affectees 4 aires entourant P , calculees respectivement aux points
249			! ( i + 1/4, j - 1/4 ) : aireij1_2d (i, j)
250			! ( i + 1/4, j + 1/4 ) : aireij2_2d (i, j)
251			! ( i - 1/4, j + 1/4 ) : aireij3_2d (i, j)
252			! ( i - 1/4, j - 1/4 ) : aireij4_2d (i, j)
253	guez	7
254	guez	25	! ,
255			! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).
256			! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme
257			! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont
258			! affectees au
259			! point (i, j) .
260			! On definit en outre les coefficients alpha comme etant egaux a
261			! (aireij / aire_2d), c.a.d par exp.
262			! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)
263	guez	7
264	guez	25	! De meme, toute aire centree en 1 point U est egale a la somme des
265			! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
266			! le point U.
267			! Idem pour airev_2d, airez .
268	guez	7
269	guez	25	! On a , pour chaque maille : dX = dY = 1
270	guez	7
271	guez	25	! . V
272	guez	7
273	guez	25	! aireij4_2d . . aireij1_2d
274	guez	7
275	guez	25	! U . . P . U
276	guez	7
277	guez	25	! aireij3_2d . . aireij2_2d
278	guez	7
279	guez	25	! . V
280	guez	3
281	guez	25	! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
282			! aireij3_2d, aireij4_2d
283			! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations
284			! elementaires
285			! cuij et les 4 elongat. cvij qui sont calculees aux memes
286			! endroits que les aireij .
287	guez	7
288	guez	25	! ....... do 35 : boucle sur les jjm + 1 latitudes .....
289	guez	7
290	guez	25	DO j = 1, jjp1
291	guez	3
292	guez	25	IF (j==1) THEN
293	guez	3
294	guez	25	yprm = yprimu1(j)
295			rlatm = rlatu1(j)
296	guez	3
297	guez	25	coslatm = cos(rlatm)
298			radclatm = 0.5radcoslatm
299	guez	3
300	guez	25	DO i = 1, iim
301			xprp = xprimp025(i)
302			xprm = xprimm025(i)
303			aireij2_2d(i, 1) = un4rad2coslatmxprp*yprm
304			aireij3_2d(i, 1) = un4rad2coslatmxprm*yprm
305			cuij2(i, 1) = radclatm*xprp
306			cuij3(i, 1) = radclatm*xprm
307			cvij2(i, 1) = 0.5radyprm
308			cvij3(i, 1) = cvij2(i, 1)
309			END DO
310	guez	3
311	guez	25	DO i = 1, iim
312			aireij1_2d(i, 1) = 0.
313			aireij4_2d(i, 1) = 0.
314			cuij1(i, 1) = 0.
315			cuij4(i, 1) = 0.
316			cvij1(i, 1) = 0.
317			cvij4(i, 1) = 0.
318			END DO
319	guez	3
320	guez	25	END IF
321	guez	3
322	guez	25	IF (j==jjp1) THEN
323			yprp = yprimu2(j-1)
324			rlatp = rlatu2(j-1)
325	guez	3
326	guez	25	coslatp = cos(rlatp)
327			radclatp = 0.5radcoslatp
328	guez	3
329	guez	25	DO i = 1, iim
330			xprp = xprimp025(i)
331			xprm = xprimm025(i)
332			aireij1_2d(i, jjp1) = un4rad2coslatpxprp*yprp
333			aireij4_2d(i, jjp1) = un4rad2coslatpxprm*yprp
334			cuij1(i, jjp1) = radclatp*xprp
335			cuij4(i, jjp1) = radclatp*xprm
336			cvij1(i, jjp1) = 0.5radyprp
337			cvij4(i, jjp1) = cvij1(i, jjp1)
338			END DO
339	guez	3
340	guez	25	DO i = 1, iim
341			aireij2_2d(i, jjp1) = 0.
342			aireij3_2d(i, jjp1) = 0.
343			cvij2(i, jjp1) = 0.
344			cvij3(i, jjp1) = 0.
345			cuij2(i, jjp1) = 0.
346			cuij3(i, jjp1) = 0.
347			END DO
348	guez	3
349	guez	25	END IF
350	guez	3
351	guez	25	IF (j>1 .AND. j<jjp1) THEN
352	guez	3
353	guez	25	rlatp = rlatu2(j-1)
354			yprp = yprimu2(j-1)
355			rlatm = rlatu1(j)
356			yprm = yprimu1(j)
357	guez	3
358	guez	25	coslatm = cos(rlatm)
359			coslatp = cos(rlatp)
360			radclatp = 0.5radcoslatp
361			radclatm = 0.5radcoslatm
362	guez	3
363	guez	25	DO i = 1, iim
364			xprp = xprimp025(i)
365			xprm = xprimm025(i)
366	guez	7
367	guez	25	ai14 = un4rad2coslatpyprp
368			ai23 = un4rad2coslatmyprm
369			aireij1_2d(i, j) = ai14*xprp
370			aireij2_2d(i, j) = ai23*xprp
371			aireij3_2d(i, j) = ai23*xprm
372			aireij4_2d(i, j) = ai14*xprm
373			cuij1(i, j) = radclatp*xprp
374			cuij2(i, j) = radclatm*xprp
375			cuij3(i, j) = radclatm*xprm
376			cuij4(i, j) = radclatp*xprm
377			cvij1(i, j) = 0.5radyprp
378			cvij2(i, j) = 0.5radyprm
379			cvij3(i, j) = cvij2(i, j)
380			cvij4(i, j) = cvij1(i, j)
381			END DO
382	guez	7
383	guez	25	END IF
384	guez	7
385	guez	25	! ........ periodicite ............
386	guez	7
387	guez	25	cvij1(iip1, j) = cvij1(1, j)
388			cvij2(iip1, j) = cvij2(1, j)
389			cvij3(iip1, j) = cvij3(1, j)
390			cvij4(iip1, j) = cvij4(1, j)
391			cuij1(iip1, j) = cuij1(1, j)
392			cuij2(iip1, j) = cuij2(1, j)
393			cuij3(iip1, j) = cuij3(1, j)
394			cuij4(iip1, j) = cuij4(1, j)
395			aireij1_2d(iip1, j) = aireij1_2d(1, j)
396			aireij2_2d(iip1, j) = aireij2_2d(1, j)
397			aireij3_2d(iip1, j) = aireij3_2d(1, j)
398			aireij4_2d(iip1, j) = aireij4_2d(1, j)
399	guez	7
400	guez	25	END DO
401	guez	7
402	guez	25	DO j = 1, jjp1
403			DO i = 1, iim
404			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
405			+ aireij3_2d(i, j) + aireij4_2d(i, j)
406			alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)
407			alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)
408			alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)
409			alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)
410			alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
411			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
412			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
413			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
414			END DO
415	guez	7
416	guez	25	aire_2d(iip1, j) = aire_2d(1, j)
417			alpha1_2d(iip1, j) = alpha1_2d(1, j)
418			alpha2_2d(iip1, j) = alpha2_2d(1, j)
419			alpha3_2d(iip1, j) = alpha3_2d(1, j)
420			alpha4_2d(iip1, j) = alpha4_2d(1, j)
421			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
422			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
423			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
424			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
425			END DO
426	guez	7
427	guez	25	DO j = 1, jjp1
428			DO i = 1, iim
429			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
430			aireij4_2d(i+1, j) + aireij3_2d(i+1, j)
431			unsaire_2d(i, j) = 1./aire_2d(i, j)
432			unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
433			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
434			airesurg_2d(i, j) = aire_2d(i, j)/g
435			END DO
436			aireu_2d(iip1, j) = aireu_2d(1, j)
437			unsaire_2d(iip1, j) = unsaire_2d(1, j)
438			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
439			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
440			airesurg_2d(iip1, j) = airesurg_2d(1, j)
441			END DO
442	guez	7
443	guez	25	DO j = 1, jjm
444	guez	7
445	guez	25	DO i = 1, iim
446			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
447			aireij1_2d(i, j+1) + aireij4_2d(i, j+1)
448			END DO
449			DO i = 1, iim
450			airez = aireij2_2d(i, j) + aireij1_2d(i, j+1) + aireij3_2d(i+1, j) &
451			+ aireij4_2d(i+1, j+1)
452			unsairez_2d(i, j) = 1./airez
453			unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
454			fext_2d(i, j) = airezsin(rlatv(j))2.*omeg
455			END DO
456			airev_2d(iip1, j) = airev_2d(1, j)
457			unsairez_2d(iip1, j) = unsairez_2d(1, j)
458			fext_2d(iip1, j) = fext_2d(1, j)
459			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
460	guez	7
461	guez	25	END DO
462	guez	7
463	guez	25	! ..... Calcul des elongations cu_2d, cv_2d, cvu .........
464	guez	7
465	guez	25	DO j = 1, jjm
466			DO i = 1, iim
467			cv_2d(i, j) = 0.5 * &
468			(cvij2(i, j) + cvij3(i, j) + cvij1(i, j+1) + cvij4(i, j+1))
469			cvu(i, j) = 0.5*(cvij1(i, j)+cvij4(i, j)+cvij2(i, j)+cvij3(i, j))
470			cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))
471			unscv2_2d(i, j) = 1./(cv_2d(i, j)*cv_2d(i, j))
472			END DO
473			DO i = 1, iim
474			cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)
475			cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)
476			cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
477			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
478			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
479			END DO
480			cv_2d(iip1, j) = cv_2d(1, j)
481			cvu(iip1, j) = cvu(1, j)
482			unscv2_2d(iip1, j) = unscv2_2d(1, j)
483			cuv(iip1, j) = cuv(1, j)
484			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
485			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
486			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
487			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
488			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
489			END DO
490	guez	7
491	guez	25	DO j = 2, jjm
492			DO i = 1, iim
493			cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i+1, j) + cuij2(i, j) &
494			+ cuij3(i+1, j))
495			unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))
496			cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)
497			cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)
498			cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
499			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
500			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
501			END DO
502			cu_2d(iip1, j) = cu_2d(1, j)
503			unscu2_2d(iip1, j) = unscu2_2d(1, j)
504			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
505			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
506			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
507			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
508			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
509			END DO
510	guez	7
511	guez	25	! .... calcul aux poles ....
512	guez	7
513	guez	25	DO i = 1, iip1
514			cu_2d(i, 1) = 0.
515			unscu2_2d(i, 1) = 0.
516			cvu(i, 1) = 0.
517	guez	7
518	guez	25	cu_2d(i, jjp1) = 0.
519			unscu2_2d(i, jjp1) = 0.
520			cvu(i, jjp1) = 0.
521			END DO
522	guez	7
523	guez	25	DO j = 1, jjm
524			DO i = 1, iim
525			airvscu2_2d(i, j) = airev_2d(i, j)/(cuv(i, j)*cuv(i, j))
526			aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
527			END DO
528			airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
529			aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
530			END DO
531	guez	7
532	guez	25	DO j = 2, jjm
533			DO i = 1, iim
534			airuscv2_2d(i, j) = aireu_2d(i, j)/(cvu(i, j)*cvu(i, j))
535			aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
536			END DO
537			airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
538			aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
539			END DO
540	guez	7
541	guez	25	! calcul des aires aux poles :
542	guez	7
543	guez	25	apoln = sum(aire_2d(:iim, 1))
544			apols = sum(aire_2d(:iim, jjp1))
545			unsapolnga1 = 1./(apoln**(-gamdi_gdiv))
546			unsapolsga1 = 1./(apols**(-gamdi_gdiv))
547			unsapolnga2 = 1./(apoln**(-gamdi_h))
548			unsapolsga2 = 1./(apols**(-gamdi_h))
549	guez	7
550	guez	25	! changement F. Hourdin calcul conservatif pour fext_2d
551			! constang_2d contient le produit a * cos ( latitude ) * omega
552	guez	7
553	guez	25	DO i = 1, iim
554			constang_2d(i, 1) = 0.
555			END DO
556			DO j = 1, jjm - 1
557			DO i = 1, iim
558			constang_2d(i, j+1) = radomegcu_2d(i, j+1)*cos(rlatu(j+1))
559			END DO
560			END DO
561			DO i = 1, iim
562			constang_2d(i, jjp1) = 0.
563			END DO
564	guez	7
565	guez	25	! periodicite en longitude
566	guez	7
567	guez	25	DO j = 1, jjm
568			fext_2d(iip1, j) = fext_2d(1, j)
569			END DO
570			DO j = 1, jjp1
571			constang_2d(iip1, j) = constang_2d(1, j)
572			END DO
573	guez	7
574	guez	25	! fin du changement
575	guez	7
576	guez	25	print , ' Coordonnees de la grille * '
577			WRITE (6, 995)
578	guez	7
579	guez	25	print *, ' LONGITUDES aux pts. V ( degres ) '
580			WRITE (6, 995)
581			DO i = 1, iip1
582			rlonvv(i) = rlonv(i)*180./pi
583			END DO
584			WRITE (6, 400) rlonvv
585	guez	7
586	guez	25	WRITE (6, 995)
587			print *, ' LATITUDES aux pts. V ( degres ) '
588			WRITE (6, 995)
589			DO i = 1, jjm
590			rlatuu(i) = rlatv(i)*180./pi
591			END DO
592			WRITE (6, 400) (rlatuu(i), i=1, jjm)
593	guez	7
594	guez	25	DO i = 1, iip1
595			rlonvv(i) = rlonu(i)*180./pi
596			END DO
597			WRITE (6, 995)
598			print *, ' LONGITUDES aux pts. U ( degres ) '
599			WRITE (6, 995)
600			WRITE (6, 400) rlonvv
601			WRITE (6, 995)
602	guez	7
603	guez	25	print *, ' LATITUDES aux pts. U ( degres ) '
604			WRITE (6, 995)
605			DO i = 1, jjp1
606			rlatuu(i) = rlatu(i)*180./pi
607			END DO
608			WRITE (6, 400) (rlatuu(i), i=1, jjp1)
609			WRITE (6, 995)
610	guez	7
611	guez	25	400 FORMAT (1X, 8F8.2)
612	guez	7	990 FORMAT (//)
613			995 FORMAT (/)
614
615	guez	25	END SUBROUTINE inigeom
616
617			end module inigeom_m