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

Diff of /trunk/Sources/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/inidissip.f90 revision 48 by guez, Tue Jul 19 12:54:20 2011 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
10  c   declarations:    integer idissip ! période de la dissipation (en pas de temps)
11  c   -------------    real tetaudiv(llm), tetaurot(llm), tetah(llm)
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      ! Initialisation de la dissipation horizontale
20    
21        LOGICAL lstardis      USE comconst, ONLY : dtvr
22        INTEGER nitergdiv,nitergrot,niterh      use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
23        REAL    tetagdiv,tetagrot,tetatemp           tetagrot, tetatemp
24        REAL fact,zvert(llm),zz      USE comvert, ONLY : preff, presnivs
25        REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)      USE conf_gcm_m, ONLY : iperiod
26        REAL ullm,vllm,umin,vmin,zhmin,zhmax      USE dimens_m, ONLY : jjm, llm
27        REAL zllm,z1llm      USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
28        use jumble, only: new_unit
29        INTEGER l,ij,idum,ii      use filtreg_m, only: filtreg
30        REAL tetamin  
31        ! Variables local to the procedure:
32        REAL ran1      REAL zvert(llm), max_zvert
33        REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
34        REAL zhmin, zhmax
35  c-----------------------------------------------------------------------      REAL zllm
36        INTEGER l, ij, idum, ii, unit
37        print *, "Call sequence information: inidissip"      REAL tetamin ! in s
38  c      REAL ran1
39  c   calcul des valeurs propres des operateurs par methode iterrative:  
40  c   -----------------------------------------------------------------      !-----------------------------------------------------------------------
41    
42        crot     = -1.      PRINT *, 'Call sequence information: inidissip'
43        cdivu    = -1.  
44        cdivh    = -1.      ! Calcul des valeurs propres des opérateurs par méthode itérative :
45    
46  c   calcul de la valeur propre de divgrad:      crot = -1.
47  c   --------------------------------------      cdivu = -1.
48        idum = 0      cdivh = -1.
49        DO l = 1, llm  
50        ! Calcul de la valeur propre de divgrad :
51    
52        deltap = 1.
53        idum = -1
54        zh(1) = ran1(idum) - 0.5
55        idum = 0
56        DO ij = 2, ip1jmp1
57           zh(ij) = ran1(idum) - 0.5
58        END DO
59    
60        CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
61    
62        CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
63        IF (zhmin >= zhmax) THEN
64           PRINT *, 'zhmin zhmax', zhmin, zhmax
65           print *, 'Problème générateur aléatoire dans inidissip'
66           STOP 1
67        END IF
68    
69        DO l = 1, 50
70           IF (lstardis) THEN
71              CALL divgrad2(1, zh, deltap, niterh, zh)
72           ELSE
73              CALL divgrad(1, zh, niterh, zh)
74           END IF
75    
76           zllm = abs(maxval(zh))
77           zh = zh / zllm
78        END DO
79    
80        IF (lstardis) THEN
81           cdivh = 1. / zllm
82        ELSE
83           cdivh = zllm**(- 1. / niterh)
84        END IF
85    
86        ! Calcul des valeurs propres de gradiv (ii = 1) et nxgrarot (ii = 2)
87    
88        PRINT *, 'Calcul des valeurs propres'
89    
90        DO ii = 1, 2
91         DO ij = 1, ip1jmp1         DO ij = 1, ip1jmp1
92          deltap(ij,l) = 1.            zu(ij) = ran1(idum) - 0.5
93         ENDDO         END DO
94        ENDDO         CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
95           DO ij = 1, ip1jm
96        idum  = -1            zv(ij) = ran1(idum) - 0.5
97        zh(1) = RAN1(idum)-.5         END DO
98        idum  = 0         CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
99        DO ij = 2, ip1jmp1  
100          zh(ij) = RAN1(idum) -.5         DO l = 1, 50
101        ENDDO            IF (ii==1) THEN
102                 IF (lstardis) THEN
103        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)                  CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
104                 ELSE
105        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )                  CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
106                 END IF
107        IF ( zhmin .GE. zhmax  )     THEN            ELSE
108           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax               IF (lstardis) THEN
109           STOP'probleme generateur alleatoire dans inidissip'                  CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
110        ENDIF               ELSE
111                    CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
112        zllm = ABS( zhmax )               END IF
113        DO l = 1,50            END IF
114           IF(lstardis) THEN  
115              CALL divgrad2(1,zh,deltap,niterh,zh)            zllm = max(abs(maxval(zu)), abs(maxval(zv)))
116           ELSE            zu = zu / zllm
117              CALL divgrad (1,zh,niterh,zh)            zv = zv / zllm
118           ENDIF         end DO
119    
120          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )         IF (ii==1) THEN
121              IF (lstardis) THEN
122           zllm  = ABS( zhmax )               cdivu = 1. / zllm
123           z1llm = 1./zllm            ELSE
124           DO ij = 1,ip1jmp1               cdivu = zllm**(- 1. / nitergdiv)
125              zh(ij) = zh(ij)* z1llm            END IF
126           ENDDO         ELSE
127        ENDDO            IF (lstardis) THEN
128                 crot = 1./zllm
129        IF(lstardis) THEN            ELSE
130           cdivh = 1./ zllm               crot = zllm**(-1. / nitergrot)
131        ELSE            END IF
132           cdivh = zllm ** ( -1./niterh )         END IF
133        ENDIF      END DO
134    
135  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)      PRINT *, 'cdivu = ', cdivu
136  c   -----------------------------------------------------------------      PRINT *, 'crot = ', crot
137        print*,'calcul des valeurs propres'      PRINT *, 'cdivh = ', cdivh
138    
139        DO  20  ii = 1, 2      ! Variation verticale du coefficient de dissipation :
140  c      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
141           DO ij = 1, ip1jmp1      ! (between 1 and 2)
142             zu(ij)  = RAN1(idum) -.5  
143           ENDDO      tetaudiv = zvert / tetagdiv
144           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)      tetaurot = zvert / tetagrot
145           DO ij = 1, ip1jm      tetah = zvert / tetatemp
146              zv(ij) = RAN1(idum) -.5      call new_unit(unit)
147           ENDDO      open(unit, file="inidissip.csv", status="replace", action="write")
148           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)      write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
149        do l = 1, llm
150           CALL minmax(iip1*jjp1,zu,umin,ullm )         write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
151           CALL minmax(iip1*jjm, zv,vmin,vllm )      end do
152        close(unit)
153           ullm = ABS ( ullm )      print *, 'Created file "inidissip.csv".'
154           vllm = ABS ( vllm )  
155        max_zvert = maxval(zvert)
156           DO  5  l = 1, 50      tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
157              IF(ii.EQ.1) THEN           tetatemp / max_zvert)
158  ccccc             CALL covcont( 1,zu,zv,zu,zv )      PRINT *, 'tetamin = ', tetamin
159                 IF(lstardis) THEN      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
160                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )      PRINT *, 'idissip = ', idissip
161                 ELSE      dtdiss = idissip * dtvr
162                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )      PRINT *, 'dtdiss = ', dtdiss
163                 ENDIF  
164              ELSE    END SUBROUTINE inidissip
165                 IF(lstardis) THEN  
166                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )  end module inidissip_m
                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  
Legend:
Removed from v.25  
changed lines
  Added in v.48
  ViewVC Help
Powered by ViewVC 1.1.21