source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dmem/disvert.F90 @ 2136

Last change on this file since 2136 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer


End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

File size: 5.3 KB
Line 
1! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
2
3SUBROUTINE 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
180END SUBROUTINE disvert
Note: See TracBrowser for help on using the repository browser.