--- trunk/dyn3d/tau2alpha.f	2014/09/18 13:36:51	112
+++ trunk/dyn3d/Guide/tau2alpha.f	2014/09/19 17:36:20	115
@@ -4,131 +4,44 @@
 
 contains
 
-  SUBROUTINE tau2alpha(type, factt, taumin, taumax, alpha)
+  SUBROUTINE tau2alpha(dxdy, rlat, taumin, taumax, alpha)
 
-    USE comgeom, ONLY: cu_2d, cv_2d, rlatu, rlatv
-    use conf_guide_m, only: lat_min_guide, lat_max_guide
-    USE dimens_m, ONLY: iim, jjm
-    USE nr_util, ONLY: pi
-    USE paramet_m, ONLY: iip1, jjp1
-    USE serre, ONLY: clat, clon, grossismx, grossismy
-    use writefield_m, only: writefield
-
-    INTEGER, intent(in):: type
-    REAL, intent(in):: factt, taumin, taumax
-    real, intent(out):: alpha(:, :)
+    use conf_guide_m, only: lat_min_guide, lat_max_guide, factt
+    use init_tau2alpha_m, only: dxdy_min, dxdy_max, gamma
+    USE nr_util, ONLY: pi, assert_eq
+
+    REAL, intent(in):: dxdy(:, :) ! (n_lon, n_lat)
+    REAL, intent(in):: rlat(:) ! (n_lat)
+    REAL, intent(in):: taumin, taumax
+    real, intent(out):: alpha(:, :) ! (n_lon, n_lat)
 
     ! Local:
-    REAL dxdy
-    REAL dxdy_min, dxdy_max
     REAL alphamin, alphamax, xi
-    REAL, SAVE:: gamma
-    INTEGER i, j, ilon, ilat
-    LOGICAL:: first = .TRUE.
-    REAL zdx(iip1, jjp1), zdy(iip1, jjp1)
+    INTEGER i, j, n_lon, n_lat
     REAL zlat
-    REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
 
     !------------------------------------------------------------
 
-    IF (first) THEN
-       DO j = 2, jjm
-          DO i = 2, iip1
-             zdx(i, j) = 0.5 * (cu_2d(i - 1, j) + cu_2d(i, j)) / cos(rlatu(j))
-          END DO
-          zdx(1, j) = zdx(iip1, j)
-       END DO
-       DO j = 2, jjm
-          DO i = 1, iip1
-             zdy(i, j) = 0.5 * (cv_2d(i, j - 1) + cv_2d(i, j))
-          END DO
-       END DO
-       DO i = 1, iip1
-          zdx(i, 1) = zdx(i, 2)
-          zdx(i, jjp1) = zdx(i, jjm)
-          zdy(i, 1) = zdy(i, 2)
-          zdy(i, jjp1) = zdy(i, jjm)
-       END DO
-
-       DO j = 1, jjp1
-          DO i = 1, iip1
-             dxdys(i, j) = sqrt(zdx(i, j)**2 + zdy(i, j)**2)
-          END DO
-       END DO
-       CALL writefield("dxdys", dxdys)
+    PRINT *, 'Call sequence information: tau2alpha'
 
-       if (type == 2) then
-          DO j = 1, jjp1
-             DO i = 1, iim
-                dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
-             END DO
-             dxdyu(iip1, j) = dxdyu(1, j)
-          END DO
-       elseif (type == 3) then
-          DO j = 1, jjm
-             DO i = 1, iip1
-                dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
-             END DO
-          END DO
-       end if
-
-       ! coordonnees du centre du zoom
-       CALL coordij(clon, clat, ilon, ilat)
-       ! aire de la maille au centre du zoom
-       dxdy_min = dxdys(ilon, ilat)
-
-       ! dxdy maximal de la maille :
-       dxdy_max = 0.
-       DO j = 1, jjp1
-          DO i = 1, iip1
-             dxdy_max = max(dxdy_max, dxdys(i, j))
-          END DO
-       END DO
-
-       IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
-          PRINT *, 'Attention : modèle peu zoomé.'
-          PRINT *, 'On prend une constante de guidage constante.'
-       ELSE
-          gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
-          PRINT *, 'gamma=', gamma
-          IF (gamma<1.E-5) THEN
-             PRINT *, 'gamma =', gamma, '<1e-5'
-             STOP
-          END IF
-          PRINT *, 'gamma=', gamma
-          gamma = log(0.5) / log(gamma)
-       END IF
-       first = .false.
-    END IF
+    n_lon = assert_eq(size(alpha, 1), size(dxdy, 1), "tau2alpha n_lon")
+    n_lat = assert_eq(size(alpha, 2), size(dxdy, 2), size(rlat), &
+         "tau2alpha n_lat")
 
     alphamin = factt / taumax
     alphamax = factt / taumin
 
-    DO j = 1, size(alpha, 2)
-       DO i = 1, size(alpha, 1)
-          IF (type==1) THEN
-             dxdy = dxdys(i, j)
-             zlat = rlatu(j) * 180. / pi
-          ELSE IF (type==2) THEN
-             dxdy = dxdyu(i, j)
-             zlat = rlatu(j) * 180. / pi
-          ELSE IF (type==3) THEN
-             dxdy = dxdyv(i, j)
-             zlat = rlatv(j) * 180. / pi
-          END IF
-          IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
-             ! grille regulière
-             alpha(i, j) = alphamin
-          ELSE
-             xi = ((dxdy_max - dxdy) / (dxdy_max - dxdy_min))**gamma
-             xi = min(xi, 1.)
-             IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN
-                alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
-             ELSE
-                alpha(i, j) = 0.
-             END IF
-          END IF
-       END DO
+    DO j = 1, n_lat
+       zlat = rlat(j) * 180. / pi
+       IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN
+          DO i = 1, n_lon
+             xi = min(((dxdy_max - dxdy(i, j)) &
+                  / (dxdy_max - dxdy_min))**gamma, 1.)
+             alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
+          END DO
+       ELSE
+          alpha(:, j) = 0.
+       END IF
     END DO
 
   END SUBROUTINE tau2alpha