libf/dyn3d/inigeom.f90

module inigeom_m

  IMPLICIT NONE

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van

    ! 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 coefficients cu_2d et 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(contravariant) = u(naturel) / cu_2d
    ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / 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 calculés aux points concernés. cv, bien que
    ! dépendant de j uniquement, sera ici indicé aussi en i pour un
    ! adressage plus facile en ij.

    ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
    ! points u et v. yprimu et yprimv sont respectivement les valeurs
    ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
    ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
    ! respectivement les valeurs de cv_2d aux points u et v.

    ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
    ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".

    USE comconst, ONLY : g, omeg, rad
    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 comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
    use conf_gcm_m, ONLY : fxyhypb, ysinus
    USE dimens_m, ONLY : iim, jjm
    use fxy_m, only: fxy
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1
    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, unit
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, un4rad2
    REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
    REAL coslatm, coslatp, radclatm, radclatp
    REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
    REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
    REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
    real yprimv(jjm), yprimu(jjp1)
    REAL gamdi_gdiv, gamdi_grot, gamdi_h
    REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
    real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
         aireij4_2d ! in m2
    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

    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) à dérivée tangente hyperbolique
       print *, 'Inigeom, Y = Latitude, dérivée tangente 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) = pi / 2.
    rlatu(jjp1) = -rlatu(1)

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

    ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
    ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
    ! chaque aire_2d(i, j), ainsi que les quatre élongations
    ! élémentaires cuij et les quatre élongations cvij qui sont
    ! calculées aux mêmes endroits que les aireij.

    coslatm = cos(rlatu1(1))
    radclatm = 0.5 * rad * coslatm

    aireij1_2d(:iim, 1) = 0.
    aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
    aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
    aireij4_2d(:iim, 1) = 0.

    cuij1(:iim, 1) = 0.
    cuij2(:iim, 1) = radclatm * xprimp025(:iim)
    cuij3(:iim, 1) = radclatm * xprimm025(:iim)
    cuij4(:iim, 1) = 0.

    cvij1(:iim, 1) = 0.
    cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
    cvij3(:iim, 1) = cvij2(:iim, 1)
    cvij4(:iim, 1) = 0.

    do j = 2, jjm
       coslatm = cos(rlatu1(j))
       coslatp = cos(rlatu2(j-1))
       radclatp = 0.5 * rad * coslatp
       radclatm = 0.5 * rad * coslatm
       ai14 = un4rad2 * coslatp * yprimu2(j-1)
       ai23 = un4rad2 * coslatm * yprimu1(j)

       aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
       aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
       aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
       aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
       cuij1(:iim, j) = radclatp * xprimp025(:iim)
       cuij2(:iim, j) = radclatm * xprimp025(:iim)
       cuij3(:iim, j) = radclatm * xprimm025(:iim)
       cuij4(:iim, j) = radclatp * xprimm025(:iim)
       cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
       cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
       cvij3(:iim, j) = cvij2(:iim, j)
       cvij4(:iim, j) = cvij1(:iim, j)
    end do

    coslatp = cos(rlatu2(jjm))
    radclatp = 0.5 * rad * coslatp

    aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
    aireij2_2d(:iim, jjp1) = 0.
    aireij3_2d(:iim, jjp1) = 0.
    aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)

    cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
    cuij2(:iim, jjp1) = 0.
    cuij3(:iim, jjp1) = 0.
    cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)

    cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
    cvij2(:iim, jjp1) = 0.
    cvij3(:iim, jjp1) = 0.
    cvij4(:iim, jjp1) = cvij1(:iim, jjp1)

    ! Périodicité :

    cvij1(iip1, :) = cvij1(1, :)
    cvij2(iip1, :) = cvij2(1, :)
    cvij3(iip1, :) = cvij3(1, :)
    cvij4(iip1, :) = cvij4(1, :)

    cuij1(iip1, :) = cuij1(1, :)
    cuij2(iip1, :) = cuij2(1, :)
    cuij3(iip1, :) = cuij3(1, :)
    cuij4(iip1, :) = cuij4(1, :)

    aireij1_2d(iip1, :) = aireij1_2d(1, :)
    aireij2_2d(iip1, :) = aireij2_2d(1, :)
    aireij3_2d(iip1, :) = aireij3_2d(1, :)
    aireij4_2d(iip1, :) = aireij4_2d(1, :)

    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 élongations 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)**2
       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)**2
          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 pôles 

    cu_2d(:, 1) = 0.
    unscu2_2d(:, 1) = 0.
    cvu(:, 1) = 0.

    cu_2d(:, jjp1) = 0.
    unscu2_2d(:, jjp1) = 0.
    cvu(:, jjp1) = 0.

    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 pôles :

    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

    ! Périodicité 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

    call new_unit(unit)
    open(unit, file="longitude_latitude.txt", status="replace", action="write")
    write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
    write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
    write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
    write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
    close(unit)

  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	guez	7
11	guez	25	! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
12			! endroits que les aires aireij1_2d, ..., aireij4_2d.
13	guez	7
14	guez	57	! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à
15			! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,
16			! cv_2d, 1. / cu_2d2, 1. / cv_2d2. Les coefficients cu_2d et cv_2d
17			! permettent de passer des vitesses naturelles aux vitesses
18			! covariantes et contravariantes, ou vice-versa.
19	guez	7
20	guez	57	! On a :
21			! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
22			! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
23	guez	7
24	guez	57	! On en tire :
25	guez	38	! u(covariant) = cu_2d * cu_2d * u(contravariant)
26			! v(covariant) = cv_2d * cv_2d * v(contravariant)
27	guez	7
28	guez	57	! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
29			! et - jm / 2 <= Y <= jm / 2
30	guez	7
31	guez	38	! x est la longitude du point en radians.
32			! y est la latitude du point en radians.
33			!
34	guez	57	! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
35			! cv(j) = rad * dy / dY
36	guez	38	! aire_2d(i, j) = cu_2d(i, j) * cv(j)
37			!
38	guez	57	! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
39			! dépendant de j uniquement, sera ici indicé aussi en i pour un
40			! adressage plus facile en ij.
41	guez	7
42	guez	57	! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
43			! points u et v. yprimu et yprimv sont respectivement les valeurs
44			! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
45			! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
46			! respectivement les valeurs de cv_2d aux points u et v.
47	guez	7
48	guez	38	! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
49	guez	57	! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
50	guez	7
51	guez	39	USE comconst, ONLY : g, omeg, rad
52	guez	25	USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
53			alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
54			alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
55			apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &
56			cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &
57			cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &
58			rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
59			unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
60			unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
61	guez	39	USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
62	guez	57	use conf_gcm_m, ONLY : fxyhypb, ysinus
63	guez	39	USE dimens_m, ONLY : iim, jjm
64	guez	57	use fxy_m, only: fxy
65			use jumble, only: new_unit
66	guez	39	use nr_util, only: pi
67			USE paramet_m, ONLY : iip1, jjp1
68	guez	25	USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
69			grossismy, pxo, pyo, taux, tauy, transx, transy
70	guez	7
71	guez	25	! Variables locales
72	guez	7
73	guez	57	INTEGER i, j, itmax, itmay, iter, unit
74	guez	25	REAL cvu(iip1, jjp1), cuv(iip1, jjm)
75	guez	57	REAL ai14, ai23, airez, un4rad2
76	guez	25	REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
77			REAL coslatm, coslatp, radclatm, radclatp
78	guez	57	REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
79	guez	60	REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
80	guez	57	REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
81			real yprimv(jjm), yprimu(jjp1)
82	guez	25	REAL gamdi_gdiv, gamdi_grot, gamdi_h
83			REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
84	guez	57	real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
85			aireij4_2d ! in m2
86	guez	25	real airuscv2_2d(iim + 1, jjm)
87	guez	38	real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
88	guez	25	real aivscu2gam_2d(iim + 1, jjm)
89	guez	7
90	guez	25	!------------------------------------------------------------------
91	guez	7
92	guez	25	PRINT *, 'Call sequence information: inigeom'
93	guez	7
94	guez	25	IF (nitergdiv/=2) THEN
95	guez	57	gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
96	guez	25	ELSE
97			gamdi_gdiv = 0.
98			END IF
99			IF (nitergrot/=2) THEN
100	guez	57	gamdi_grot = coefdis / (real(nitergrot)-2.)
101	guez	25	ELSE
102			gamdi_grot = 0.
103			END IF
104			IF (niterh/=2) THEN
105	guez	57	gamdi_h = coefdis / (real(niterh)-2.)
106	guez	25	ELSE
107			gamdi_h = 0.
108			END IF
109	guez	7
110	guez	25	print *, 'gamdi_gdiv = ', gamdi_gdiv
111			print *, "gamdi_grot = ", gamdi_grot
112			print *, "gamdi_h = ", gamdi_h
113	guez	7
114	guez	38	IF (.NOT. fxyhypb) THEN
115	guez	25	IF (ysinus) THEN
116	guez	38	print *, ' Inigeom, Y = Sinus (Latitude) '
117			! utilisation de f(x, y) avec y = sinus de la latitude
118	guez	25	CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
119			rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
120			xprimm025, rlonp025, xprimp025)
121			ELSE
122	guez	38	print *, 'Inigeom, Y = Latitude, der. sinusoid .'
123			! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
124	guez	7
125	guez	57	pxo = clon * pi / 180.
126			pyo = 2. * clat * pi / 180.
127	guez	7
128	guez	38	! determination de transx (pour le zoom) par Newton-Raphson
129	guez	7
130	guez	25	itmax = 10
131			eps = .1E-7
132	guez	7
133	guez	25	xo1 = 0.
134			DO iter = 1, itmax
135			x1 = xo1
136	guez	57	f = x1 + alphax * sin(x1-pxo)
137			df = 1. + alphax * cos(x1-pxo)
138			x1 = x1 - f / df
139	guez	25	xdm = abs(x1-xo1)
140			IF (xdm<=eps) EXIT
141			xo1 = x1
142			END DO
143	guez	7
144	guez	25	transx = xo1
145	guez	7
146	guez	25	itmay = 10
147			eps = .1E-7
148	guez	7
149	guez	25	yo1 = 0.
150			DO iter = 1, itmay
151			y1 = yo1
152	guez	57	f = y1 + alphay * sin(y1-pyo)
153			df = 1. + alphay * cos(y1-pyo)
154			y1 = y1 - f / df
155	guez	25	ydm = abs(y1-yo1)
156			IF (ydm<=eps) EXIT
157			yo1 = y1
158			END DO
159	guez	7
160	guez	25	transy = yo1
161	guez	7
162	guez	25	CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
163			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
164			rlonp025, xprimp025)
165			END IF
166			ELSE
167	guez	38	! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
168			print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
169	guez	25	CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
170			taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
171			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
172			rlonp025, xprimp025)
173			END IF
174	guez	7
175	guez	57	rlatu(1) = pi / 2.
176	guez	25	rlatu(jjp1) = -rlatu(1)
177	guez	7
178	guez	57	! Calcul aux pôles
179	guez	7
180	guez	25	yprimu(1) = 0.
181			yprimu(jjp1) = 0.
182	guez	7
183	guez	57	un4rad2 = 0.25 * rad * rad
184	guez	7
185	guez	57	! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
186			! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
187			! chaque aire_2d(i, j), ainsi que les quatre élongations
188			! élémentaires cuij et les quatre élongations cvij qui sont
189			! calculées aux mêmes endroits que les aireij.
190	guez	7
191	guez	57	coslatm = cos(rlatu1(1))
192			radclatm = 0.5 * rad * coslatm
193	guez	7
194	guez	57	aireij1_2d(:iim, 1) = 0.
195			aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
196			aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
197			aireij4_2d(:iim, 1) = 0.
198	guez	7
199	guez	57	cuij1(:iim, 1) = 0.
200			cuij2(:iim, 1) = radclatm * xprimp025(:iim)
201			cuij3(:iim, 1) = radclatm * xprimm025(:iim)
202			cuij4(:iim, 1) = 0.
203	guez	7
204	guez	57	cvij1(:iim, 1) = 0.
205			cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
206			cvij3(:iim, 1) = cvij2(:iim, 1)
207			cvij4(:iim, 1) = 0.
208	guez	7
209	guez	57	do j = 2, jjm
210			coslatm = cos(rlatu1(j))
211			coslatp = cos(rlatu2(j-1))
212			radclatp = 0.5 * rad * coslatp
213			radclatm = 0.5 * rad * coslatm
214			ai14 = un4rad2 * coslatp * yprimu2(j-1)
215			ai23 = un4rad2 * coslatm * yprimu1(j)
216	guez	7
217	guez	57	aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
218			aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
219			aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
220			aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
221			cuij1(:iim, j) = radclatp * xprimp025(:iim)
222			cuij2(:iim, j) = radclatm * xprimp025(:iim)
223			cuij3(:iim, j) = radclatm * xprimm025(:iim)
224			cuij4(:iim, j) = radclatp * xprimm025(:iim)
225			cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
226			cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
227			cvij3(:iim, j) = cvij2(:iim, j)
228			cvij4(:iim, j) = cvij1(:iim, j)
229			end do
230	guez	7
231	guez	57	coslatp = cos(rlatu2(jjm))
232			radclatp = 0.5 * rad * coslatp
233	guez	7
234	guez	57	aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
235			aireij2_2d(:iim, jjp1) = 0.
236			aireij3_2d(:iim, jjp1) = 0.
237			aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
238	guez	7
239	guez	57	cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
240			cuij2(:iim, jjp1) = 0.
241			cuij3(:iim, jjp1) = 0.
242			cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
243	guez	3
244	guez	57	cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
245			cvij2(:iim, jjp1) = 0.
246			cvij3(:iim, jjp1) = 0.
247			cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
248	guez	7
249	guez	57	! Périodicité :
250	guez	61
251	guez	57	cvij1(iip1, :) = cvij1(1, :)
252			cvij2(iip1, :) = cvij2(1, :)
253			cvij3(iip1, :) = cvij3(1, :)
254			cvij4(iip1, :) = cvij4(1, :)
255	guez	7
256	guez	57	cuij1(iip1, :) = cuij1(1, :)
257			cuij2(iip1, :) = cuij2(1, :)
258			cuij3(iip1, :) = cuij3(1, :)
259			cuij4(iip1, :) = cuij4(1, :)
260	guez	3
261	guez	57	aireij1_2d(iip1, :) = aireij1_2d(1, :)
262			aireij2_2d(iip1, :) = aireij2_2d(1, :)
263			aireij3_2d(iip1, :) = aireij3_2d(1, :)
264			aireij4_2d(iip1, :) = aireij4_2d(1, :)
265	guez	3
266	guez	25	DO j = 1, jjp1
267			DO i = 1, iim
268			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
269			+ aireij3_2d(i, j) + aireij4_2d(i, j)
270	guez	57	alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
271			alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
272			alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
273			alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
274	guez	25	alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
275			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
276			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
277			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
278			END DO
279	guez	7
280	guez	25	aire_2d(iip1, j) = aire_2d(1, j)
281			alpha1_2d(iip1, j) = alpha1_2d(1, j)
282			alpha2_2d(iip1, j) = alpha2_2d(1, j)
283			alpha3_2d(iip1, j) = alpha3_2d(1, j)
284			alpha4_2d(iip1, j) = alpha4_2d(1, j)
285			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
286			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
287			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
288			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
289			END DO
290	guez	7
291	guez	25	DO j = 1, jjp1
292			DO i = 1, iim
293			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
294	guez	57	aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
295			unsaire_2d(i, j) = 1. / aire_2d(i, j)
296	guez	25	unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
297			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
298	guez	57	airesurg_2d(i, j) = aire_2d(i, j) / g
299	guez	25	END DO
300			aireu_2d(iip1, j) = aireu_2d(1, j)
301			unsaire_2d(iip1, j) = unsaire_2d(1, j)
302			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
303			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
304			airesurg_2d(iip1, j) = airesurg_2d(1, j)
305			END DO
306	guez	7
307	guez	25	DO j = 1, jjm
308			DO i = 1, iim
309			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
310	guez	57	aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
311	guez	25	END DO
312			DO i = 1, iim
313	guez	57	airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
314			+ aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
315			unsairez_2d(i, j) = 1. / airez
316	guez	25	unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
317	guez	57	fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
318	guez	25	END DO
319			airev_2d(iip1, j) = airev_2d(1, j)
320			unsairez_2d(iip1, j) = unsairez_2d(1, j)
321			fext_2d(iip1, j) = fext_2d(1, j)
322			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
323			END DO
324	guez	7
325	guez	57	! Calcul des élongations cu_2d, cv_2d, cvu
326	guez	7
327	guez	25	DO j = 1, jjm
328			DO i = 1, iim
329			cv_2d(i, j) = 0.5 * &
330	guez	57	(cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
331			cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
332			+ cvij3(i, j))
333			cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
334			+ cuij4(i, j + 1))
335	guez	60	unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
336	guez	25	END DO
337			DO i = 1, iim
338	guez	57	cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
339			cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
340	guez	25	cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
341			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
342			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
343			END DO
344			cv_2d(iip1, j) = cv_2d(1, j)
345			cvu(iip1, j) = cvu(1, j)
346			unscv2_2d(iip1, j) = unscv2_2d(1, j)
347			cuv(iip1, j) = cuv(1, j)
348			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
349			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
350			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
351			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
352			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
353			END DO
354	guez	7
355	guez	25	DO j = 2, jjm
356			DO i = 1, iim
357	guez	57	cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
358			+ cuij3(i + 1, j))
359	guez	60	unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
360	guez	57	cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
361			cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
362	guez	25	cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
363			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
364			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
365			END DO
366			cu_2d(iip1, j) = cu_2d(1, j)
367			unscu2_2d(iip1, j) = unscu2_2d(1, j)
368			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
369			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
370			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
371			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
372			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
373			END DO
374	guez	7
375	guez	57	! Calcul aux pôles
376	guez	7
377	guez	61	cu_2d(:, 1) = 0.
378			unscu2_2d(:, 1) = 0.
379			cvu(:, 1) = 0.
380	guez	7
381	guez	61	cu_2d(:, jjp1) = 0.
382			unscu2_2d(:, jjp1) = 0.
383			cvu(:, jjp1) = 0.
384	guez	7
385	guez	25	DO j = 1, jjm
386			DO i = 1, iim
387	guez	57	airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
388	guez	25	aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
389			END DO
390			airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
391			aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
392			END DO
393	guez	7
394	guez	25	DO j = 2, jjm
395			DO i = 1, iim
396	guez	57	airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
397	guez	25	aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
398			END DO
399			airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
400			aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
401			END DO
402	guez	7
403	guez	57	! Calcul des aires aux pôles :
404	guez	7
405	guez	25	apoln = sum(aire_2d(:iim, 1))
406			apols = sum(aire_2d(:iim, jjp1))
407	guez	57	unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
408			unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
409			unsapolnga2 = 1. / (apoln**(-gamdi_h))
410			unsapolsga2 = 1. / (apols**(-gamdi_h))
411	guez	7
412	guez	57	! Changement F. Hourdin calcul conservatif pour fext_2d
413	guez	38	! constang_2d contient le produit a * cos (latitude) * omega
414	guez	7
415	guez	25	DO i = 1, iim
416			constang_2d(i, 1) = 0.
417			END DO
418			DO j = 1, jjm - 1
419			DO i = 1, iim
420	guez	57	constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
421			* cos(rlatu(j + 1))
422	guez	25	END DO
423			END DO
424			DO i = 1, iim
425			constang_2d(i, jjp1) = 0.
426			END DO
427	guez	7
428	guez	57	! Périodicité en longitude
429	guez	7
430	guez	25	DO j = 1, jjm
431			fext_2d(iip1, j) = fext_2d(1, j)
432			END DO
433			DO j = 1, jjp1
434			constang_2d(iip1, j) = constang_2d(1, j)
435			END DO
436	guez	7
437	guez	57	call new_unit(unit)
438			open(unit, file="longitude_latitude.txt", status="replace", action="write")
439			write(unit, fmt=) '"longitudes at V points (degrees)"', rlonv 180. / pi
440			write(unit, fmt=) '"latitudes at V points (degrees)"', rlatv 180. / pi
441			write(unit, fmt=) '"longitudes at U points (degrees)"', rlonu 180. / pi
442			write(unit, fmt=) '"latitudes at U points (degrees)"', rlatu 180. / pi
443			close(unit)
444	guez	7
445	guez	25	END SUBROUTINE inigeom
446
447			end module inigeom_m