source: LMDZ5/trunk/libf/dyn3d_common/disvert.F90 @ 1992

Last change on this file since 1992 was 1959, checked in by fhourdin, 11 years ago

Nouvelles options dans disvert.F90 :
=------------------------------------

  1. avec vert_sampling="sigma" -> niveaux sigma plutôt qu'hybride.

P_i,j,k= Bhyb_k * Ps , Ahyb_k=0.

  1. disgmin, l'épaisseur approximative du premier niveau en coordonnée

sigma peut être défini dans les fichiers .def.
Cette option est valable pour vert_sampling=sigma, tropo ou strato.
Valeur par défaut : 1, sauf pour (llm==39 .and. vert_sampling=="strato")

ce qui correspond aux configurations CMIP5

=======================================================================
New options in disvert.F90 :
=---------------------------

  1. with vert_samplig="sigma" -> sigma levels instead of hybrid

P_i,j,k= Bhyb_k * Ps , Ahyb_k=0.

  1. disgmin, the approximate thickness of the first layer in sigma

coordinate can be specified from the .def files.
Active for vert_sampling=sigma, tropo or strato.
Defaut value : 1, except for (llm==39 .and. vert_sampling=="strato")

which corresponds to the CMIP5 configurations.

--Cette ligne, et les suivantes ci-dessous, seront ignorées-

M dyn3d_common/disvert.F90

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.6 KB
Line 
1! $Id: disvert.F90 1959 2014-02-07 09:37:38Z 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  if (llm==39 .and. vert_sampling=="strato") then
53     dsigmin=0.3 ! Vieille option par défaut pour CMIP5
54  else
55     dsigmin=1.
56  endif
57  call getin('dsigmin', dsigmin)
58  WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
59
60
61  select case (vert_sampling)
62  case ("param")
63     ! On lit les options dans sigma.def:
64     OPEN(99, file='sigma.def', status='old', form='formatted')
65     READ(99, *) scaleheight ! hauteur d'echelle 8.
66     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
67     READ(99, *) beta ! facteur d'acroissement en haut 1.3
68     READ(99, *) k0 ! nombre de couches dans la transition surf
69     READ(99, *) k1 ! nombre de couches dans la transition haute
70     CLOSE(99)
71     alpha=deltaz/(llm*scaleheight)
72     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
73                               scaleheight, alpha, k0, k1, beta
74
75     alpha=deltaz/tanh(1./k0)*2.
76     zkm1=0.
77     sig(1)=1.
78     do l=1, llm
79        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
80             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
81                  *beta**(l-(llm-k1))/log(beta))
82        zk=-scaleheight*log(sig(l+1))
83
84        dzk1=alpha*tanh(l/k0)
85        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
86        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
87        zkm1=zk
88     enddo
89
90     sig(llm+1)=0.
91
92     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
93     bp(llmp1) = 0.
94
95     ap = pa * (sig - bp)
96  case("sigma")
97     DO l = 1, llm
98        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
99        dsig(l) = dsigmin + 7.0 * SIN(x)**2
100     ENDDO
101     dsig = dsig / sum(dsig)
102     sig(llm+1) = 0.
103     DO l = llm, 1, -1
104        sig(l) = sig(l+1) + dsig(l)
105     ENDDO
106
107     bp(1)=1.
108     bp(2: llm) = sig(2:llm)
109     bp(llmp1) = 0.
110     ap(:)=0.
111  case("tropo")
112     DO l = 1, llm
113        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
114        dsig(l) = dsigmin + 7.0 * SIN(x)**2
115     ENDDO
116     dsig = dsig / sum(dsig)
117     sig(llm+1) = 0.
118     DO l = llm, 1, -1
119        sig(l) = sig(l+1) + dsig(l)
120     ENDDO
121
122     bp(1)=1.
123     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
124     bp(llmp1) = 0.
125
126     ap(1)=0.
127     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
128  case("strato")
129     DO l = 1, llm
130        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
131        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
132             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
133     ENDDO
134     dsig = dsig / sum(dsig)
135     sig(llm+1) = 0.
136     DO l = llm, 1, -1
137        sig(l) = sig(l+1) + dsig(l)
138     ENDDO
139
140     bp(1)=1.
141     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
142     bp(llmp1) = 0.
143
144     ap(1)=0.
145     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
146  case("read")
147     ! Read "ap" and "bp". First line is skipped (title line). "ap"
148     ! should be in Pa. First couple of values should correspond to
149     ! the surface, that is : "bp" should be in descending order.
150     call new_unit(unit)
151     open(unit, file="hybrid.txt", status="old", action="read", &
152          position="rewind")
153     read(unit, fmt=*) ! skip title line
154     do l = 1, llm + 1
155        read(unit, fmt=*) ap(l), bp(l)
156     end do
157     close(unit)
158     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
159          bp(llm + 1) == 0., "disvert: bad ap or bp values")
160  case default
161     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
162  END select
163
164  DO l=1, llm
165     nivsigs(l) = REAL(l)
166  ENDDO
167
168  DO l=1, llmp1
169     nivsig(l)= REAL(l)
170  ENDDO
171
172  write(lunout, *)  trim(modname),': BP '
173  write(lunout, *) bp
174  write(lunout, *)  trim(modname),': AP '
175  write(lunout, *) ap
176
177  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
178  write(lunout, *)'couches calcules pour une pression de surface =', preff
179  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
180  write(lunout, *) scaleheight,' km'
181  DO l = 1, llm
182     dpres(l) = bp(l) - bp(l+1)
183     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
184     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
185          log(preff/presnivs(l))*scaleheight &
186          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
187          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
188  ENDDO
189
190  write(lunout, *) trim(modname),': PRESNIVS '
191  write(lunout, *) presnivs
192
193END SUBROUTINE disvert
Note: See TracBrowser for help on using the repository browser.