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