--- trunk/libf/dyn3d/inidissip.f	2010/03/05 16:43:45	25
+++ trunk/libf/dyn3d/inidissip.f90	2011/12/06 15:07:04	54
@@ -1,218 +1,154 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/inidissip.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $
-!
-      SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
-     *                       tetagdiv,tetagrot,tetatemp             )
-c=======================================================================
-c   initialisation de la dissipation horizontale
-c=======================================================================
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      use dimens_m
-      use paramet_m
-      use comconst
-      use comvert
-      use conf_gcm_m
-            use comdissipn
-      IMPLICIT NONE
-
-      LOGICAL lstardis
-      INTEGER nitergdiv,nitergrot,niterh
-      REAL    tetagdiv,tetagrot,tetatemp
-      REAL fact,zvert(llm),zz
-      REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
-      REAL ullm,vllm,umin,vmin,zhmin,zhmax
-      REAL zllm,z1llm
-
-      INTEGER l,ij,idum,ii
-      REAL tetamin
-
-      REAL ran1
-
-
-c-----------------------------------------------------------------------
-
-      print *, "Call sequence information: inidissip"
-c
-c   calcul des valeurs propres des operateurs par methode iterrative:
-c   -----------------------------------------------------------------
-
-      crot     = -1.
-      cdivu    = -1.
-      cdivh    = -1.
-
-c   calcul de la valeur propre de divgrad:
-c   --------------------------------------
-      idum = 0
-      DO l = 1, llm
-       DO ij = 1, ip1jmp1
-        deltap(ij,l) = 1.
-       ENDDO
-      ENDDO
-
-      idum  = -1
-      zh(1) = RAN1(idum)-.5
-      idum  = 0
-      DO ij = 2, ip1jmp1
-        zh(ij) = RAN1(idum) -.5
-      ENDDO
-
-      CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
-
-      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-      IF ( zhmin .GE. zhmax  )     THEN
-         PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
-         STOP'probleme generateur alleatoire dans inidissip'
-      ENDIF
-
-      zllm = ABS( zhmax )
-      DO l = 1,50
-         IF(lstardis) THEN
-            CALL divgrad2(1,zh,deltap,niterh,zh)
-         ELSE
-            CALL divgrad (1,zh,niterh,zh)
-         ENDIF
-
-        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-         zllm  = ABS( zhmax )
-         z1llm = 1./zllm
-         DO ij = 1,ip1jmp1
-            zh(ij) = zh(ij)* z1llm
-         ENDDO
-      ENDDO
-
-      IF(lstardis) THEN
-         cdivh = 1./ zllm
-      ELSE
-         cdivh = zllm ** ( -1./niterh )
-      ENDIF
-
-c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
-c   -----------------------------------------------------------------
-      print*,'calcul des valeurs propres'
-
-      DO  20  ii = 1, 2
-c
-         DO ij = 1, ip1jmp1
-           zu(ij)  = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
-         DO ij = 1, ip1jm
-            zv(ij) = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
-
-         CALL minmax(iip1*jjp1,zu,umin,ullm )
-         CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-         ullm = ABS ( ullm )
-         vllm = ABS ( vllm )
-
-         DO  5  l = 1, 50
-            IF(ii.EQ.1) THEN
-ccccc             CALL covcont( 1,zu,zv,zu,zv )
-               IF(lstardis) THEN
-                  CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
-               ELSE
-                  CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
-               ENDIF
-            ELSE
-               IF(lstardis) THEN
-                  CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
-               ELSE
-                  CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
-               ENDIF
-            ENDIF
-
-            CALL minmax(iip1*jjp1,zu,umin,ullm )
-            CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-            ullm = ABS  ( ullm )
-            vllm = ABS  ( vllm )
-
-            zllm  = MAX( ullm,vllm )
-            z1llm = 1./ zllm
-            DO ij = 1, ip1jmp1
-              zu(ij) = zu(ij)* z1llm
-            ENDDO
-            DO ij = 1, ip1jm
-               zv(ij) = zv(ij)* z1llm
-            ENDDO
- 5       CONTINUE
-
-         IF ( ii.EQ.1 ) THEN
-            IF(lstardis) THEN
-               cdivu  = 1./zllm
-            ELSE
-               cdivu  = zllm **( -1./nitergdiv )
-            ENDIF
-         ELSE
-            IF(lstardis) THEN
-               crot   = 1./ zllm
-            ELSE
-               crot   = zllm **( -1./nitergrot )
-            ENDIF
-         ENDIF
-
- 20   CONTINUE
-
-c   petit test pour les operateurs non star:
-c   ----------------------------------------
-
-c     IF(.NOT.lstardis) THEN
-         fact    = rad*24./float(jjm)
-         fact    = fact*fact
-         PRINT*,'coef u ', fact/cdivu, 1./cdivu
-         PRINT*,'coef r ', fact/crot , 1./crot
-         PRINT*,'coef h ', fact/cdivh, 1./cdivh
-c     ENDIF
-
-c-----------------------------------------------------------------------
-c   variation verticale du coefficient de dissipation:
-c   --------------------------------------------------
-
-      DO l=1,llm
-         zvert(l)=1.
-      ENDDO
-
-      fact=2.
-c
-      DO l = 1, llm
-         zz      = 1. - preff/presnivs(l)
-         zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
-      ENDDO
-
-
-      PRINT*,'Constantes de temps de la diffusion horizontale'
-
-      tetamin =  1.e+6
-
-      DO l=1,llm
-        tetaudiv(l)   = zvert(l)/tetagdiv
-        tetaurot(l)   = zvert(l)/tetagrot
-        tetah(l)      = zvert(l)/tetatemp
-
-        IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
-        IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
-        IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
-      ENDDO
-
-      PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
-      idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
-      PRINT *,' tetamin = ',tetamin
-      idissip = MAX(iperiod,idissip)
-      dtdiss  = idissip * dtvr
-      PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
-
-      DO l = 1,llm
-         PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
-     *                   dtdiss*tetah(l)
-      ENDDO
-
-c
-      RETURN
-      END
+module inidissip_m
+
+  use dimens_m, only: llm
+
+  IMPLICIT NONE
+
+  private llm
+
+  REAL dtdiss
+  integer idissip ! période de la dissipation (en pas de temps)
+  real tetaudiv(llm), tetaurot(llm), tetah(llm) 
+  real cdivu, crot, cdivh
+
+contains
+
+  SUBROUTINE inidissip
+
+    ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
+
+    ! 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 : 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, dimension(iim + 1, jjm + 1):: zh, zu
+    real zv(iim + 1, jjm), deltap(iim + 1, jjm + 1, llm)
+    REAL zllm
+    INTEGER l, seed_size, ii, unit
+    REAL tetamin ! in s
+
+    !-----------------------------------------------------------------------
+
+    PRINT *, 'Call sequence information: inidissip'
+    call random_seed(size=seed_size)
+    call random_seed(put=(/(0, ii = 1, seed_size)/))
+
+    PRINT *, 'Calcul des valeurs propres de divgrad'
+    deltap = 1.
+    call random_number(zh)
+    zh = zh - 0.5
+    CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
+
+    DO l = 1, 50
+       IF (lstardis) THEN
+          CALL divgrad2(1, zh, deltap, niterh, zh, -1.)
+       ELSE
+          CALL divgrad(1, zh, niterh, zh, -1.)
+       END IF
+
+       zllm = abs(maxval(zh))
+       zh = zh / zllm
+    END DO
+
+    IF (lstardis) THEN
+       cdivh = 1. / zllm
+    ELSE
+       cdivh = zllm**(- 1. / niterh)
+    END IF
+    PRINT *, 'cdivh = ', cdivh
+
+    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 l = 1, 50
+       IF (lstardis) THEN
+          CALL gradiv2(1, zu, zv, nitergdiv, zu, zv, -1.)
+       ELSE
+          CALL gradiv(1, zu, zv, nitergdiv, zu, zv, -1.)
+       END IF
+
+       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
+
+    ! Variation verticale du coefficient de dissipation :
+    zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
+    ! (between 1 and 2)
+
+    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
+    do l = 1, llm
+       write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
+    end do
+    close(unit)
+    print *, 'Created file "inidissip.csv".'
+
+    max_zvert = maxval(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
+    PRINT *, 'idissip = ', idissip
+    dtdiss = idissip * dtvr
+    PRINT *, 'dtdiss = ', dtdiss
+
+  END SUBROUTINE inidissip
+
+end module inidissip_m