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 |