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