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