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