1 | ! |
---|
2 | ! $Id: inidissip.F90 1611 2012-01-25 14:31:54Z idelkadi $ |
---|
3 | ! |
---|
4 | SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & |
---|
5 | tetagdiv,tetagrot,tetatemp ) |
---|
6 | !======================================================================= |
---|
7 | ! initialisation de la dissipation horizontale |
---|
8 | !======================================================================= |
---|
9 | !----------------------------------------------------------------------- |
---|
10 | ! declarations: |
---|
11 | ! ------------- |
---|
12 | |
---|
13 | USE control_mod, only : dissip_period,iperiod |
---|
14 | |
---|
15 | IMPLICIT NONE |
---|
16 | include "dimensions.h" |
---|
17 | include "paramet.h" |
---|
18 | include "comdissipn.h" |
---|
19 | include "comconst.h" |
---|
20 | include "comvert.h" |
---|
21 | include "logic.h" |
---|
22 | include "iniprint.h" |
---|
23 | |
---|
24 | LOGICAL,INTENT(in) :: lstardis |
---|
25 | INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh |
---|
26 | REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp |
---|
27 | |
---|
28 | ! Local variables: |
---|
29 | REAL fact,zvert(llm),zz |
---|
30 | REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) |
---|
31 | real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm) |
---|
32 | REAL ullm,vllm,umin,vmin,zhmin,zhmax |
---|
33 | REAL zllm |
---|
34 | |
---|
35 | INTEGER l,ij,idum,ii |
---|
36 | REAL tetamin |
---|
37 | REAL pseudoz |
---|
38 | character (len=80) :: abort_message |
---|
39 | |
---|
40 | REAL ran1 |
---|
41 | |
---|
42 | |
---|
43 | !----------------------------------------------------------------------- |
---|
44 | ! |
---|
45 | ! calcul des valeurs propres des operateurs par methode iterrative: |
---|
46 | ! ----------------------------------------------------------------- |
---|
47 | |
---|
48 | crot = -1. |
---|
49 | cdivu = -1. |
---|
50 | cdivh = -1. |
---|
51 | |
---|
52 | ! calcul de la valeur propre de divgrad: |
---|
53 | ! -------------------------------------- |
---|
54 | idum = 0 |
---|
55 | DO l = 1, llm |
---|
56 | DO ij = 1, ip1jmp1 |
---|
57 | deltap(ij,l) = 1. |
---|
58 | ENDDO |
---|
59 | ENDDO |
---|
60 | |
---|
61 | idum = -1 |
---|
62 | zh(1) = RAN1(idum)-.5 |
---|
63 | idum = 0 |
---|
64 | DO ij = 2, ip1jmp1 |
---|
65 | zh(ij) = RAN1(idum) -.5 |
---|
66 | ENDDO |
---|
67 | |
---|
68 | CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) |
---|
69 | |
---|
70 | CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) |
---|
71 | |
---|
72 | IF ( zhmin .GE. zhmax ) THEN |
---|
73 | write(lunout,*)' Inidissip zh min max ',zhmin,zhmax |
---|
74 | abort_message='probleme generateur alleatoire dans inidissip' |
---|
75 | call abort_gcm('inidissip',abort_message,1) |
---|
76 | ENDIF |
---|
77 | |
---|
78 | zllm = ABS( zhmax ) |
---|
79 | DO l = 1,50 |
---|
80 | IF(lstardis) THEN |
---|
81 | CALL divgrad2(1,zh,deltap,niterh,divgra) |
---|
82 | ELSE |
---|
83 | CALL divgrad (1,zh,niterh,divgra) |
---|
84 | ENDIF |
---|
85 | |
---|
86 | zllm = ABS(maxval(divgra)) |
---|
87 | zh = divgra / zllm |
---|
88 | ENDDO |
---|
89 | |
---|
90 | IF(lstardis) THEN |
---|
91 | cdivh = 1./ zllm |
---|
92 | ELSE |
---|
93 | cdivh = zllm ** ( -1./niterh ) |
---|
94 | ENDIF |
---|
95 | |
---|
96 | ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) |
---|
97 | ! ----------------------------------------------------------------- |
---|
98 | write(lunout,*)'inidissip: calcul des valeurs propres' |
---|
99 | |
---|
100 | DO ii = 1, 2 |
---|
101 | ! |
---|
102 | DO ij = 1, ip1jmp1 |
---|
103 | zu(ij) = RAN1(idum) -.5 |
---|
104 | ENDDO |
---|
105 | CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) |
---|
106 | DO ij = 1, ip1jm |
---|
107 | zv(ij) = RAN1(idum) -.5 |
---|
108 | ENDDO |
---|
109 | CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) |
---|
110 | |
---|
111 | CALL minmax(iip1*jjp1,zu,umin,ullm ) |
---|
112 | CALL minmax(iip1*jjm, zv,vmin,vllm ) |
---|
113 | |
---|
114 | ullm = ABS ( ullm ) |
---|
115 | vllm = ABS ( vllm ) |
---|
116 | |
---|
117 | DO l = 1, 50 |
---|
118 | IF(ii.EQ.1) THEN |
---|
119 | !cccc CALL covcont( 1,zu,zv,zu,zv ) |
---|
120 | IF(lstardis) THEN |
---|
121 | CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) |
---|
122 | ELSE |
---|
123 | CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) |
---|
124 | ENDIF |
---|
125 | ELSE |
---|
126 | IF(lstardis) THEN |
---|
127 | CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) |
---|
128 | ELSE |
---|
129 | CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) |
---|
130 | ENDIF |
---|
131 | ENDIF |
---|
132 | |
---|
133 | zllm = max(abs(maxval(gx)), abs(maxval(gy))) |
---|
134 | zu = gx / zllm |
---|
135 | zv = gy / zllm |
---|
136 | end DO |
---|
137 | |
---|
138 | IF ( ii.EQ.1 ) THEN |
---|
139 | IF(lstardis) THEN |
---|
140 | cdivu = 1./zllm |
---|
141 | ELSE |
---|
142 | cdivu = zllm **( -1./nitergdiv ) |
---|
143 | ENDIF |
---|
144 | ELSE |
---|
145 | IF(lstardis) THEN |
---|
146 | crot = 1./ zllm |
---|
147 | ELSE |
---|
148 | crot = zllm **( -1./nitergrot ) |
---|
149 | ENDIF |
---|
150 | ENDIF |
---|
151 | |
---|
152 | end DO |
---|
153 | |
---|
154 | ! petit test pour les operateurs non star: |
---|
155 | ! ---------------------------------------- |
---|
156 | |
---|
157 | ! IF(.NOT.lstardis) THEN |
---|
158 | fact = rad*24./REAL(jjm) |
---|
159 | fact = fact*fact |
---|
160 | write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu |
---|
161 | write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot |
---|
162 | write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh |
---|
163 | ! ENDIF |
---|
164 | |
---|
165 | !----------------------------------------------------------------------- |
---|
166 | ! variation verticale du coefficient de dissipation: |
---|
167 | ! -------------------------------------------------- |
---|
168 | |
---|
169 | if (ok_strato .and. llm==39) then |
---|
170 | do l=1,llm |
---|
171 | pseudoz=8.*log(preff/presnivs(l)) |
---|
172 | zvert(l)=1+ & |
---|
173 | (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. & |
---|
174 | *(dissip_factz-1.) |
---|
175 | enddo |
---|
176 | else |
---|
177 | DO l=1,llm |
---|
178 | zvert(l)=1. |
---|
179 | ENDDO |
---|
180 | fact=2. |
---|
181 | DO l = 1, llm |
---|
182 | zz = 1. - preff/presnivs(l) |
---|
183 | zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) |
---|
184 | ENDDO |
---|
185 | endif |
---|
186 | |
---|
187 | |
---|
188 | write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' |
---|
189 | |
---|
190 | tetamin = 1.e+6 |
---|
191 | |
---|
192 | DO l=1,llm |
---|
193 | tetaudiv(l) = zvert(l)/tetagdiv |
---|
194 | tetaurot(l) = zvert(l)/tetagrot |
---|
195 | tetah(l) = zvert(l)/tetatemp |
---|
196 | |
---|
197 | IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) |
---|
198 | IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) |
---|
199 | IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) |
---|
200 | ENDDO |
---|
201 | |
---|
202 | ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def |
---|
203 | IF (dissip_period == 0) THEN |
---|
204 | dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod |
---|
205 | write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period |
---|
206 | dissip_period = MAX(iperiod,dissip_period) |
---|
207 | END IF |
---|
208 | |
---|
209 | dtdiss = dissip_period * dtvr |
---|
210 | write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr |
---|
211 | |
---|
212 | DO l = 1,llm |
---|
213 | write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & |
---|
214 | dtdiss*tetah(l) |
---|
215 | ENDDO |
---|
216 | |
---|
217 | END SUBROUTINE inidissip |
---|