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	guez	7	SUBROUTINE inigeom
2	guez	3
3	guez	7	! Auteur : P. Le Van
4	guez	3
5	guez	7	! ............ Version du 01/04/2001 ...................
6	guez	3
7	guez	7	! Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en-
8			! endroits que les aires aireij1_2d,..aireij4_2d .
9	guez	3
10	guez	7	! 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	guez	3
15
16	guez	7	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	guez	3
25
26	guez	7	!------------------------------------------------------------------
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	guez	3	itmax = 10
213	guez	7	eps = .1E-7
214
215	guez	3	xo1 = 0.
216	guez	7	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	guez	3	transx = xo1
227
228			itmay = 10
229	guez	7	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	guez	3	transy = yo1
243
244	guez	7	CALL fxy(rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,rlatu2,yprimu2, &
245			rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
246	guez	3
247	guez	7	END IF
248	guez	3
249	guez	7	ELSE
250	guez	3
251	guez	7	! .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
252			! ..................................................................
253	guez	3
254	guez	7	WRITE (6,) ' Inigeom , Y = Latitude , der.tg. hyperbolique *'
255	guez	3
256	guez	7	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	guez	3
260
261	guez	7	END IF
262	guez	3
263	guez	7	! -------------------------------------------------------------------
264	guez	3
265
266	guez	7	rlatu(1) = asin(1.)
267			rlatu(jjp1) = -rlatu(1)
268	guez	3
269
270	guez	7	! .... calcul aux poles ....
271	guez	3
272	guez	7	yprimu(1) = 0.
273			yprimu(jjp1) = 0.
274	guez	3
275	guez	7
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