/[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 27 by guez, Thu Mar 25 14:29:07 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
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         DO ij = 1, ip1jmp1
55          deltap(ij,l) = 1.            deltap(ij, l) = 1.
56         ENDDO         END DO
57        ENDDO      END DO
58    
59        idum  = -1      idum = -1
60        zh(1) = RAN1(idum)-.5      zh(1) = ran1(idum) - .5
61        idum  = 0      idum = 0
62        DO ij = 2, ip1jmp1      DO ij = 2, ip1jmp1
63          zh(ij) = RAN1(idum) -.5         zh(ij) = ran1(idum) - .5
64        ENDDO      END DO
65    
66        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)      CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
67    
68        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )      CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
69    
70        IF ( zhmin .GE. zhmax  )     THEN      IF (zhmin>=zhmax) THEN
71           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax         PRINT *, '  Inidissip  zh min max  ', zhmin, zhmax
72           STOP'probleme generateur alleatoire dans inidissip'         STOP 'probleme generateur alleatoire dans inidissip'
73        ENDIF      END IF
74    
75        zllm = ABS( zhmax )      zllm = abs(zhmax)
76        DO l = 1,50      DO l = 1, 50
77           IF(lstardis) THEN         IF (lstardis) THEN
78              CALL divgrad2(1,zh,deltap,niterh,zh)            CALL divgrad2(1, zh, deltap, niterh, zh)
79           ELSE         ELSE
80              CALL divgrad (1,zh,niterh,zh)            CALL divgrad(1, zh, niterh, zh)
81           ENDIF         END IF
82    
83          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )         CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
84    
85           zllm  = ABS( zhmax )         zllm = abs(zhmax)
86           z1llm = 1./zllm         z1llm = 1./zllm
87           DO ij = 1,ip1jmp1         DO ij = 1, ip1jmp1
88              zh(ij) = zh(ij)* z1llm            zh(ij) = zh(ij)*z1llm
89           ENDDO         END DO
90        ENDDO      END DO
91    
92        IF(lstardis) THEN      IF (lstardis) THEN
93           cdivh = 1./ zllm         cdivh = 1./zllm
94        ELSE      ELSE
95           cdivh = zllm ** ( -1./niterh )         cdivh = zllm**(-1./niterh)
96        ENDIF      END IF
97    
98  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)    
99  c   -----------------------------------------------------------------  
100        print*,'calcul des valeurs propres'      PRINT *, 'calcul des valeurs propres'
101    
102        DO  20  ii = 1, 2      DO  ii = 1, 2
103  c  
104           DO ij = 1, ip1jmp1         DO ij = 1, ip1jmp1
105             zu(ij)  = RAN1(idum) -.5            zu(ij) = ran1(idum) - .5
106           ENDDO         END DO
107           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)         CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
108           DO ij = 1, ip1jm         DO ij = 1, ip1jm
109              zv(ij) = RAN1(idum) -.5            zv(ij) = ran1(idum) - .5
110           ENDDO         END DO
111           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)         CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
112    
113           CALL minmax(iip1*jjp1,zu,umin,ullm )         CALL minmax(iip1*jjp1, zu, umin, ullm)
114           CALL minmax(iip1*jjm, zv,vmin,vllm )         CALL minmax(iip1*jjm, zv, vmin, vllm)
115    
116           ullm = ABS ( ullm )         ullm = abs(ullm)
117           vllm = ABS ( vllm )         vllm = abs(vllm)
118    
119           DO  5  l = 1, 50         DO  l = 1, 50
120              IF(ii.EQ.1) THEN            IF (ii==1) THEN
121  ccccc             CALL covcont( 1,zu,zv,zu,zv )               IF (lstardis) THEN
122                 IF(lstardis) THEN                  CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
123                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )               ELSE
124                 ELSE                  CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
125                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )               END IF
126                 ENDIF            ELSE
127              ELSE               IF (lstardis) THEN
128                 IF(lstardis) THEN                  CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
129                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )               ELSE
130                 ELSE                  CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
131                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )               END IF
132                 ENDIF            END IF
133              ENDIF  
134              CALL minmax(iip1*jjp1, zu, umin, ullm)
135              CALL minmax(iip1*jjp1,zu,umin,ullm )            CALL minmax(iip1*jjm, zv, vmin, vllm)
136              CALL minmax(iip1*jjm, zv,vmin,vllm )  
137              ullm = abs(ullm)
138              ullm = ABS  ( ullm )            vllm = abs(vllm)
139              vllm = ABS  ( vllm )  
140              zllm = max(ullm, vllm)
141              zllm  = MAX( ullm,vllm )            z1llm = 1./zllm
142              z1llm = 1./ zllm            DO ij = 1, ip1jmp1
143              DO ij = 1, ip1jmp1               zu(ij) = zu(ij)*z1llm
144                zu(ij) = zu(ij)* z1llm            END DO
145              ENDDO            DO ij = 1, ip1jm
146              DO ij = 1, ip1jm               zv(ij) = zv(ij)*z1llm
147                 zv(ij) = zv(ij)* z1llm            END DO
148              ENDDO         end DO
149   5       CONTINUE  
150           IF (ii==1) THEN
151           IF ( ii.EQ.1 ) THEN            IF (lstardis) THEN
152              IF(lstardis) THEN               cdivu = 1./zllm
153                 cdivu  = 1./zllm            ELSE
154              ELSE               cdivu = zllm**(-1./nitergdiv)
155                 cdivu  = zllm **( -1./nitergdiv )            END IF
156              ENDIF         ELSE
157           ELSE            IF (lstardis) THEN
158              IF(lstardis) THEN               crot = 1./zllm
159                 crot   = 1./ zllm            ELSE
160              ELSE               crot = zllm**(-1./nitergrot)
161                 crot   = zllm **( -1./nitergrot )            END IF
162              ENDIF         END IF
163           ENDIF  
164        END DO
165   20   CONTINUE  
166        PRINT *, 'cdivu = ', cdivu
167  c   petit test pour les operateurs non star:      PRINT *, 'crot = ', crot
168  c   ----------------------------------------      PRINT *, 'cdivh = ', cdivh
169    
170  c     IF(.NOT.lstardis) THEN      ! Variation verticale du coefficient de dissipation :
171           fact    = rad*24./float(jjm)      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
172           fact    = fact*fact      ! (between 1 and 2)
173           PRINT*,'coef u ', fact/cdivu, 1./cdivu  
174           PRINT*,'coef r ', fact/crot , 1./crot      tetaudiv = zvert / tetagdiv
175           PRINT*,'coef h ', fact/cdivh, 1./cdivh      tetaurot = zvert / tetagrot
176  c     ENDIF      tetah = zvert / tetatemp
177        call new_unit(unit)
178  c-----------------------------------------------------------------------      open(unit, file="inidissip.csv", status="replace", action="write")
179  c   variation verticale du coefficient de dissipation:      write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
180  c   --------------------------------------------------      do l = 1, llm
181           write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
182        DO l=1,llm      end do
183           zvert(l)=1.      close(unit)
184        ENDDO      print *, 'Created file "inidissip.csv".'
185    
186        fact=2.      max_zvert = maxval(zvert)
187  c      tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
188        DO l = 1, llm           tetatemp / max_zvert)
189           zz      = 1. - preff/presnivs(l)      PRINT *, 'tetamin = ', tetamin
190           zvert(l)= fact -( fact-1.)/( 1.+zz*zz )      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
191        ENDDO      PRINT *, 'idissip = ', idissip
192        dtdiss = idissip * dtvr
193        PRINT *, 'dtdiss = ', dtdiss
194        PRINT*,'Constantes de temps de la diffusion horizontale'  
195      END SUBROUTINE inidissip
196        tetamin =  1.e+6  
197    end module inidissip_m
       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.27
  ViewVC Help
Powered by ViewVC 1.1.21