--- trunk/libf/dyn3d/inidissip.f90	2011/10/07 13:11:58	53
+++ trunk/libf/dyn3d/inidissip.f90	2011/12/06 15:07:04	54
@@ -16,61 +16,46 @@
   SUBROUTINE inidissip
 
     ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
-    ! Initialisation de la dissipation horizontale 
+
+    ! Initialisation de la dissipation horizontale. Calcul des valeurs
+    ! propres des opérateurs par méthode itérative.
 
     USE comconst, ONLY : dtvr
     use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
          tetagrot, tetatemp
     USE comvert, ONLY : preff, presnivs
     USE conf_gcm_m, ONLY : iperiod
-    USE dimens_m, ONLY : jjm, llm
-    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
+    USE dimens_m, ONLY : iim, jjm, llm
+    USE paramet_m, ONLY : jjp1
     use jumble, only: new_unit
     use filtreg_m, only: filtreg
+    use gradiv2_m, only: gradiv2
 
     ! Variables local to the procedure:
     REAL zvert(llm), max_zvert
-    REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
-    REAL zhmin, zhmax
+    REAL, dimension(iim + 1, jjm + 1):: zh, zu
+    real zv(iim + 1, jjm), deltap(iim + 1, jjm + 1, llm)
     REAL zllm
-    INTEGER l, ij, idum, ii, unit
+    INTEGER l, seed_size, ii, unit
     REAL tetamin ! in s
-    REAL ran1
 
     !-----------------------------------------------------------------------
 
     PRINT *, 'Call sequence information: inidissip'
+    call random_seed(size=seed_size)
+    call random_seed(put=(/(0, ii = 1, seed_size)/))
 
-    ! Calcul des valeurs propres des opérateurs par méthode itérative :
-
-    crot = -1.
-    cdivu = -1.
-    cdivh = -1.
-
-    ! Calcul de la valeur propre de divgrad :
-
+    PRINT *, 'Calcul des valeurs propres de divgrad'
     deltap = 1.
-    idum = -1
-    zh(1) = ran1(idum) - 0.5
-    idum = 0
-    DO ij = 2, ip1jmp1
-       zh(ij) = ran1(idum) - 0.5
-    END DO
-
+    call random_number(zh)
+    zh = zh - 0.5
     CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
 
-    CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
-    IF (zhmin >= zhmax) THEN
-       PRINT *, 'zhmin zhmax', zhmin, zhmax
-       print *, 'Problème générateur aléatoire dans inidissip'
-       STOP 1
-    END IF
-
     DO l = 1, 50
        IF (lstardis) THEN
-          CALL divgrad2(1, zh, deltap, niterh, zh)
+          CALL divgrad2(1, zh, deltap, niterh, zh, -1.)
        ELSE
-          CALL divgrad(1, zh, niterh, zh)
+          CALL divgrad(1, zh, niterh, zh, -1.)
        END IF
 
        zllm = abs(maxval(zh))
@@ -82,59 +67,61 @@
     ELSE
        cdivh = zllm**(- 1. / niterh)
     END IF
+    PRINT *, 'cdivh = ', cdivh
 
-    ! Calcul des valeurs propres de gradiv (ii = 1) et nxgrarot (ii = 2) 
-
-    PRINT *, 'Calcul des valeurs propres'
+    PRINT *, 'Calcul des valeurs propres de gradiv'
+    call random_number(zu)
+    zu = zu - 0.5
+    CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
+    call random_number(zv)
+    zv = zv - 0.5
+    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
 
-    DO ii = 1, 2
-       DO ij = 1, ip1jmp1
-          zu(ij) = ran1(idum) - 0.5
-       END DO
-       CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
-       DO ij = 1, ip1jm
-          zv(ij) = ran1(idum) - 0.5
-       END DO
-       CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
-
-       DO l = 1, 50
-          IF (ii==1) THEN
-             IF (lstardis) THEN
-                CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
-             ELSE
-                CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
-             END IF
-          ELSE
-             IF (lstardis) THEN
-                CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
-             ELSE
-                CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
-             END IF
-          END IF
-
-          zllm = max(abs(maxval(zu)), abs(maxval(zv)))
-          zu = zu / zllm
-          zv = zv / zllm
-       end DO
-
-       IF (ii==1) THEN
-          IF (lstardis) THEN
-             cdivu = 1. / zllm
-          ELSE
-             cdivu = zllm**(- 1. / nitergdiv)
-          END IF
+    DO l = 1, 50
+       IF (lstardis) THEN
+          CALL gradiv2(1, zu, zv, nitergdiv, zu, zv, -1.)
        ELSE
-          IF (lstardis) THEN
-             crot = 1./zllm
-          ELSE
-             crot = zllm**(-1. / nitergrot)
-          END IF
+          CALL gradiv(1, zu, zv, nitergdiv, zu, zv, -1.)
        END IF
-    END DO
 
+       zllm = max(abs(maxval(zu)), abs(maxval(zv)))
+       zu = zu / zllm
+       zv = zv / zllm
+    end DO
+
+    IF (lstardis) THEN
+       cdivu = 1. / zllm
+    ELSE
+       cdivu = zllm**(- 1. / nitergdiv)
+    END IF
     PRINT *, 'cdivu = ', cdivu
+
+    PRINT *, 'Calcul des valeurs propres de nxgrarot'
+    call random_number(zu)
+    zu = zu - 0.5
+    CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
+    call random_number(zv)
+    zv = zv - 0.5
+    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
+
+    DO l = 1, 50
+       IF (lstardis) THEN
+          CALL nxgraro2(1, zu, zv, nitergrot, zu, zv, -1.)
+       ELSE
+          CALL nxgrarot(1, zu, zv, nitergrot, zu, zv, -1.)
+       END IF
+
+       zllm = max(abs(maxval(zu)), abs(maxval(zv)))
+       zu = zu / zllm
+       zv = zv / zllm
+    end DO
+
+    IF (lstardis) THEN
+       crot = 1. / zllm
+    ELSE
+       crot = zllm**(-1. / nitergrot)
+    END IF
     PRINT *, 'crot = ', crot
-    PRINT *, 'cdivh = ', cdivh
 
     ! Variation verticale du coefficient de dissipation :
     zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
@@ -143,6 +130,7 @@
     tetaudiv = zvert / tetagdiv
     tetaurot = zvert / tetagrot
     tetah = zvert / tetatemp
+
     call new_unit(unit)
     open(unit, file="inidissip.csv", status="replace", action="write")
     write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
@@ -153,7 +141,7 @@
     print *, 'Created file "inidissip.csv".'
 
     max_zvert = maxval(zvert)
-    tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
+    tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
          tetatemp / max_zvert)
     PRINT *, 'tetamin = ', tetamin
     idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod