/[lmdze]/trunk/dyn3d/Dissipation/inidissip.f
ViewVC logotype

Diff of /trunk/dyn3d/Dissipation/inidissip.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/inidissip.f revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC trunk/libf/dyn3d/Dissipation/inidissip.f90 revision 64 by guez, Wed Aug 29 14:47:17 2012 UTC
# Line 1  Line 1 
1  !  module inidissip_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/inidissip.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $  
3  !    use dimens_m, only: llm
4        SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,  
5       *                       tetagdiv,tetagrot,tetatemp             )    IMPLICIT NONE
6  c=======================================================================  
7  c   initialisation de la dissipation horizontale    private llm
8  c=======================================================================  
9  c-----------------------------------------------------------------------    REAL dtdiss ! in s
10  c   declarations:    integer idissip ! période de la dissipation (en pas de temps)
11  c   -------------    real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s
12      real cdivu, crot, cdivh
13        use dimens_m  
14        use paramet_m  contains
15        use comconst  
16        use comvert    SUBROUTINE inidissip
17        use conf_gcm_m  
18              use comdissipn      ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
19        IMPLICIT NONE  
20        ! Initialisation de la dissipation horizontale. Calcul des valeurs
21        LOGICAL lstardis      ! propres des opérateurs par méthode itérative.
22        INTEGER nitergdiv,nitergrot,niterh  
23        REAL    tetagdiv,tetagrot,tetatemp      USE comconst, ONLY: dtvr
24        REAL fact,zvert(llm),zz      use comdissnew, only: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
25        REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)           tetatemp
26        REAL ullm,vllm,umin,vmin,zhmin,zhmax      USE comvert, ONLY: preff, presnivs
27        REAL zllm,z1llm      USE conf_gcm_m, ONLY: iperiod
28        USE dimens_m, ONLY: iim, jjm
29        INTEGER l,ij,idum,ii      use filtreg_m, only: filtreg
30        REAL tetamin      use gradiv2_m, only: gradiv2
31        use jumble, only: new_unit
32        REAL ran1      USE paramet_m, ONLY: jjp1
33    
34        ! Variables local to the procedure:
35  c-----------------------------------------------------------------------      REAL zvert(llm), max_zvert ! no dimension
36        REAL, dimension(iim + 1, jjm + 1, 1):: zh, zu, gx, divgra, deltap
37        print *, "Call sequence information: inidissip"      real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1)
38  c      REAL zllm
39  c   calcul des valeurs propres des operateurs par methode iterrative:      INTEGER l, seed_size, ii, unit
40  c   -----------------------------------------------------------------      REAL tetamin ! in s
41    
42        crot     = -1.      !-----------------------------------------------------------------------
43        cdivu    = -1.  
44        cdivh    = -1.      PRINT *, 'Call sequence information: inidissip'
45        call random_seed(size=seed_size)
46  c   calcul de la valeur propre de divgrad:      call random_seed(put=(/(0, ii = 1, seed_size)/))
47  c   --------------------------------------  
48        idum = 0      PRINT *, 'Calcul des valeurs propres de divgrad'
49        DO l = 1, llm      deltap = 1.
50         DO ij = 1, ip1jmp1      call random_number(zh)
51          deltap(ij,l) = 1.      zh = zh - 0.5
52         ENDDO      CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE.)
53        ENDDO  
54        DO l = 1, 50
55        idum  = -1         CALL divgrad2(1, zh, deltap, niterh, divgra, -1.)
56        zh(1) = RAN1(idum)-.5         zllm = abs(maxval(divgra))
57        idum  = 0         zh = divgra / zllm
58        DO ij = 2, ip1jmp1      END DO
59          zh(ij) = RAN1(idum) -.5  
60        ENDDO      cdivh = 1. / zllm
61        PRINT *, 'cdivh = ', cdivh
62        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)  
63        PRINT *, 'Calcul des valeurs propres de gradiv'
64        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )      call random_number(zu)
65        zu = zu - 0.5
66        IF ( zhmin .GE. zhmax  )     THEN      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
67           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax      call random_number(zv)
68           STOP'probleme generateur alleatoire dans inidissip'      zv = zv - 0.5
69        ENDIF      CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
70    
71        zllm = ABS( zhmax )      DO l = 1, 50
72        DO l = 1,50         CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.)
73           IF(lstardis) THEN         zllm = max(abs(maxval(gx)), abs(maxval(gy)))
74              CALL divgrad2(1,zh,deltap,niterh,zh)         zu = gx / zllm
75           ELSE         zv = gy / zllm
76              CALL divgrad (1,zh,niterh,zh)      end DO
77           ENDIF  
78        cdivu = 1. / zllm
79          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )      PRINT *, 'cdivu = ', cdivu
80    
81           zllm  = ABS( zhmax )      PRINT *, 'Calcul des valeurs propres de nxgrarot'
82           z1llm = 1./zllm      call random_number(zu)
83           DO ij = 1,ip1jmp1      zu = zu - 0.5
84              zh(ij) = zh(ij)* z1llm      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
85           ENDDO      call random_number(zv)
86        ENDDO      zv = zv - 0.5
87        CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
88        IF(lstardis) THEN  
89           cdivh = 1./ zllm      DO l = 1, 50
90        ELSE         CALL nxgraro2(1, zu, zv, nitergrot, gx, gy, -1.)
91           cdivh = zllm ** ( -1./niterh )         zllm = max(abs(maxval(gx)), abs(maxval(gy)))
92        ENDIF         zu = gx / zllm
93           zv = gy / zllm
94  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)      end DO
95  c   -----------------------------------------------------------------  
96        print*,'calcul des valeurs propres'      crot = 1. / zllm
97        PRINT *, 'crot = ', crot
98        DO  20  ii = 1, 2  
99  c      ! Variation verticale du coefficient de dissipation :
100           DO ij = 1, ip1jmp1      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
101             zu(ij)  = RAN1(idum) -.5      ! (between 1 and 2)
102           ENDDO  
103           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)      tetaudiv = zvert / tetagdiv
104           DO ij = 1, ip1jm      tetaurot = zvert / tetagrot
105              zv(ij) = RAN1(idum) -.5      tetah = zvert / tetatemp
106           ENDDO  
107           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)      call new_unit(unit)
108        open(unit, file="inidissip.csv", status="replace", action="write")
109           CALL minmax(iip1*jjp1,zu,umin,ullm )      write(unit, fmt=*) '"tetaudiv (s)" "tetaurot (s)" "tetah (s)"' ! title line
110           CALL minmax(iip1*jjm, zv,vmin,vllm )      do l = 1, llm
111           write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
112           ullm = ABS ( ullm )      end do
113           vllm = ABS ( vllm )      close(unit)
114        print *, 'Created file "inidissip.csv".'
115           DO  5  l = 1, 50  
116              IF(ii.EQ.1) THEN      max_zvert = maxval(zvert)
117  ccccc             CALL covcont( 1,zu,zv,zu,zv )      tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
118                 IF(lstardis) THEN           tetatemp / max_zvert)
119                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )      PRINT *, 'tetamin = ', tetamin
120                 ELSE      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
121                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )      PRINT *, 'idissip = ', idissip
122                 ENDIF      dtdiss = idissip * dtvr
123              ELSE      PRINT *, 'dtdiss = ', dtdiss, "s"
124                 IF(lstardis) THEN  
125                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )    END SUBROUTINE inidissip
126                 ELSE  
127                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )  end module inidissip_m
                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  
Legend:
Removed from v.25  
changed lines
  Added in v.64
  ViewVC Help
Powered by ViewVC 1.1.21