/[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/libf/dyn3d/inidissip.f90 revision 54 by guez, Tue Dec 6 15:07:04 2011 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
10  c   declarations:    integer idissip ! période de la dissipation (en pas de temps)
11  c   -------------    real tetaudiv(llm), tetaurot(llm), tetah(llm)
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: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
25        REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)           tetagrot, tetatemp
26        REAL ullm,vllm,umin,vmin,zhmin,zhmax      USE comvert, ONLY : preff, presnivs
27        REAL zllm,z1llm      USE conf_gcm_m, ONLY : iperiod
28        USE dimens_m, ONLY : iim, jjm, llm
29        INTEGER l,ij,idum,ii      USE paramet_m, ONLY : jjp1
30        REAL tetamin      use jumble, only: new_unit
31        use filtreg_m, only: filtreg
32        REAL ran1      use gradiv2_m, only: gradiv2
33    
34        ! Variables local to the procedure:
35  c-----------------------------------------------------------------------      REAL zvert(llm), max_zvert
36        REAL, dimension(iim + 1, jjm + 1):: zh, zu
37        print *, "Call sequence information: inidissip"      real zv(iim + 1, jjm), deltap(iim + 1, jjm + 1, llm)
38  c      REAL zllm
39  c   calcul des valeurs propres des operateurs par methode iterrative:      INTEGER l, seed_size, ii, unit
40  c   -----------------------------------------------------------------      REAL tetamin ! in s
41    
42        crot     = -1.      !-----------------------------------------------------------------------
43        cdivu    = -1.  
44        cdivh    = -1.      PRINT *, 'Call sequence information: inidissip'
45        call random_seed(size=seed_size)
46  c   calcul de la valeur propre de divgrad:      call random_seed(put=(/(0, ii = 1, seed_size)/))
47  c   --------------------------------------  
48        idum = 0      PRINT *, 'Calcul des valeurs propres de divgrad'
49        DO l = 1, llm      deltap = 1.
50         DO ij = 1, ip1jmp1      call random_number(zh)
51          deltap(ij,l) = 1.      zh = zh - 0.5
52         ENDDO      CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
53        ENDDO  
54        DO l = 1, 50
55        idum  = -1         IF (lstardis) THEN
56        zh(1) = RAN1(idum)-.5            CALL divgrad2(1, zh, deltap, niterh, zh, -1.)
57        idum  = 0         ELSE
58        DO ij = 2, ip1jmp1            CALL divgrad(1, zh, niterh, zh, -1.)
59          zh(ij) = RAN1(idum) -.5         END IF
60        ENDDO  
61           zllm = abs(maxval(zh))
62        CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)         zh = zh / zllm
63        END DO
64        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
65        IF (lstardis) THEN
66        IF ( zhmin .GE. zhmax  )     THEN         cdivh = 1. / zllm
67           PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax      ELSE
68           STOP'probleme generateur alleatoire dans inidissip'         cdivh = zllm**(- 1. / niterh)
69        ENDIF      END IF
70        PRINT *, 'cdivh = ', cdivh
71        zllm = ABS( zhmax )  
72        DO l = 1,50      PRINT *, 'Calcul des valeurs propres de gradiv'
73           IF(lstardis) THEN      call random_number(zu)
74              CALL divgrad2(1,zh,deltap,niterh,zh)      zu = zu - 0.5
75           ELSE      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
76              CALL divgrad (1,zh,niterh,zh)      call random_number(zv)
77           ENDIF      zv = zv - 0.5
78        CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
79          CALL minmax(iip1*jjp1,zh,zhmin,zhmax )  
80        DO l = 1, 50
81           zllm  = ABS( zhmax )         IF (lstardis) THEN
82           z1llm = 1./zllm            CALL gradiv2(1, zu, zv, nitergdiv, zu, zv, -1.)
83           DO ij = 1,ip1jmp1         ELSE
84              zh(ij) = zh(ij)* z1llm            CALL gradiv(1, zu, zv, nitergdiv, zu, zv, -1.)
85           ENDDO         END IF
86        ENDDO  
87           zllm = max(abs(maxval(zu)), abs(maxval(zv)))
88        IF(lstardis) THEN         zu = zu / zllm
89           cdivh = 1./ zllm         zv = zv / zllm
90        ELSE      end DO
91           cdivh = zllm ** ( -1./niterh )  
92        ENDIF      IF (lstardis) THEN
93           cdivu = 1. / zllm
94  c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)      ELSE
95  c   -----------------------------------------------------------------         cdivu = zllm**(- 1. / nitergdiv)
96        print*,'calcul des valeurs propres'      END IF
97        PRINT *, 'cdivu = ', cdivu
98        DO  20  ii = 1, 2  
99  c      PRINT *, 'Calcul des valeurs propres de nxgrarot'
100           DO ij = 1, ip1jmp1      call random_number(zu)
101             zu(ij)  = RAN1(idum) -.5      zu = zu - 0.5
102           ENDDO      CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
103           CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)      call random_number(zv)
104           DO ij = 1, ip1jm      zv = zv - 0.5
105              zv(ij) = RAN1(idum) -.5      CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
106           ENDDO  
107           CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)      DO l = 1, 50
108           IF (lstardis) THEN
109           CALL minmax(iip1*jjp1,zu,umin,ullm )            CALL nxgraro2(1, zu, zv, nitergrot, zu, zv, -1.)
110           CALL minmax(iip1*jjm, zv,vmin,vllm )         ELSE
111              CALL nxgrarot(1, zu, zv, nitergrot, zu, zv, -1.)
112           ullm = ABS ( ullm )         END IF
113           vllm = ABS ( vllm )  
114           zllm = max(abs(maxval(zu)), abs(maxval(zv)))
115           DO  5  l = 1, 50         zu = zu / zllm
116              IF(ii.EQ.1) THEN         zv = zv / zllm
117  ccccc             CALL covcont( 1,zu,zv,zu,zv )      end DO
118                 IF(lstardis) THEN  
119                    CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )      IF (lstardis) THEN
120                 ELSE         crot = 1. / zllm
121                    CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )      ELSE
122                 ENDIF         crot = zllm**(-1. / nitergrot)
123              ELSE      END IF
124                 IF(lstardis) THEN      PRINT *, 'crot = ', crot
125                    CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )  
126                 ELSE      ! Variation verticale du coefficient de dissipation :
127                    CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )      zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
128                 ENDIF      ! (between 1 and 2)
129              ENDIF  
130        tetaudiv = zvert / tetagdiv
131              CALL minmax(iip1*jjp1,zu,umin,ullm )      tetaurot = zvert / tetagrot
132              CALL minmax(iip1*jjm, zv,vmin,vllm )      tetah = zvert / tetatemp
133    
134              ullm = ABS  ( ullm )      call new_unit(unit)
135              vllm = ABS  ( vllm )      open(unit, file="inidissip.csv", status="replace", action="write")
136        write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
137              zllm  = MAX( ullm,vllm )      do l = 1, llm
138              z1llm = 1./ zllm         write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
139              DO ij = 1, ip1jmp1      end do
140                zu(ij) = zu(ij)* z1llm      close(unit)
141              ENDDO      print *, 'Created file "inidissip.csv".'
142              DO ij = 1, ip1jm  
143                 zv(ij) = zv(ij)* z1llm      max_zvert = maxval(zvert)
144              ENDDO      tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
145   5       CONTINUE           tetatemp / max_zvert)
146        PRINT *, 'tetamin = ', tetamin
147           IF ( ii.EQ.1 ) THEN      idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
148              IF(lstardis) THEN      PRINT *, 'idissip = ', idissip
149                 cdivu  = 1./zllm      dtdiss = idissip * dtvr
150              ELSE      PRINT *, 'dtdiss = ', dtdiss
151                 cdivu  = zllm **( -1./nitergdiv )  
152              ENDIF    END SUBROUTINE inidissip
153           ELSE  
154              IF(lstardis) THEN  end module inidissip_m
                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.54
  ViewVC Help
Powered by ViewVC 1.1.21