libf/dyn3d/inigeom.f90

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_2d,..aireij4_2d .

  !  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
  ! Possibilité d'appeler une fonction "f(y)" à
  ! dérivée tangente hyperbolique à la place de la fonction à dérivée
  ! sinusoïdale.


  USE dimens_m
  USE paramet_m
  USE comconst
  USE comdissnew
  USE logic
  USE comgeom
  USE serre
  IMPLICIT NONE


  !------------------------------------------------------------------
  !   ....  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

  !   calcul des coeff. ( 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)


  PRINT *, 'Call sequence information: inigeom'
  PRINT 3
3 FORMAT ('Calcul des elongations cu_2d et cv_2d  comme sommes ', &
       'des 4 '/5X, &
       ' elong. cuij1, .. 4  , cvij1,.. 4  qui les entourent , aux '/5X, &
       ' memes endroits que les aires aireij1_2d,...j4   . '/)


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

  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  ...

        CALL fxysinus(rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,rlatu2, &
             yprimu2,rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025, &
             xprimp025)

     ELSE

        WRITE (6,*) '*** 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.
     !   ..................................................................

     WRITE (6,*) '*** 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)
        !cc       yprp             = fyprim( FLOAT(j) - 0.25 )
        !cc       rlatp            = fy    ( FLOAT(j) - 0.25 )

        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)
        !c         rlatp    = fy    ( FLOAT(j) - 0.25 )
        !c         yprp     = fyprim( FLOAT(j) - 0.25 )
        !c         rlatm    = fy    ( FLOAT(j) + 0.25 )
        !c         yprm     = fyprim( FLOAT(j) + 0.25 )

        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))

  !----------------------------------------------------------------
  !     gtitre='Coriolis version ancienne'
  !     gfichier='fext1'
  !     CALL writestd(fext_2d,iip1*jjm)

  !   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


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

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

  WRITE (6,*) '   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)
  WRITE (6,*) '   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)
  WRITE (6,*) '   LONGITUDES  aux pts.   U  ( degres )  '
  WRITE (6,995)
  WRITE (6,400) rlonvv
  WRITE (6,995)

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