source: LMDZ.3.3/trunk/libf/dyn3d/inidissip.F @ 507

Last change on this file since 507 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 5.7 KB
Line 
1      SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
2     *                       tetagdiv,tetagrot,tetatemp             )
3c=======================================================================
4c   initialisation de la dissipation horizontale
5c=======================================================================
6c-----------------------------------------------------------------------
7c   declarations:
8c   -------------
9
10      IMPLICIT NONE
11#include "dimensions.h"
12#include "paramet.h"
13#include "comdissipn.h"
14#include "comconst.h"
15#include "comvert.h"
16#include "control.h"
17
18      LOGICAL lstardis
19      INTEGER nitergdiv,nitergrot,niterh
20      REAL    tetagdiv,tetagrot,tetatemp
21      REAL fact,zvert(llm),zz
22      REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
23      REAL ullm,vllm,umin,vmin,zhmin,zhmax
24      REAL zllm,z1llm
25
26      INTEGER l,ij,idum,ii
27      REAL tetamin
28
29      EXTERNAL ran1
30      REAL ran1
31
32
33c-----------------------------------------------------------------------
34c
35c   calcul des valeurs propres des operateurs par methode iterrative:
36c   -----------------------------------------------------------------
37
38      crot     = -1.
39      cdivu    = -1.
40      cdivh    = -1.
41
42c   calcul de la valeur propre de divgrad:
43c   --------------------------------------
44      idum = 0
45      DO l = 1, llm
46       DO ij = 1, ip1jmp1
47        deltap(ij,l) = 1.
48       ENDDO
49      ENDDO
50
51      idum  = -1
52      zh(1) = RAN1(idum)-.5
53      idum  = 0
54      DO ij = 2, ip1jmp1
55        zh(ij) = RAN1(idum) -.5
56      ENDDO
57
58      CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
59
60      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
61
62      IF ( zhmin .GE. zhmax  )     THEN
63         PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
64         STOP'probleme generateur alleatoire dans inidissip'
65      ENDIF
66
67      zllm = ABS( zhmax )
68      DO l = 1,50
69         IF(lstardis) THEN
70            CALL divgrad2(1,zh,deltap,niterh,zh)
71         ELSE
72            CALL divgrad (1,zh,niterh,zh)
73         ENDIF
74
75        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
76
77         zllm  = ABS( zhmax )
78         z1llm = 1./zllm
79         DO ij = 1,ip1jmp1
80            zh(ij) = zh(ij)* z1llm
81         ENDDO
82      ENDDO
83
84      IF(lstardis) THEN
85         cdivh = 1./ zllm
86      ELSE
87         cdivh = zllm ** ( -1./niterh )
88      ENDIF
89
90c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
91c   -----------------------------------------------------------------
92      print*,'calcul des valeurs propres'
93
94      DO  20  ii = 1, 2
95c
96         DO ij = 1, ip1jmp1
97           zu(ij)  = RAN1(idum) -.5
98         ENDDO
99         CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
100         DO ij = 1, ip1jm
101            zv(ij) = RAN1(idum) -.5
102         ENDDO
103         CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
104
105         CALL minmax(iip1*jjp1,zu,umin,ullm )
106         CALL minmax(iip1*jjm, zv,vmin,vllm )
107
108         ullm = ABS ( ullm )
109         vllm = ABS ( vllm )
110
111         DO  5  l = 1, 50
112            IF(ii.EQ.1) THEN
113ccccc             CALL covcont( 1,zu,zv,zu,zv )
114               IF(lstardis) THEN
115                  CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
116               ELSE
117                  CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
118               ENDIF
119            ELSE
120               IF(lstardis) THEN
121                  CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
122               ELSE
123                  CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
124               ENDIF
125            ENDIF
126
127            CALL minmax(iip1*jjp1,zu,umin,ullm )
128            CALL minmax(iip1*jjm, zv,vmin,vllm )
129
130            ullm = ABS  ( ullm )
131            vllm = ABS  ( vllm )
132
133            zllm  = MAX( ullm,vllm )
134            z1llm = 1./ zllm
135            DO ij = 1, ip1jmp1
136              zu(ij) = zu(ij)* z1llm
137            ENDDO
138            DO ij = 1, ip1jm
139               zv(ij) = zv(ij)* z1llm
140            ENDDO
141 5       CONTINUE
142
143         IF ( ii.EQ.1 ) THEN
144            IF(lstardis) THEN
145               cdivu  = 1./zllm
146            ELSE
147               cdivu  = zllm **( -1./nitergdiv )
148            ENDIF
149         ELSE
150            IF(lstardis) THEN
151               crot   = 1./ zllm
152            ELSE
153               crot   = zllm **( -1./nitergrot )
154            ENDIF
155         ENDIF
156
157 20   CONTINUE
158
159c   petit test pour les operateurs non star:
160c   ----------------------------------------
161
162c     IF(.NOT.lstardis) THEN
163         fact    = rad*24./float(jjm)
164         fact    = fact*fact
165         PRINT*,'coef u ', fact/cdivu, 1./cdivu
166         PRINT*,'coef r ', fact/crot , 1./crot
167         PRINT*,'coef h ', fact/cdivh, 1./cdivh
168c     ENDIF
169
170c-----------------------------------------------------------------------
171c   variation verticale du coefficient de dissipation:
172c   --------------------------------------------------
173
174      DO l=1,llm
175         zvert(l)=1.
176      ENDDO
177
178      fact=2.
179c
180      DO l = 1, llm
181         zz      = 1. - preff/presnivs(l)
182         zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
183      ENDDO
184
185
186      PRINT*,'Constantes de temps de la diffusion horizontale'
187
188      tetamin =  1.e+6
189
190      DO l=1,llm
191        tetaudiv(l)   = zvert(l)/tetagdiv
192        tetaurot(l)   = zvert(l)/tetagrot
193        tetah(l)      = zvert(l)/tetatemp
194
195        IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
196        IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
197        IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
198      ENDDO
199
200      PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
201      idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
202      PRINT *,' INIDI tetamin idissip ',tetamin,idissip
203      idissip = MAX(iperiod,idissip)
204      dtdiss  = idissip * dtvr
205      PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
206
207      DO l = 1,llm
208         PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
209     *                   dtdiss*tetah(l)
210      ENDDO
211
212c
213      RETURN
214      END
Note: See TracBrowser for help on using the repository browser.