1 | ! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ |
---|
2 | |
---|
3 | SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) |
---|
4 | |
---|
5 | ! Auteur : P. Le Van |
---|
6 | |
---|
7 | IMPLICIT NONE |
---|
8 | |
---|
9 | include "dimensions.h" |
---|
10 | include "paramet.h" |
---|
11 | include "iniprint.h" |
---|
12 | include "logic.h" |
---|
13 | |
---|
14 | ! s = sigma ** kappa : coordonnee verticale |
---|
15 | ! dsig(l) : epaisseur de la couche l ds la coord. s |
---|
16 | ! sig(l) : sigma a l'interface des couches l et l-1 |
---|
17 | ! ds(l) : distance entre les couches l et l-1 en coord.s |
---|
18 | |
---|
19 | REAL pa, preff |
---|
20 | REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) |
---|
21 | REAL presnivs(llm) |
---|
22 | |
---|
23 | REAL sig(llm+1), dsig(llm) |
---|
24 | real zk, zkm1, dzk1, dzk2, k0, k1 |
---|
25 | |
---|
26 | INTEGER l |
---|
27 | REAL dsigmin |
---|
28 | REAL alpha, beta, deltaz, h |
---|
29 | INTEGER iostat |
---|
30 | REAL pi, x |
---|
31 | |
---|
32 | !----------------------------------------------------------------------- |
---|
33 | |
---|
34 | pi = 2 * ASIN(1.) |
---|
35 | |
---|
36 | OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) |
---|
37 | |
---|
38 | IF (iostat == 0) THEN |
---|
39 | ! cas 1 on lit les options dans sigma.def: |
---|
40 | READ(99, *) h ! hauteur d'echelle 8. |
---|
41 | READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 |
---|
42 | READ(99, *) beta ! facteur d'acroissement en haut 1.3 |
---|
43 | READ(99, *) k0 ! nombre de couches dans la transition surf |
---|
44 | READ(99, *) k1 ! nombre de couches dans la transition haute |
---|
45 | CLOSE(99) |
---|
46 | alpha=deltaz/(llm*h) |
---|
47 | write(lunout, *)'disvert: h, alpha, k0, k1, beta' |
---|
48 | |
---|
49 | alpha=deltaz/tanh(1./k0)*2. |
---|
50 | zkm1=0. |
---|
51 | sig(1)=1. |
---|
52 | do l=1, llm |
---|
53 | sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & |
---|
54 | *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) |
---|
55 | zk=-h*log(sig(l+1)) |
---|
56 | |
---|
57 | dzk1=alpha*tanh(l/k0) |
---|
58 | dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) |
---|
59 | write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 |
---|
60 | zkm1=zk |
---|
61 | enddo |
---|
62 | |
---|
63 | sig(llm+1)=0. |
---|
64 | |
---|
65 | DO l = 1, llm |
---|
66 | dsig(l) = sig(l)-sig(l+1) |
---|
67 | end DO |
---|
68 | ELSE |
---|
69 | if (ok_strato) then |
---|
70 | if (llm==39) then |
---|
71 | dsigmin=0.3 |
---|
72 | else if (llm==50) then |
---|
73 | dsigmin=1. |
---|
74 | else |
---|
75 | WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster' |
---|
76 | dsigmin=1. |
---|
77 | endif |
---|
78 | WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin |
---|
79 | endif |
---|
80 | |
---|
81 | DO l = 1, llm |
---|
82 | x = 2*asin(1.) * (l - 0.5) / (llm + 1) |
---|
83 | |
---|
84 | IF (ok_strato) THEN |
---|
85 | dsig(l) =(dsigmin + 7. * SIN(x)**2) & |
---|
86 | *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 |
---|
87 | ELSE |
---|
88 | dsig(l) = 1.0 + 7.0 * SIN(x)**2 |
---|
89 | ENDIF |
---|
90 | ENDDO |
---|
91 | dsig = dsig / sum(dsig) |
---|
92 | sig(llm+1) = 0. |
---|
93 | DO l = llm, 1, -1 |
---|
94 | sig(l) = sig(l+1) + dsig(l) |
---|
95 | ENDDO |
---|
96 | ENDIF |
---|
97 | |
---|
98 | DO l=1, llm |
---|
99 | nivsigs(l) = REAL(l) |
---|
100 | ENDDO |
---|
101 | |
---|
102 | DO l=1, llmp1 |
---|
103 | nivsig(l)= REAL(l) |
---|
104 | ENDDO |
---|
105 | |
---|
106 | ! .... Calculs de ap(l) et de bp(l) .... |
---|
107 | ! ..... pa et preff sont lus sur les fichiers start par lectba ..... |
---|
108 | |
---|
109 | bp(llmp1) = 0. |
---|
110 | |
---|
111 | DO l = 1, llm |
---|
112 | bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) |
---|
113 | ap(l) = pa * ( sig(l) - bp(l) ) |
---|
114 | ENDDO |
---|
115 | |
---|
116 | bp(1)=1. |
---|
117 | ap(1)=0. |
---|
118 | |
---|
119 | ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) |
---|
120 | |
---|
121 | write(lunout, *)'disvert: BP ' |
---|
122 | write(lunout, *) bp |
---|
123 | write(lunout, *)'disvert: AP ' |
---|
124 | write(lunout, *) ap |
---|
125 | |
---|
126 | write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' |
---|
127 | write(lunout, *)'couches calcules pour une pression de surface =', preff |
---|
128 | write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' |
---|
129 | write(lunout, *)'8km' |
---|
130 | DO l = 1, llm |
---|
131 | dpres(l) = bp(l) - bp(l+1) |
---|
132 | presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) |
---|
133 | write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & |
---|
134 | log(preff/presnivs(l))*8. & |
---|
135 | , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & |
---|
136 | max(ap(l+1)+bp(l+1)*preff, 1.e-10)) |
---|
137 | ENDDO |
---|
138 | |
---|
139 | write(lunout, *) 'disvert: PRESNIVS ' |
---|
140 | write(lunout, *) presnivs |
---|
141 | |
---|
142 | END SUBROUTINE disvert |
---|