1 |
|
|
! $Id: $ |
2 |
|
|
SUBROUTINE disvert_noterre |
3 |
|
|
|
4 |
|
|
c Auteur : F. Forget Y. Wanherdrick, P. Levan |
5 |
|
|
c Nouvelle version 100% Mars !! |
6 |
|
|
c On l'utilise aussi pour Venus et Titan, legerment modifiee. |
7 |
|
|
|
8 |
|
|
#ifdef CPP_IOIPSL |
9 |
|
|
use IOIPSL |
10 |
|
|
#else |
11 |
|
|
! if not using IOIPSL, we still need to use (a local version of) getin |
12 |
|
|
use ioipsl_getincom |
13 |
|
|
#endif |
14 |
|
|
USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt, |
15 |
|
|
& nivsig,nivsigs,pa,preff,scaleheight |
16 |
|
|
USE comconst_mod, ONLY: kappa |
17 |
|
|
USE logic_mod, ONLY: hybrid |
18 |
|
|
|
19 |
|
|
IMPLICIT NONE |
20 |
|
|
|
21 |
|
|
include "dimensions.h" |
22 |
|
|
include "paramet.h" |
23 |
|
|
include "iniprint.h" |
24 |
|
|
c |
25 |
|
|
c======================================================================= |
26 |
|
|
c Discretisation verticale en coordonn�e hybride (ou sigma) |
27 |
|
|
c |
28 |
|
|
c======================================================================= |
29 |
|
|
c |
30 |
|
|
c declarations: |
31 |
|
|
c ------------- |
32 |
|
|
c |
33 |
|
|
c |
34 |
|
|
INTEGER l,ll |
35 |
|
|
REAL snorm |
36 |
|
|
REAL alpha,beta,gama,delta,deltaz |
37 |
|
|
real quoi,quand |
38 |
|
|
REAL zsig(llm),sig(llm+1) |
39 |
|
|
INTEGER np,ierr |
40 |
|
|
integer :: ierr1,ierr2,ierr3,ierr4 |
41 |
|
|
REAL x |
42 |
|
|
|
43 |
|
|
REAL SSUM |
44 |
|
|
EXTERNAL SSUM |
45 |
|
|
real newsig |
46 |
|
|
REAL dz0,dz1,nhaut,sig1,esig,csig,zz |
47 |
|
|
real tt,rr,gg, prevz |
48 |
|
|
real s(llm),dsig(llm) |
49 |
|
|
|
50 |
|
|
integer iz |
51 |
|
|
real z, ps,p |
52 |
|
|
character(len=*),parameter :: modname="disvert_noterre" |
53 |
|
|
|
54 |
|
|
c |
55 |
|
|
c----------------------------------------------------------------------- |
56 |
|
|
c |
57 |
|
|
! Initializations: |
58 |
|
|
! pi=2.*ASIN(1.) ! already done in iniconst |
59 |
|
|
|
60 |
|
|
hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates) |
61 |
|
|
CALL getin('hybrid',hybrid) |
62 |
|
|
write(lunout,*) trim(modname),': hybrid=',hybrid |
63 |
|
|
|
64 |
|
|
! Ouverture possible de fichiers typiquement E.T. |
65 |
|
|
|
66 |
|
|
open(99,file="esasig.def",status='old',form='formatted', |
67 |
|
|
s iostat=ierr2) |
68 |
|
|
if(ierr2.ne.0) then |
69 |
|
|
close(99) |
70 |
|
|
open(99,file="z2sig.def",status='old',form='formatted', |
71 |
|
|
s iostat=ierr4) |
72 |
|
|
endif |
73 |
|
|
|
74 |
|
|
c----------------------------------------------------------------------- |
75 |
|
|
c cas 1 on lit les options dans esasig.def: |
76 |
|
|
c ---------------------------------------- |
77 |
|
|
|
78 |
|
|
IF(ierr2.eq.0) then |
79 |
|
|
|
80 |
|
|
c Lecture de esasig.def : |
81 |
|
|
c Systeme peu souple, mais qui respecte en theorie |
82 |
|
|
c La conservation de l'energie (conversion Energie potentielle |
83 |
|
|
c <-> energie cinetique, d'apres la note de Frederic Hourdin... |
84 |
|
|
|
85 |
|
|
write(lunout,*)'*****************************' |
86 |
|
|
write(lunout,*)'WARNING reading esasig.def' |
87 |
|
|
write(lunout,*)'*****************************' |
88 |
|
|
READ(99,*) scaleheight |
89 |
|
|
READ(99,*) dz0 |
90 |
|
|
READ(99,*) dz1 |
91 |
|
|
READ(99,*) nhaut |
92 |
|
|
CLOSE(99) |
93 |
|
|
|
94 |
|
|
dz0=dz0/scaleheight |
95 |
|
|
dz1=dz1/scaleheight |
96 |
|
|
|
97 |
|
|
sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut) |
98 |
|
|
|
99 |
|
|
esig=1. |
100 |
|
|
|
101 |
|
|
do l=1,20 |
102 |
|
|
esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.) |
103 |
|
|
enddo |
104 |
|
|
csig=(1./sig1-1.)/(exp(esig)-1.) |
105 |
|
|
|
106 |
|
|
DO L = 2, llm |
107 |
|
|
zz=csig*(exp(esig*(l-1.))-1.) |
108 |
|
|
sig(l) =1./(1.+zz) |
109 |
|
|
& * tanh(.5*(llm+1-l)/nhaut) |
110 |
|
|
ENDDO |
111 |
|
|
sig(1)=1. |
112 |
|
|
sig(llm+1)=0. |
113 |
|
|
quoi = 1. + 2.* kappa |
114 |
|
|
s( llm ) = 1. |
115 |
|
|
s(llm-1) = quoi |
116 |
|
|
IF( llm.gt.2 ) THEN |
117 |
|
|
DO ll = 2, llm-1 |
118 |
|
|
l = llm+1 - ll |
119 |
|
|
quand = sig(l+1)/ sig(l) |
120 |
|
|
s(l-1) = quoi * (1.-quand) * s(l) + quand * s(l+1) |
121 |
|
|
ENDDO |
122 |
|
|
END IF |
123 |
|
|
c |
124 |
|
|
snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2) |
125 |
|
|
DO l = 1, llm |
126 |
|
|
s(l) = s(l)/ snorm |
127 |
|
|
ENDDO |
128 |
|
|
|
129 |
|
|
c----------------------------------------------------------------------- |
130 |
|
|
c cas 2 on lit les options dans z2sig.def: |
131 |
|
|
c ---------------------------------------- |
132 |
|
|
|
133 |
|
|
ELSE IF(ierr4.eq.0) then |
134 |
|
|
write(lunout,*)'****************************' |
135 |
|
|
write(lunout,*)'Reading z2sig.def' |
136 |
|
|
write(lunout,*)'****************************' |
137 |
|
|
|
138 |
|
|
READ(99,*) scaleheight |
139 |
|
|
do l=1,llm |
140 |
|
|
read(99,*) zsig(l) |
141 |
|
|
end do |
142 |
|
|
CLOSE(99) |
143 |
|
|
|
144 |
|
|
sig(1) =1 |
145 |
|
|
do l=2,llm |
146 |
|
|
sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + |
147 |
|
|
& exp(-zsig(l-1)/scaleheight) ) |
148 |
|
|
end do |
149 |
|
|
sig(llm+1) =0 |
150 |
|
|
|
151 |
|
|
c----------------------------------------------------------------------- |
152 |
|
|
ELSE |
153 |
|
|
write(lunout,*) 'didn t you forget something ??? ' |
154 |
|
|
write(lunout,*) 'We need file z2sig.def ! (OR esasig.def)' |
155 |
|
|
stop |
156 |
|
|
ENDIF |
157 |
|
|
c----------------------------------------------------------------------- |
158 |
|
|
|
159 |
|
|
DO l=1,llm |
160 |
|
|
nivsigs(l) = REAL(l) |
161 |
|
|
ENDDO |
162 |
|
|
|
163 |
|
|
DO l=1,llmp1 |
164 |
|
|
nivsig(l)= REAL(l) |
165 |
|
|
ENDDO |
166 |
|
|
|
167 |
|
|
|
168 |
|
|
c----------------------------------------------------------------------- |
169 |
|
|
c .... Calculs de ap(l) et de bp(l) .... |
170 |
|
|
c ......................................... |
171 |
|
|
c |
172 |
|
|
c ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... |
173 |
|
|
c----------------------------------------------------------------------- |
174 |
|
|
c |
175 |
|
|
|
176 |
|
|
if (hybrid) then ! use hybrid coordinates |
177 |
|
|
write(lunout,*) "*********************************" |
178 |
|
|
write(lunout,*) "Using hybrid vertical coordinates" |
179 |
|
|
write(lunout,*) |
180 |
|
|
c Coordonnees hybrides avec mod |
181 |
|
|
DO l = 1, llm |
182 |
|
|
|
183 |
|
|
call sig_hybrid(sig(l),pa,preff,newsig) |
184 |
|
|
bp(l) = EXP( 1. - 1./(newsig**2) ) |
185 |
|
|
ap(l) = pa * (newsig - bp(l) ) |
186 |
|
|
enddo |
187 |
|
|
bp(llmp1) = 0. |
188 |
|
|
ap(llmp1) = 0. |
189 |
|
|
else ! use sigma coordinates |
190 |
|
|
write(lunout,*) "********************************" |
191 |
|
|
write(lunout,*) "Using sigma vertical coordinates" |
192 |
|
|
write(lunout,*) |
193 |
|
|
c Pour ne pas passer en coordonnees hybrides |
194 |
|
|
DO l = 1, llm |
195 |
|
|
ap(l) = 0. |
196 |
|
|
bp(l) = sig(l) |
197 |
|
|
ENDDO |
198 |
|
|
ap(llmp1) = 0. |
199 |
|
|
endif |
200 |
|
|
|
201 |
|
|
bp(llmp1) = 0. |
202 |
|
|
|
203 |
|
|
write(lunout,*) trim(modname),': BP ' |
204 |
|
|
write(lunout,*) bp |
205 |
|
|
write(lunout,*) trim(modname),': AP ' |
206 |
|
|
write(lunout,*) ap |
207 |
|
|
|
208 |
|
|
c Calcul au milieu des couches : |
209 |
|
|
c WARNING : le choix de placer le milieu des couches au niveau de |
210 |
|
|
c pression interm�diaire est arbitraire et pourrait etre modifi�. |
211 |
|
|
c Le calcul du niveau pour la derniere couche |
212 |
|
|
c (on met la meme distance (en log pression) entre P(llm) |
213 |
|
|
c et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est |
214 |
|
|
c Specifique. Ce choix est sp�cifi� ici ET dans exner_milieu.F |
215 |
|
|
|
216 |
|
|
DO l = 1, llm-1 |
217 |
|
|
aps(l) = 0.5 *( ap(l) +ap(l+1)) |
218 |
|
|
bps(l) = 0.5 *( bp(l) +bp(l+1)) |
219 |
|
|
ENDDO |
220 |
|
|
|
221 |
|
|
if (hybrid) then |
222 |
|
|
aps(llm) = aps(llm-1)**2 / aps(llm-2) |
223 |
|
|
bps(llm) = 0.5*(bp(llm) + bp(llm+1)) |
224 |
|
|
else |
225 |
|
|
bps(llm) = bps(llm-1)**2 / bps(llm-2) |
226 |
|
|
aps(llm) = 0. |
227 |
|
|
end if |
228 |
|
|
|
229 |
|
|
write(lunout,*) trim(modname),': BPs ' |
230 |
|
|
write(lunout,*) bps |
231 |
|
|
write(lunout,*) trim(modname),': APs' |
232 |
|
|
write(lunout,*) aps |
233 |
|
|
|
234 |
|
|
DO l = 1, llm |
235 |
|
|
presnivs(l) = aps(l)+bps(l)*preff |
236 |
|
|
pseudoalt(l) = -scaleheight*log(presnivs(l)/preff) |
237 |
|
|
ENDDO |
238 |
|
|
|
239 |
|
|
write(lunout,*)trim(modname),' : PRESNIVS' |
240 |
|
|
write(lunout,*)presnivs |
241 |
|
|
write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', |
242 |
|
|
& 'height of ',scaleheight,' km)' |
243 |
|
|
write(lunout,*)pseudoalt |
244 |
|
|
|
245 |
|
|
c -------------------------------------------------- |
246 |
|
|
c This can be used to plot the vertical discretization |
247 |
|
|
c (> xmgrace -nxy testhybrid.tab ) |
248 |
|
|
c -------------------------------------------------- |
249 |
|
|
c open (53,file='testhybrid.tab') |
250 |
|
|
c scaleheight=15.5 |
251 |
|
|
c do iz=0,34 |
252 |
|
|
c z = -5 + min(iz,34-iz) |
253 |
|
|
c approximation of scale height for Venus |
254 |
|
|
c scaleheight = 15.5 - z/55.*10. |
255 |
|
|
c ps = preff*exp(-z/scaleheight) |
256 |
|
|
c zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff) |
257 |
|
|
c do l=2,llm |
258 |
|
|
c approximation of scale height for Venus |
259 |
|
|
c if (zsig(l-1).le.55.) then |
260 |
|
|
c scaleheight = 15.5 - zsig(l-1)/55.*10. |
261 |
|
|
c else |
262 |
|
|
c scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2. |
263 |
|
|
c endif |
264 |
|
|
c zsig(l)= zsig(l-1)-scaleheight* |
265 |
|
|
c . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) |
266 |
|
|
c end do |
267 |
|
|
c write(53,'(I3,50F10.5)') iz, zsig |
268 |
|
|
c end do |
269 |
|
|
c close(53) |
270 |
|
|
c -------------------------------------------------- |
271 |
|
|
|
272 |
|
|
|
273 |
|
|
RETURN |
274 |
|
|
END |
275 |
|
|
|
276 |
|
|
c ************************************************************ |
277 |
|
|
subroutine sig_hybrid(sig,pa,preff,newsig) |
278 |
|
|
c ---------------------------------------------- |
279 |
|
|
c Subroutine utilisee pour calculer des valeurs de sigma modifie |
280 |
|
|
c pour conserver les coordonnees verticales decrites dans |
281 |
|
|
c esasig.def/z2sig.def lors du passage en coordonnees hybrides |
282 |
|
|
c F. Forget 2002 |
283 |
|
|
c Connaissant sig (niveaux "sigma" ou on veut mettre les couches) |
284 |
|
|
c L'objectif est de calculer newsig telle que |
285 |
|
|
c (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig |
286 |
|
|
c Cela ne se r�soud pas analytiquement: |
287 |
|
|
c => on r�soud par iterration bourrine |
288 |
|
|
c ---------------------------------------------- |
289 |
|
|
c Information : where exp(1-1./x**2) become << x |
290 |
|
|
c x exp(1-1./x**2) /x |
291 |
|
|
c 1 1 |
292 |
|
|
c 0.68 0.5 |
293 |
|
|
c 0.5 1.E-1 |
294 |
|
|
c 0.391 1.E-2 |
295 |
|
|
c 0.333 1.E-3 |
296 |
|
|
c 0.295 1.E-4 |
297 |
|
|
c 0.269 1.E-5 |
298 |
|
|
c 0.248 1.E-6 |
299 |
|
|
c => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25 |
300 |
|
|
|
301 |
|
|
|
302 |
|
|
implicit none |
303 |
|
|
real x1, x2, sig,pa,preff, newsig, F |
304 |
|
|
integer j |
305 |
|
|
|
306 |
|
|
newsig = sig |
307 |
|
|
x1=0 |
308 |
|
|
x2=1 |
309 |
|
|
if (sig.ge.1) then |
310 |
|
|
newsig= sig |
311 |
|
|
else if (sig*preff/pa.ge.0.25) then |
312 |
|
|
DO J=1,9999 ! nombre d''iteration max |
313 |
|
|
F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig |
314 |
|
|
c write(0,*) J, ' newsig =', newsig, ' F= ', F |
315 |
|
|
if (F.gt.1) then |
316 |
|
|
X2 = newsig |
317 |
|
|
newsig=(X1+newsig)*0.5 |
318 |
|
|
else |
319 |
|
|
X1 = newsig |
320 |
|
|
newsig=(X2+newsig)*0.5 |
321 |
|
|
end if |
322 |
|
|
c Test : on arete lorsque on approxime sig � moins de 0.01 m pr�s |
323 |
|
|
c (en pseudo altitude) : |
324 |
|
|
IF(abs(10.*log(F)).LT.1.E-5) goto 999 |
325 |
|
|
END DO |
326 |
|
|
else ! if (sig*preff/pa.le.0.25) then |
327 |
|
|
newsig= sig*preff/pa |
328 |
|
|
end if |
329 |
|
|
999 continue |
330 |
|
|
Return |
331 |
|
|
END |