/[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 26 by guez, Tue Mar 9 15:27:15 2010 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(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
17        use conf_gcm_m         tetagrot, tetatemp)
18              use comdissipn  
19        IMPLICIT NONE      ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
20        ! Initialisation de la dissipation horizontale                        
21        LOGICAL lstardis  
22        INTEGER nitergdiv,nitergrot,niterh      USE comconst, ONLY : dtvr
23        REAL    tetagdiv,tetagrot,tetatemp      USE comvert, ONLY : preff, presnivs
24        REAL fact,zvert(llm),zz      USE conf_gcm_m, ONLY : iperiod
25        REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)      USE dimens_m, ONLY : jjm, llm
26        REAL ullm,vllm,umin,vmin,zhmin,zhmax      USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
27        REAL zllm,z1llm      use new_unit_m, only: new_unit
28    
29        INTEGER l,ij,idum,ii      LOGICAL, intent(in):: lstardis
30        REAL tetamin      INTEGER, intent(in):: nitergdiv, nitergrot, niterh
31        REAL, intent(in):: tetagdiv, tetagrot, tetatemp
32        REAL ran1  
33        ! Variables local to the procedure:
34        REAL zvert(llm)
35  c-----------------------------------------------------------------------      REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
36        REAL ullm, vllm, umin, vmin, zhmin, zhmax
37        print *, "Call sequence information: inidissip"      REAL zllm, z1llm
38  c      INTEGER l, ij, idum, ii, unit
39  c   calcul des valeurs propres des operateurs par methode iterrative:      REAL tetamin
40  c   -----------------------------------------------------------------      REAL ran1
41    
42        crot     = -1.      !-----------------------------------------------------------------------
43        cdivu    = -1.  
44        cdivh    = -1.      PRINT *, 'Call sequence information: inidissip'
45    
46  c   calcul de la valeur propre de divgrad:      !   calcul des valeurs propres des operateurs par methode iterrative:  
47  c   --------------------------------------  
48        idum = 0      crot = -1.
49        DO l = 1, llm      cdivu = -1.
50        cdivh = -1.
51    
52        !   calcul de la valeur propre de divgrad:                              
53    
54        idum = 0
55        DO l = 1, llm
56         DO ij = 1, ip1jmp1         DO ij = 1, ip1jmp1
57          deltap(ij,l) = 1.            deltap(ij, l) = 1.
58         ENDDO         END DO
59        ENDDO      END DO
60    
61        idum  = -1      idum = -1
62        zh(1) = RAN1(idum)-.5      zh(1) = ran1(idum) - .5
63        idum  = 0      idum = 0
64        DO ij = 2, ip1jmp1      DO ij = 2, ip1jmp1
65          zh(ij) = RAN1(idum) -.5         zh(ij) = ran1(idum) - .5
66        ENDDO      END DO
67    
68        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)      CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
69    
70        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )      CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
71    
72        IF ( zhmin .GE. zhmax  )     THEN      IF (zhmin>=zhmax) THEN
73           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax         PRINT *, '  Inidissip  zh min max  ', zhmin, zhmax
74           STOP'probleme generateur alleatoire dans inidissip'         STOP 'probleme generateur alleatoire dans inidissip'
75        ENDIF      END IF
76    
77        zllm = ABS( zhmax )      zllm = abs(zhmax)
78        DO l = 1,50      DO l = 1, 50
79           IF(lstardis) THEN         IF (lstardis) THEN
80              CALL divgrad2(1,zh,deltap,niterh,zh)            CALL divgrad2(1, zh, deltap, niterh, zh)
81           ELSE         ELSE
82              CALL divgrad (1,zh,niterh,zh)            CALL divgrad(1, zh, niterh, zh)
83           ENDIF         END IF
84    
85          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )         CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
86    
87           zllm  = ABS( zhmax )         zllm = abs(zhmax)
88           z1llm = 1./zllm         z1llm = 1./zllm
89           DO ij = 1,ip1jmp1         DO ij = 1, ip1jmp1
90              zh(ij) = zh(ij)* z1llm            zh(ij) = zh(ij)*z1llm
91           ENDDO         END DO
92        ENDDO      END DO
93    
94        IF(lstardis) THEN      IF (lstardis) THEN
95           cdivh = 1./ zllm         cdivh = 1./zllm
96        ELSE      ELSE
97           cdivh = zllm ** ( -1./niterh )         cdivh = zllm**(-1./niterh)
98        ENDIF      END IF
99    
100  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)      !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)    
101  c   -----------------------------------------------------------------  
102        print*,'calcul des valeurs propres'      PRINT *, 'calcul des valeurs propres'
103    
104        DO  20  ii = 1, 2      DO  ii = 1, 2
105  c  
106           DO ij = 1, ip1jmp1         DO ij = 1, ip1jmp1
107             zu(ij)  = RAN1(idum) -.5            zu(ij) = ran1(idum) - .5
108           ENDDO         END DO
109           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)         CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
110           DO ij = 1, ip1jm         DO ij = 1, ip1jm
111              zv(ij) = RAN1(idum) -.5            zv(ij) = ran1(idum) - .5
112           ENDDO         END DO
113           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)         CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
114    
115           CALL minmax(iip1*jjp1,zu,umin,ullm )         CALL minmax(iip1*jjp1, zu, umin, ullm)
116           CALL minmax(iip1*jjm, zv,vmin,vllm )         CALL minmax(iip1*jjm, zv, vmin, vllm)
117    
118           ullm = ABS ( ullm )         ullm = abs(ullm)
119           vllm = ABS ( vllm )         vllm = abs(vllm)
120    
121           DO  5  l = 1, 50         DO  l = 1, 50
122              IF(ii.EQ.1) THEN            IF (ii==1) THEN
123  ccccc             CALL covcont( 1,zu,zv,zu,zv )               IF (lstardis) THEN
124                 IF(lstardis) THEN                  CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
125                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )               ELSE
126                 ELSE                  CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
127                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )               END IF
128                 ENDIF            ELSE
129              ELSE               IF (lstardis) THEN
130                 IF(lstardis) THEN                  CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
131                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )               ELSE
132                 ELSE                  CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
133                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )               END IF
134                 ENDIF            END IF
135              ENDIF  
136              CALL minmax(iip1*jjp1, zu, umin, ullm)
137              CALL minmax(iip1*jjp1,zu,umin,ullm )            CALL minmax(iip1*jjm, zv, vmin, vllm)
138              CALL minmax(iip1*jjm, zv,vmin,vllm )  
139              ullm = abs(ullm)
140              ullm = ABS  ( ullm )            vllm = abs(vllm)
141              vllm = ABS  ( vllm )  
142              zllm = max(ullm, vllm)
143              zllm  = MAX( ullm,vllm )            z1llm = 1./zllm
144              z1llm = 1./ zllm            DO ij = 1, ip1jmp1
145              DO ij = 1, ip1jmp1               zu(ij) = zu(ij)*z1llm
146                zu(ij) = zu(ij)* z1llm            END DO
147              ENDDO            DO ij = 1, ip1jm
148              DO ij = 1, ip1jm               zv(ij) = zv(ij)*z1llm
149                 zv(ij) = zv(ij)* z1llm            END DO
150              ENDDO         end DO
151   5       CONTINUE  
152           IF (ii==1) THEN
153           IF ( ii.EQ.1 ) THEN            IF (lstardis) THEN
154              IF(lstardis) THEN               cdivu = 1./zllm
155                 cdivu  = 1./zllm            ELSE
156              ELSE               cdivu = zllm**(-1./nitergdiv)
157                 cdivu  = zllm **( -1./nitergdiv )            END IF
158              ENDIF         ELSE
159           ELSE            IF (lstardis) THEN
160              IF(lstardis) THEN               crot = 1./zllm
161                 crot   = 1./ zllm            ELSE
162              ELSE               crot = zllm**(-1./nitergrot)
163                 crot   = zllm **( -1./nitergrot )            END IF
164              ENDIF         END IF
165           ENDIF  
166        END DO
167   20   CONTINUE  
168        PRINT *, 'cdivu = ', cdivu
169  c   petit test pour les operateurs non star:      PRINT *, 'crot = ', crot
170  c   ----------------------------------------      PRINT *, 'cdivh = ', cdivh
171    
172  c     IF(.NOT.lstardis) THEN      ! Variation verticale du coefficient de dissipation :
173           fact    = rad*24./float(jjm)      zvert = 2. - 1. / (1. + (1. - preff / presnivs)**2)
174           fact    = fact*fact  
175           PRINT*,'coef u ', fact/cdivu, 1./cdivu      tetaudiv = zvert / tetagdiv
176           PRINT*,'coef r ', fact/crot , 1./crot      tetaurot = zvert / tetagrot
177           PRINT*,'coef h ', fact/cdivh, 1./cdivh      tetah = zvert / tetatemp
178  c     ENDIF      call new_unit(unit)
179        open(unit, file="inidissip.csv", status="replace", action="write")
180  c-----------------------------------------------------------------------      write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
181  c   variation verticale du coefficient de dissipation:      do l = 1, llm
182  c   --------------------------------------------------         write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
183        end do
184        DO l=1,llm      close(unit)
185           zvert(l)=1.      print *, 'Created file "inidissip.csv".'
186        ENDDO  
187        tetamin = min(1E6, minval(1. / tetaudiv), minval(1. / tetaurot), &
188        fact=2.           minval(1. / tetah))
189  c      PRINT *, 'tetamin = ', tetamin
190        DO l = 1, llm      idissip = max(iperiod, int(tetamin / (2 * dtvr * iperiod)) * iperiod)
191           zz      = 1. - preff/presnivs(l)      PRINT *, 'idissip = ', idissip
192           zvert(l)= fact -( fact-1.)/( 1.+zz*zz )      dtdiss = idissip * dtvr
193        ENDDO      PRINT *, 'dtdiss = ', dtdiss
194    
195      END SUBROUTINE inidissip
196        PRINT*,'Constantes de temps de la diffusion horizontale'  
197    end module inidissip_m
       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.26
  ViewVC Help
Powered by ViewVC 1.1.21