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