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