My Project
 All Classes Files Functions Variables Macros
disvert.F90
Go to the documentation of this file.
1 ! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
2 
3 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
4 
5  ! Auteur : P. Le Van
6 
7  use new_unit_m, only: new_unit
8  use ioipsl, only: getin
9  use assert_m, only: assert
10 
11  IMPLICIT NONE
12 
13  include "dimensions.h"
14  include "paramet.h"
15  include "iniprint.h"
16  include "logic.h"
17 
18  ! s = sigma ** kappa : coordonnee verticale
19  ! dsig(l) : epaisseur de la couche l ds la coord. s
20  ! sig(l) : sigma a l'interface des couches l et l-1
21  ! ds(l) : distance entre les couches l et l-1 en coord.s
22 
23  real,intent(in) :: pa, preff
24  real,intent(out) :: ap(llmp1) ! in Pa
25  real, intent(out):: bp(llmp1)
26  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
27  real,intent(out) :: presnivs(llm)
28  real,intent(out) :: scaleheight
29 
30  REAL sig(llm+1), dsig(llm)
31  real zk, zkm1, dzk1, dzk2, k0, k1
32 
33  INTEGER l, unit
34  REAL dsigmin
35  REAL alpha, beta, deltaz
36  REAL x
37  character(len=*),parameter :: modname="disvert"
38 
39  character(len=6):: vert_sampling
40  ! (allowed values are "param", "tropo", "strato" and "read")
41 
42  !-----------------------------------------------------------------------
43 
44  print *, "Call sequence information: disvert"
45 
46  ! default scaleheight is 8km for earth
47  scaleheight=8.
48 
49  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
50  call getin('vert_sampling', vert_sampling)
51  print *, 'vert_sampling = ' // vert_sampling
52 
53  select case (vert_sampling)
54  case ("param")
55  ! On lit les options dans sigma.def:
56  OPEN(99, file='sigma.def', status='old', form='formatted')
57  READ(99, *) scaleheight ! hauteur d'echelle 8.
58  READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
59  READ(99, *) beta ! facteur d'acroissement en haut 1.3
60  READ(99, *) k0 ! nombre de couches dans la transition surf
61  READ(99, *) k1 ! nombre de couches dans la transition haute
62  CLOSE(99)
63  alpha=deltaz/(llm*scaleheight)
64  write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
65  scaleheight, alpha, k0, k1, beta
66 
67  alpha=deltaz/tanh(1./k0)*2.
68  zkm1=0.
69  sig(1)=1.
70  do l=1, llm
71  sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
72  *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
73  *beta**(l-(llm-k1))/log(beta))
74  zk=-scaleheight*log(sig(l+1))
75 
76  dzk1=alpha*tanh(l/k0)
77  dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
78  write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
79  zkm1=zk
80  enddo
81 
82  sig(llm+1)=0.
83 
84  bp(: llm) = exp(1. - 1. / sig(: llm)**2)
85  bp(llmp1) = 0.
86 
87  ap = pa * (sig - bp)
88  case("tropo")
89  DO l = 1, llm
90  x = 2*asin(1.) * (l - 0.5) / (llm + 1)
91  dsig(l) = 1.0 + 7.0 * sin(x)**2
92  ENDDO
93  dsig = dsig / sum(dsig)
94  sig(llm+1) = 0.
95  DO l = llm, 1, -1
96  sig(l) = sig(l+1) + dsig(l)
97  ENDDO
98 
99  bp(1)=1.
100  bp(2: llm) = exp(1. - 1. / sig(2: llm)**2)
101  bp(llmp1) = 0.
102 
103  ap(1)=0.
104  ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
105  case("strato")
106  if (llm==39) then
107  dsigmin=0.3
108  else if (llm==50) then
109  dsigmin=1.
110  else
111  write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
112  dsigmin=1.
113  endif
114  WRITE(lunout,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
115 
116  DO l = 1, llm
117  x = 2*asin(1.) * (l - 0.5) / (llm + 1)
118  dsig(l) =(dsigmin + 7. * sin(x)**2) &
119  *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
120  ENDDO
121  dsig = dsig / sum(dsig)
122  sig(llm+1) = 0.
123  DO l = llm, 1, -1
124  sig(l) = sig(l+1) + dsig(l)
125  ENDDO
126 
127  bp(1)=1.
128  bp(2: llm) = exp(1. - 1. / sig(2: llm)**2)
129  bp(llmp1) = 0.
130 
131  ap(1)=0.
132  ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
133  case("read")
134  ! Read "ap" and "bp". First line is skipped (title line). "ap"
135  ! should be in Pa. First couple of values should correspond to
136  ! the surface, that is : "bp" should be in descending order.
137  call new_unit(unit)
138  open(unit, file="hybrid.txt", status="old", action="read", &
139  position="rewind")
140  read(unit, fmt=*) ! skip title line
141  do l = 1, llm + 1
142  read(unit, fmt=*) ap(l), bp(l)
143  end do
144  close(unit)
145  call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
146  bp(llm + 1) == 0., "disvert: bad ap or bp values")
147  case default
148  call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
149  END select
150 
151  DO l=1, llm
152  nivsigs(l) = REAL(l)
153  ENDDO
154 
155  DO l=1, llmp1
156  nivsig(l)= REAL(l)
157  ENDDO
158 
159  write(lunout, *) trim(modname),': BP '
160  write(lunout, *) bp
161  write(lunout, *) trim(modname),': AP '
162  write(lunout, *) ap
163 
164  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
165  write(lunout, *)'couches calcules pour une pression de surface =', preff
166  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
167  write(lunout, *) scaleheight,' km'
168  DO l = 1, llm
169  dpres(l) = bp(l) - bp(l+1)
170  presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
171  write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
172  log(preff/presnivs(l))*scaleheight &
173  , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
174  max(ap(l+1)+bp(l+1)*preff, 1.e-10))
175  ENDDO
176 
177  write(lunout, *) trim(modname),': PRESNIVS '
178  write(lunout, *) presnivs
179 
180 END SUBROUTINE disvert