/[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 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/libf/dyn3d/inidissip.f90 revision 39 by guez, Tue Jan 25 15:11:05 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 ! periode de la dissipation (en pas)
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        IMPLICIT NONE      ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
19        include "comdissipn.h"      ! 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 new_unit_m, 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 ullm, vllm, umin, vmin, zhmin, zhmax
35  c-----------------------------------------------------------------------      REAL zllm, z1llm
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 operateurs par methode iterrative:  
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        idum = 0
53        DO l = 1, llm
54           DO ij = 1, ip1jmp1
55              deltap(ij, l) = 1.
56           END DO
57        END DO
58    
59        idum = -1
60        zh(1) = ran1(idum) - .5
61        idum = 0
62        DO ij = 2, ip1jmp1
63           zh(ij) = ran1(idum) - .5
64        END DO
65    
66        CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
67    
68        CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
69        IF (zhmin >= zhmax) THEN
70           PRINT *, 'zhmin zhmax', zhmin, zhmax
71           print *, 'Problème générateur aléatoire dans inidissip'
72           STOP 1
73        END IF
74    
75        zllm = abs(zhmax)
76        DO l = 1, 50
77           IF (lstardis) THEN
78              CALL divgrad2(1, zh, deltap, niterh, zh)
79           ELSE
80              CALL divgrad(1, zh, niterh, zh)
81           END IF
82    
83           CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
84    
85           zllm = abs(zhmax)
86           z1llm = 1./zllm
87           DO ij = 1, ip1jmp1
88              zh(ij) = zh(ij)*z1llm
89           END DO
90        END DO
91    
92        IF (lstardis) THEN
93           cdivh = 1./zllm
94        ELSE
95           cdivh = zllm**(-1./niterh)
96        END IF
97    
98        !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)    
99    
100        PRINT *, 'calcul des valeurs propres'
101    
102        DO  ii = 1, 2
103         DO ij = 1, ip1jmp1         DO ij = 1, ip1jmp1
104          deltap(ij,l) = 1.            zu(ij) = ran1(idum) - .5
105         ENDDO         END DO
106        ENDDO         CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
107           DO ij = 1, ip1jm
108        idum  = -1            zv(ij) = ran1(idum) - .5
109        zh(1) = RAN1(idum)-.5         END DO
110        idum  = 0         CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
111        DO ij = 2, ip1jmp1  
112          zh(ij) = RAN1(idum) -.5         CALL minmax(iip1*jjp1, zu, umin, ullm)
113        ENDDO         CALL minmax(iip1*jjm, zv, vmin, vllm)
114    
115        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)         ullm = abs(ullm)
116           vllm = abs(vllm)
117        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
118           DO  l = 1, 50
119        IF ( zhmin .GE. zhmax  )     THEN            IF (ii==1) THEN
120           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax               IF (lstardis) THEN
121           STOP'probleme generateur alleatoire dans inidissip'                  CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
122        ENDIF               ELSE
123                    CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
124        zllm = ABS( zhmax )               END IF
125        DO l = 1,50            ELSE
126           IF(lstardis) THEN               IF (lstardis) THEN
127              CALL divgrad2(1,zh,deltap,niterh,zh)                  CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
128           ELSE               ELSE
129              CALL divgrad (1,zh,niterh,zh)                  CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
130           ENDIF               END IF
131              END IF
132          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
133              CALL minmax(iip1*jjp1, zu, umin, ullm)
134           zllm  = ABS( zhmax )            CALL minmax(iip1*jjm, zv, vmin, vllm)
135           z1llm = 1./zllm  
136           DO ij = 1,ip1jmp1            ullm = abs(ullm)
137              zh(ij) = zh(ij)* z1llm            vllm = abs(vllm)
138           ENDDO  
139        ENDDO            zllm = max(ullm, vllm)
140              z1llm = 1./zllm
141        IF(lstardis) THEN            DO ij = 1, ip1jmp1
142           cdivh = 1./ zllm               zu(ij) = zu(ij)*z1llm
143        ELSE            END DO
144           cdivh = zllm ** ( -1./niterh )            DO ij = 1, ip1jm
145        ENDIF               zv(ij) = zv(ij)*z1llm
146              END DO
147  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)         end DO
148  c   -----------------------------------------------------------------  
149        print*,'calcul des valeurs propres'         IF (ii==1) THEN
150              IF (lstardis) THEN
151        DO  20  ii = 1, 2               cdivu = 1./zllm
152  c            ELSE
153           DO ij = 1, ip1jmp1               cdivu = zllm**(-1./nitergdiv)
154             zu(ij)  = RAN1(idum) -.5            END IF
155           ENDDO         ELSE
156           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)            IF (lstardis) THEN
157           DO ij = 1, ip1jm               crot = 1./zllm
158              zv(ij) = RAN1(idum) -.5            ELSE
159           ENDDO               crot = zllm**(-1./nitergrot)
160           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)            END IF
161           END IF
162           CALL minmax(iip1*jjp1,zu,umin,ullm )      END DO
163           CALL minmax(iip1*jjm, zv,vmin,vllm )  
164        PRINT *, 'cdivu = ', cdivu
165           ullm = ABS ( ullm )      PRINT *, 'crot = ', crot
166           vllm = ABS ( vllm )      PRINT *, 'cdivh = ', cdivh
167    
168           DO  5  l = 1, 50      ! Variation verticale du coefficient de dissipation :
169              IF(ii.EQ.1) THEN      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
170  ccccc             CALL covcont( 1,zu,zv,zu,zv )      ! (between 1 and 2)
171                 IF(lstardis) THEN  
172                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )      tetaudiv = zvert / tetagdiv
173                 ELSE      tetaurot = zvert / tetagrot
174                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )      tetah = zvert / tetatemp
175                 ENDIF      call new_unit(unit)
176              ELSE      open(unit, file="inidissip.csv", status="replace", action="write")
177                 IF(lstardis) THEN      write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
178                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )      do l = 1, llm
179                 ELSE         write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
180                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )      end do
181                 ENDIF      close(unit)
182              ENDIF      print *, 'Created file "inidissip.csv".'
183    
184              CALL minmax(iip1*jjp1,zu,umin,ullm )      max_zvert = maxval(zvert)
185              CALL minmax(iip1*jjm, zv,vmin,vllm )      tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
186             tetatemp / max_zvert)
187              ullm = ABS  ( ullm )      PRINT *, 'tetamin = ', tetamin
188              vllm = ABS  ( vllm )      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
189        PRINT *, 'idissip = ', idissip
190              zllm  = MAX( ullm,vllm )      dtdiss = idissip * dtvr
191              z1llm = 1./ zllm      PRINT *, 'dtdiss = ', dtdiss
192              DO ij = 1, ip1jmp1  
193                zu(ij) = zu(ij)* z1llm    END SUBROUTINE inidissip
194              ENDDO  
195              DO ij = 1, ip1jm  end module inidissip_m
                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 *,' INIDI tetamin idissip ',tetamin,idissip  
       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.3  
changed lines
  Added in v.39
  ViewVC Help
Powered by ViewVC 1.1.21