/[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 26 by guez, Tue Mar 9 15:27:15 2010 UTC trunk/libf/dyn3d/Dissipation/inidissip.f90 revision 65 by guez, Thu Sep 20 09:57:03 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-1
12    real cdivu, crot, cdivh    real cdivu, crot, cdivh
13    
14  contains  contains
15    
16    SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &    SUBROUTINE inidissip
        tetagrot, tetatemp)  
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 comvert, ONLY : preff, presnivs      ! propres des opérateurs par méthode itérative.
22      USE conf_gcm_m, ONLY : iperiod  
23      USE dimens_m, ONLY : jjm, llm      USE comconst, ONLY: dtvr
24      USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1      use comdissnew, only: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
25      use new_unit_m, only: new_unit           tetatemp
26        USE comvert, ONLY: preff, presnivs
27      LOGICAL, intent(in):: lstardis      USE conf_gcm_m, ONLY: iperiod
28      INTEGER, intent(in):: nitergdiv, nitergrot, niterh      USE dimens_m, ONLY: iim, jjm
29      REAL, intent(in):: tetagdiv, tetagrot, tetatemp      use divgrad2_m, only: divgrad2
30        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)      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 ullm, vllm, umin, vmin, zhmin, zhmax      real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1)
40      REAL zllm, z1llm      REAL zllm
41      INTEGER l, ij, idum, ii, unit      INTEGER l, seed_size, ii, unit
42      REAL tetamin      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=(/(0, ii = 1, seed_size)/))
49    
50      !   calcul des valeurs propres des operateurs par methode iterrative:        PRINT *, 'Calcul des valeurs propres de divgrad'
51        deltap = 1.
52      crot = -1.      call random_number(zh)
53      cdivu = -1.      zh = zh - 0.5
54      cdivh = -1.      CALL filtreg(zh, jjp1, 1, 2, 1, .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  
55    
56      idum = -1      DO l = 1, 50
57      zh(1) = ran1(idum) - .5         CALL divgrad2(1, zh, deltap, niterh, divgra, -1.)
58      idum = 0         zllm = abs(maxval(divgra))
59      DO ij = 2, ip1jmp1         zh = divgra / zllm
        zh(ij) = ran1(idum) - .5  
60      END DO      END DO
61    
62      CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)      cdivh = 1. / zllm
63        PRINT *, 'cdivh = ', cdivh
     CALL minmax(iip1*jjp1, zh, zhmin, zhmax)  
64    
65      IF (zhmin>=zhmax) THEN      PRINT *, 'Calcul des valeurs propres de gradiv'
66         PRINT *, '  Inidissip  zh min max  ', zhmin, zhmax      call random_number(zu)
67         STOP 'probleme generateur alleatoire dans inidissip'      zu = zu - 0.5
68      END IF      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
69        call random_number(zv)
70        zv = zv - 0.5
71        CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
72    
     zllm = abs(zhmax)  
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
   
        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  
79    
80      IF (lstardis) THEN      cdivu = 1. / zllm
81         cdivh = 1./zllm      PRINT *, 'cdivu = ', cdivu
     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  
82    
83      END DO      PRINT *, 'Calcul des valeurs propres de nxgrarot'
84        call random_number(zu)
85        zu = zu - 0.5
86        CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
87        call random_number(zv)
88        zv = zv - 0.5
89        CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
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. + (1. - preff / presnivs)**2)      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
103        ! (between 1 and 2)
104    
105      tetaudiv = zvert / tetagdiv      tetaudiv = zvert / tetagdiv
106      tetaurot = zvert / tetagrot      tetaurot = zvert / tetagrot
107      tetah = zvert / tetatemp      tetah = zvert / tetatemp
108    
109        max_zvert = maxval(zvert)
110        tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
111             tetatemp / max_zvert)
112        PRINT *, 'tetamin = ', tetamin
113        idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
114        PRINT *, 'idissip = ', idissip
115        dtdiss = idissip * dtvr
116        PRINT *, 'dtdiss = ', dtdiss, "s"
117    
118      call new_unit(unit)      call new_unit(unit)
119      open(unit, file="inidissip.csv", status="replace", action="write")      open(unit, file="inidissip.csv", status="replace", action="write")
120      write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line  
121        ! Title line:
122        write(unit, fmt=*) '"presnivs (hPa)" "dtdiss * tetaudiv" ' &
123             // '"dtdiss * tetaurot" "dtdiss * tetah"'
124    
125      do l = 1, llm      do l = 1, llm
126         write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)         write(unit, fmt=*) presnivs(l) / 100., dtdiss * tetaudiv(l), &
127                dtdiss * tetaurot(l), dtdiss * tetah(l)
128      end do      end do
129      close(unit)      close(unit)
130      print *, 'Created file "inidissip.csv".'      print *, 'Created file "inidissip.csv".'
131    
     tetamin = min(1E6, minval(1. / tetaudiv), minval(1. / tetaurot), &  
          minval(1. / tetah))  
     PRINT *, 'tetamin = ', tetamin  
     idissip = max(iperiod, int(tetamin / (2 * dtvr * iperiod)) * iperiod)  
     PRINT *, 'idissip = ', idissip  
     dtdiss = idissip * dtvr  
     PRINT *, 'dtdiss = ', dtdiss  
   
132    END SUBROUTINE inidissip    END SUBROUTINE inidissip
133    
134  end module inidissip_m  end module inidissip_m
Legend:
Removed from v.26  
changed lines
  Added in v.65
  ViewVC Help
Powered by ViewVC 1.1.21