--- trunk/libf/dyn3d/fxy.f	2008/03/04 14:00:42	6
+++ trunk/libf/dyn3d/fxy.f90	2008/03/31 12:24:17	7
@@ -1,69 +1,128 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/fxy.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $
-!
-      SUBROUTINE fxy (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
-     ,                    rlatu2,yprimu2,
-     , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
-
-      use dimens_m
-      use paramet_m
-      use comconst
-      use serre
-      IMPLICIT NONE
-
-c     Auteur  :  P. Le Van
-c
-c     Calcul  des longitudes et des latitudes  pour une fonction f(x,y)
-c           a tangente sinusoidale et eventuellement avec zoom  .
-c
-c
-
-       INTEGER i,j
-
-       REAL rlatu(jjp1), yprimu(jjp1),rlatv(jjm), yprimv(jjm),
-     , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
-       REAL rlonu(iip1),xprimu(iip1),rlonv(iip1),xprimv(iip1),
-     , rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),xprimp025(iip1)
-
-       include "fxy_new.h"
-
-
-c    ......  calcul  des  latitudes  et de y'   .....
-c
-       DO j = 1, jjm + 1 
-          rlatu(j) = fy    ( FLOAT( j )        )
-         yprimu(j) = fyprim( FLOAT( j )        )
-       ENDDO
-
-
-       DO j = 1, jjm
-
-         rlatv(j)  = fy    ( FLOAT( j ) + 0.5  )
-         rlatu1(j) = fy    ( FLOAT( j ) + 0.25 ) 
-         rlatu2(j) = fy    ( FLOAT( j ) + 0.75 ) 
-
-        yprimv(j)  = fyprim( FLOAT( j ) + 0.5  ) 
-        yprimu1(j) = fyprim( FLOAT( j ) + 0.25 )
-        yprimu2(j) = fyprim( FLOAT( j ) + 0.75 )
-
-       ENDDO
-
-c
-c     .....  calcul   des  longitudes et de  x'   .....
-c
-       DO i = 1, iim + 1
-           rlonv(i)     = fx    (   FLOAT( i )          )
-           rlonu(i)     = fx    (   FLOAT( i ) + 0.5    )
-        rlonm025(i)     = fx    (   FLOAT( i ) - 0.25  )
-        rlonp025(i)     = fx    (   FLOAT( i ) + 0.25  )
-
-         xprimv  (i)    = fxprim (  FLOAT( i )          )
-         xprimu  (i)    = fxprim (  FLOAT( i ) + 0.5    )
-        xprimm025(i)    = fxprim (  FLOAT( i ) - 0.25   )
-        xprimp025(i)    = fxprim (  FLOAT( i ) + 0.25   )
-       ENDDO
-
-c
-       RETURN
-       END
+SUBROUTINE fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
+     yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, &
+     xprimp025)
 
+  ! From dyn3d/fxy.F, v 1.1.1.1 2004/05/19 12:53:06
+  ! Auteur : P. Le Van
+  ! Calcul des longitudes et des latitudes pour une fonction f(x, y)
+  ! à tangente sinusoïdale et éventuellement avec zoom.
+
+  USE dimens_m, ONLY : iim, jjm
+
+  IMPLICIT NONE
+
+  REAL, INTENT (OUT) :: rlatu(jjm + 1), yprimu(jjm + 1), rlatv(jjm)
+  REAL, INTENT (OUT) :: yprimv(jjm)
+  REAL, INTENT (OUT) :: rlatu1(jjm)
+  REAL, INTENT (OUT) :: yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
+  REAL, INTENT (OUT) :: rlonu(iim + 1), xprimu(iim + 1), rlonv(iim + 1)
+  REAL, INTENT (OUT) :: xprimv(iim + 1)
+  REAL, INTENT (OUT) :: rlonm025(iim + 1), xprimm025(iim + 1)
+  REAL, INTENT (OUT) :: rlonp025(iim + 1)
+  REAL, INTENT (OUT) :: xprimp025(iim + 1)
+
+  ! Variables local to the procedure:
+
+  INTEGER i, j
+
+  !------------------------------------------------------------
+
+  ! Calcul des latitudes et de y'
+
+  DO j = 1, jjm + 1
+     rlatu(j) = fy(real(j))
+     yprimu(j) = fyprim(real(j))
+  END DO
+
+  DO j = 1, jjm
+     rlatv(j) = fy(real(j) + 0.5)
+     rlatu1(j) = fy(real(j) + 0.25)
+     rlatu2(j) = fy(real(j) + 0.75)
+
+     yprimv(j) = fyprim(real(j) + 0.5)
+     yprimu1(j) = fyprim(real(j) + 0.25)
+     yprimu2(j) = fyprim(real(j) + 0.75)
+  END DO
+
+  ! Calcul des longitudes et de x'
+
+  DO i = 1, iim + 1
+     rlonv(i) = fx(real(i))
+     rlonu(i) = fx(real(i) + 0.5)
+     rlonm025(i) = fx(real(i) - 0.25)
+     rlonp025(i) = fx(real(i) + 0.25)
+
+     xprimv(i) = fxprim(real(i))
+     xprimu(i) = fxprim(real(i) + 0.5)
+     xprimm025(i) = fxprim(real(i) - 0.25)
+     xprimp025(i) = fxprim(real(i) + 0.25)
+  END DO
+
+CONTAINS
+
+  ! From grid/fxy_new.h, v 1.1.1.1 2004/05/19 12:53:05
+
+  REAL FUNCTION ripx(ri)
+    ! stretching in x
+    USE comconst, ONLY : pi
+    REAL, INTENT (IN) :: ri
+
+    ripx = (ri - 1.) * 2 * pi / REAL(iim)
+  end function ripx
+
+  !******************************************************
+
+  REAL FUNCTION fx(ri)
+    ! stretching in x
+    USE comconst, ONLY : pi
+    USE serre, ONLY : alphax, pxo, transx
+    REAL, INTENT (IN) :: ri
+
+    fx = ripx(ri) + transx + alphax * SIN(ripx(ri) + transx - pxo) - pi
+  end function fx
+
+  !******************************************************
+
+  REAL FUNCTION fxprim(ri)
+    ! stretching in x
+    USE comconst, ONLY : pi
+    USE serre, ONLY : alphax, pxo, transx
+    REAL, INTENT (IN) :: ri
+
+    fxprim = 2 * pi / REAL(iim) * (1. + alphax * COS(ripx(ri) + transx - pxo))
+  end function fxprim
+
+  !******************************************************
+
+  REAL FUNCTION bigy(rj)
+    ! stretching in y
+    USE comconst, ONLY : pi
+    REAL, INTENT (IN) :: rj
+
+    bigy = 2 * (REAL(jjm + 1) - rj) * pi / jjm
+  end function bigy
+
+  !******************************************************
+
+  REAL FUNCTION fy(rj)
+    ! stretching in y
+    USE comconst, ONLY : pi
+    USE serre, ONLY : alphay, pyo, transy
+    REAL, INTENT (IN) :: rj
+
+    fy = (bigy(rj) + transy + alphay * SIN(bigy(rj) + transy - pyo)) / 2 &
+         - pi / 2
+  end function fy
+
+  !******************************************************
+
+  REAL FUNCTION fyprim(rj)
+    ! stretching in y
+    USE comconst, ONLY : pi
+    USE serre, ONLY : alphay, pyo, transy
+    REAL, INTENT (IN) :: rj
+
+    fyprim = (pi / jjm) * (1. + alphay * COS(bigy(rj) + transy - pyo))
+  end function fyprim
+
+END SUBROUTINE fxy