/[lmdze]/trunk/libf/dyn3d/inidissip.f
ViewVC logotype

Contents of /trunk/libf/dyn3d/inidissip.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21