/[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/dyn3d/Dissipation/inidissip.f revision 82 by guez, Wed Mar 5 14:57:53 2014 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 ! in s
10  c   declarations:    integer idissip ! période de la dissipation (en pas de temps)
11  c   -------------    real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s-1
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"  
20        ! Initialisation de la dissipation horizontale. Calcul des valeurs
21        LOGICAL lstardis      ! propres des opérateurs par méthode itérative.
22        INTEGER nitergdiv,nitergrot,niterh  
23        REAL    tetagdiv,tetagrot,tetatemp      USE comconst, ONLY: dtvr
24        REAL fact,zvert(llm),zz      use comdissnew, only: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
25        REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)           tetatemp
26        REAL ullm,vllm,umin,vmin,zhmin,zhmax      USE disvert_m, ONLY: preff, presnivs
27        REAL zllm,z1llm      USE conf_gcm_m, ONLY: iperiod
28        USE dimens_m, ONLY: iim, jjm
29        INTEGER l,ij,idum,ii      use divgrad2_m, only: divgrad2
30        REAL tetamin      use filtreg_m, only: filtreg
31        use gradiv2_m, only: gradiv2
32        REAL ran1      use jumble, only: new_unit
33        use nxgraro2_m, only: nxgraro2
34        USE paramet_m, ONLY: jjp1
35  c-----------------------------------------------------------------------  
36        ! Variables local to the procedure:
37        print *, "Call sequence information: inidissip"      REAL zvert(llm), max_zvert ! no dimension
38  c      REAL, dimension(iim + 1, jjm + 1, 1):: zh, zu, gx, divgra, deltap
39  c   calcul des valeurs propres des operateurs par methode iterrative:      real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1)
40  c   -----------------------------------------------------------------      REAL zllm
41        INTEGER l, seed_size, ii, unit
42        crot     = -1.      REAL tetamin ! in s
43        cdivu    = -1.  
44        cdivh    = -1.      !-----------------------------------------------------------------------
45    
46  c   calcul de la valeur propre de divgrad:      PRINT *, 'Call sequence information: inidissip'
47  c   --------------------------------------      call random_seed(size=seed_size)
48        idum = 0      call random_seed(put=(/(1, ii = 1, seed_size)/))
49        DO l = 1, llm  
50         DO ij = 1, ip1jmp1      PRINT *, 'Calcul des valeurs propres de divgrad'
51          deltap(ij,l) = 1.      deltap = 1.
52         ENDDO      call random_number(zh)
53        ENDDO      zh = zh - 0.5
54        CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE.)
55        idum  = -1  
56        zh(1) = RAN1(idum)-.5      DO l = 1, 50
57        idum  = 0         CALL divgrad2(1, zh, deltap, niterh, divgra, -1.)
58        DO ij = 2, ip1jmp1         zllm = abs(maxval(divgra))
59          zh(ij) = RAN1(idum) -.5         zh = divgra / zllm
60        ENDDO      END DO
61    
62        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)      cdivh = 1. / zllm
63        PRINT *, 'cdivh = ', cdivh
64        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
65        PRINT *, 'Calcul des valeurs propres de gradiv'
66        IF ( zhmin .GE. zhmax  )     THEN      call random_number(zu)
67           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax      zu = zu - 0.5
68           STOP'probleme generateur alleatoire dans inidissip'      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
69        ENDIF      call random_number(zv)
70        zv = zv - 0.5
71        zllm = ABS( zhmax )      CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
72        DO l = 1,50  
73           IF(lstardis) THEN      DO l = 1, 50
74              CALL divgrad2(1,zh,deltap,niterh,zh)         CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.)
75           ELSE         zllm = max(abs(maxval(gx)), abs(maxval(gy)))
76              CALL divgrad (1,zh,niterh,zh)         zu = gx / zllm
77           ENDIF         zv = gy / zllm
78        end DO
79          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
80        cdivu = 1. / zllm
81           zllm  = ABS( zhmax )      PRINT *, 'cdivu = ', cdivu
82           z1llm = 1./zllm  
83           DO ij = 1,ip1jmp1      PRINT *, 'Calcul des valeurs propres de nxgrarot'
84              zh(ij) = zh(ij)* z1llm      call random_number(zu)
85           ENDDO      zu = zu - 0.5
86        ENDDO      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.)
87        call random_number(zv)
88        IF(lstardis) THEN      zv = zv - 0.5
89           cdivh = 1./ zllm      CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)
90        ELSE  
91           cdivh = zllm ** ( -1./niterh )      DO l = 1, 50
92        ENDIF         CALL nxgraro2(zu, zv, nitergrot, gx, gy, -1.)
93           zllm = max(abs(maxval(gx)), abs(maxval(gy)))
94  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)         zu = gx / zllm
95  c   -----------------------------------------------------------------         zv = gy / zllm
96        print*,'calcul des valeurs propres'      end DO
97    
98        DO  20  ii = 1, 2      crot = 1. / zllm
99  c      PRINT *, 'crot = ', crot
100           DO ij = 1, ip1jmp1  
101             zu(ij)  = RAN1(idum) -.5      ! Variation verticale du coefficient de dissipation :
102           ENDDO      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
103           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)      ! (between 1 and 2)
104           DO ij = 1, ip1jm  
105              zv(ij) = RAN1(idum) -.5      tetaudiv = zvert / tetagdiv
106           ENDDO      tetaurot = zvert / tetagrot
107           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)      tetah = zvert / tetatemp
108    
109           CALL minmax(iip1*jjp1,zu,umin,ullm )      max_zvert = maxval(zvert)
110           CALL minmax(iip1*jjm, zv,vmin,vllm )      tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
111             tetatemp / max_zvert)
112           ullm = ABS ( ullm )      PRINT *, 'tetamin = ', tetamin
113           vllm = ABS ( vllm )      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
114        PRINT *, 'idissip = ', idissip
115           DO  5  l = 1, 50      dtdiss = idissip * dtvr
116              IF(ii.EQ.1) THEN      PRINT *, 'dtdiss = ', dtdiss, "s"
117  ccccc             CALL covcont( 1,zu,zv,zu,zv )  
118                 IF(lstardis) THEN      call new_unit(unit)
119                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )      open(unit, file="inidissip.csv", status="replace", action="write")
120                 ELSE  
121                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )      ! Title line:
122                 ENDIF      write(unit, fmt=*) '"presnivs (hPa)" "dtdiss * tetaudiv" ' &
123              ELSE           // '"dtdiss * tetaurot" "dtdiss * tetah"'
124                 IF(lstardis) THEN  
125                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )      do l = 1, llm
126                 ELSE         write(unit, fmt=*) presnivs(l) / 100., dtdiss * tetaudiv(l), &
127                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )              dtdiss * tetaurot(l), dtdiss * tetah(l)
128                 ENDIF      end do
129              ENDIF      close(unit)
130        print *, 'Created file "inidissip.csv".'
131              CALL minmax(iip1*jjp1,zu,umin,ullm )  
132              CALL minmax(iip1*jjm, zv,vmin,vllm )    END SUBROUTINE inidissip
133    
134              ullm = ABS  ( ullm )  end module inidissip_m
             vllm = ABS  ( vllm )  
   
             zllm  = MAX( ullm,vllm )  
             z1llm = 1./ zllm  
             DO ij = 1, ip1jmp1  
               zu(ij) = zu(ij)* z1llm  
             ENDDO  
             DO ij = 1, ip1jm  
                zv(ij) = zv(ij)* z1llm  
             ENDDO  
  5       CONTINUE  
   
          IF ( ii.EQ.1 ) THEN  
             IF(lstardis) THEN  
                cdivu  = 1./zllm  
             ELSE  
                cdivu  = zllm **( -1./nitergdiv )  
             ENDIF  
          ELSE  
             IF(lstardis) THEN  
                crot   = 1./ zllm  
             ELSE  
                crot   = zllm **( -1./nitergrot )  
             ENDIF  
          ENDIF  
   
  20   CONTINUE  
   
 c   petit test pour les operateurs non star:  
 c   ----------------------------------------  
   
 c     IF(.NOT.lstardis) THEN  
          fact    = rad*24./float(jjm)  
          fact    = fact*fact  
          PRINT*,'coef u ', fact/cdivu, 1./cdivu  
          PRINT*,'coef r ', fact/crot , 1./crot  
          PRINT*,'coef h ', fact/cdivh, 1./cdivh  
 c     ENDIF  
   
 c-----------------------------------------------------------------------  
 c   variation verticale du coefficient de dissipation:  
 c   --------------------------------------------------  
   
       DO l=1,llm  
          zvert(l)=1.  
       ENDDO  
   
       fact=2.  
 c  
       DO l = 1, llm  
          zz      = 1. - preff/presnivs(l)  
          zvert(l)= fact -( fact-1.)/( 1.+zz*zz )  
       ENDDO  
   
   
       PRINT*,'Constantes de temps de la diffusion horizontale'  
   
       tetamin =  1.e+6  
   
       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.82
  ViewVC Help
Powered by ViewVC 1.1.21