LMDZ
disvert.F90
Go to the documentation of this file.
1 ! $Id: disvert.F90 2153 2014-11-24 15:18:11Z lguez $
2 
3 SUBROUTINE disvert()
4 
5 #ifdef CPP_IOIPSL
6  use ioipsl, only: getin
7 #else
8  USE ioipsl_getincom, only: getin
9 #endif
10  use new_unit_m, only: new_unit
11  use assert_m, only: assert
12 
13  IMPLICIT NONE
14 
15  include "dimensions.h"
16  include "paramet.h"
17  include "comvert.h"
18  include "comconst.h"
19  include "iniprint.h"
20  include "logic.h"
21 
22 !-------------------------------------------------------------------------------
23 ! Purpose: Vertical distribution functions for LMDZ.
24 ! Triggered by the levels number llm.
25 !-------------------------------------------------------------------------------
26 ! Read in "comvert.h":
27 
28 ! pa !--- vertical coordinate is close to a PRESSURE COORDINATE FOR P
29 ! < 0.3 * pa (relative variation of p on a model level is < 0.1 %)
30 
31 ! preff !--- REFERENCE PRESSURE (101325 Pa)
32 ! Written in "comvert.h":
33 ! ap(llm+1), bp(llm+1) !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
34 ! aps(llm), bps(llm) !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
35 ! dpres(llm) !--- PRESSURE DIFFERENCE FOR EACH LAYER
36 ! presnivs(llm) !--- PRESSURE AT EACH MID-LAYER
37 ! scaleheight !--- VERTICAL SCALE HEIGHT (Earth: 8kms)
38 ! nivsig(llm+1) !--- SIGMA INDEX OF EACH LAYER INTERFACE
39 ! nivsigs(llm) !--- SIGMA INDEX OF EACH MID-LAYER
40 !-------------------------------------------------------------------------------
41 ! Local variables:
42  REAL sig(llm+1), dsig(llm)
43  REAL sig0(llm+1), zz(llm+1)
44  REAL zk, zkm1, dzk1, dzk2, z, k0, k1
45 
46  INTEGER l, unit
47  REAL dsigmin
48  REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
49 
50  REAL alpha, beta, deltaz
51  REAL x
52  character(len=*),parameter :: modname="disvert"
53 
54  character(len=24):: vert_sampling
55  ! (allowed values are "param", "tropo", "strato" and "read")
56 
57  !-----------------------------------------------------------------------
58 
59  WRITE(lunout,*) trim(modname)//" starts"
60 
61  ! default scaleheight is 8km for earth
62  scaleheight=8.
63 
64  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
65  call getin('vert_sampling', vert_sampling)
66  WRITE(lunout,*) trim(modname)//' vert_sampling = ' // vert_sampling
67  if (llm==39 .and. vert_sampling=="strato") then
68  dsigmin=0.3 ! Vieille option par défaut pour CMIP5
69  else
70  dsigmin=1.
71  endif
72  call getin('dsigmin', dsigmin)
73  WRITE(lunout,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
74 
75 
76  select case (vert_sampling)
77  case ("param")
78  ! On lit les options dans sigma.def:
79  OPEN(99, file='sigma.def', status='old', form='formatted')
80  READ(99, *) scaleheight ! hauteur d'echelle 8.
81  READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
82  READ(99, *) beta ! facteur d'acroissement en haut 1.3
83  READ(99, *) k0 ! nombre de couches dans la transition surf
84  READ(99, *) k1 ! nombre de couches dans la transition haute
85  CLOSE(99)
86  alpha=deltaz/(llm*scaleheight)
87  write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
88  scaleheight, alpha, k0, k1, beta
89 
90  alpha=deltaz/tanh(1./k0)*2.
91  zkm1=0.
92  sig(1)=1.
93  do l=1, llm
94  sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
95  *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
96  *beta**(l-(llm-k1))/log(beta))
97  zk=-scaleheight*log(sig(l+1))
98 
99  dzk1=alpha*tanh(l/k0)
100  dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
101  write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
102  zkm1=zk
103  enddo
104 
105  sig(llm+1)=0.
106 
107  bp(: llm) = exp(1. - 1. / sig(: llm)**2)
108  bp(llmp1) = 0.
109 
110  ap = pa * (sig - bp)
111  case("sigma")
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) = sig(2:llm)
124  bp(llmp1) = 0.
125  ap(:)=0.
126  case("tropo")
127  DO l = 1, llm
128  x = 2*asin(1.) * (l - 0.5) / (llm + 1)
129  dsig(l) = dsigmin + 7.0 * sin(x)**2
130  ENDDO
131  dsig = dsig / sum(dsig)
132  sig(llm+1) = 0.
133  DO l = llm, 1, -1
134  sig(l) = sig(l+1) + dsig(l)
135  ENDDO
136 
137  bp(1)=1.
138  bp(2: llm) = exp(1. - 1. / sig(2: llm)**2)
139  bp(llmp1) = 0.
140 
141  ap(1)=0.
142  ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
143  case("strato")
144  DO l = 1, llm
145  x = 2*asin(1.) * (l - 0.5) / (llm + 1)
146  dsig(l) =(dsigmin + 7. * sin(x)**2) &
147  *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
148  ENDDO
149  dsig = dsig / sum(dsig)
150  sig(llm+1) = 0.
151  DO l = llm, 1, -1
152  sig(l) = sig(l+1) + dsig(l)
153  ENDDO
154 
155  bp(1)=1.
156  bp(2: llm) = exp(1. - 1. / sig(2: llm)**2)
157  bp(llmp1) = 0.
158 
159  ap(1)=0.
160  ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
161  case("strato_correct")
162 !==================================================================
163 ! Fredho 2014/05/18, Saint-Louis du Senegal
164 ! Cette version de la discretisation strato est corrige au niveau
165 ! du passage des sig aux ap, bp
166 ! la version precedente donne un coude dans l'epaisseur des couches
167 ! vers la tropopause
168 !==================================================================
169 
170 
171  DO l = 1, llm
172  x = 2*asin(1.) * (l - 0.5) / (llm + 1)
173  dsig(l) =(dsigmin + 7. * sin(x)**2) &
174  *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
175  ENDDO
176  dsig = dsig / sum(dsig)
177  sig0(llm+1) = 0.
178  DO l = llm, 1, -1
179  sig0(l) = sig0(l+1) + dsig(l)
180  ENDDO
181  sig=racinesig(sig0)
182 
183  bp(1)=1.
184  bp(2:llm)=exp(1.-1./sig(2: llm)**2)
185  bp(llmp1)=0.
186 
187  ap(1)=0.
188  ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
189  ap(llm+1)=0.
190 
191  CASE("strato_custom0")
192 !=======================================================
193 ! Version Transitoire
194  ! custumize strato distribution with specific alpha & beta values and function
195  ! depending on llm (experimental and temporary)!
196  SELECT CASE (llm)
197  CASE(55)
198  alpha=0.45
199  beta=4.0
200  CASE(63)
201  alpha=0.45
202  beta=5.0
203  CASE(71)
204  alpha=3.05
205  beta=65.
206  CASE(79)
207  alpha=3.20
208  ! alpha=2.05 ! FLOTT 79 (PLANTE)
209  beta=70.
210  END SELECT
211  ! Or used values provided by user in def file:
212  CALL getin("strato_alpha",alpha)
213  CALL getin("strato_beta",beta)
214 
215  ! Build geometrical distribution
216  scaleheight=7.
217  zz(1)=0.
218  IF (llm==55.OR.llm==63) THEN
219  DO l=1,llm
220  z=zz(l)/scaleheight
221  zz(l+1)=zz(l)+0.03+z*1.5*(1.-tanh(z-0.5))+alpha*(1.+tanh(z-1.5)) &
222  +5.0*exp((l-llm)/beta)
223  ENDDO
224  ELSEIF (llm==71) THEN !.OR.llm==79) THEN ! FLOTT 79 (PLANTE)
225  DO l=1,llm
226  z=zz(l)
227  zz(l+1)=zz(l)+0.02+0.88*tanh(z/2.5)+alpha*(1.+tanh((z-beta)/15.))
228  ENDDO
229  ELSEIF (llm==79) THEN
230  DO l=1,llm
231  z=zz(l)
232  zz(l+1)=zz(l)+0.02+0.80*tanh(z/3.8)+alpha*(1+tanh((z-beta)/17.)) &
233  +0.03*tanh(z/.25)
234  ENDDO
235  ENDIF ! of IF (llm==55.OR.llm==63) ...
236 
237 
238  ! Build sigma distribution
239  sig0=exp(-zz(:)/scaleheight)
240  sig0(llm+1)=0.
241 ! sig=ridders(sig0)
242  sig=racinesig(sig0)
243 
244  ! Compute ap() and bp()
245  bp(1)=1.
246  bp(2:llm)=exp(1.-1./sig(2:llm)**2)
247  bp(llm+1)=0.
248  ap=pa*(sig-bp)
249 
250  CASE("strato_custom")
251 !===================================================================
252 ! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
253 ! 2014/05
254 ! custumize strato distribution
255 ! Al the parameter are given in km assuming a given scalehigh
256  vert_scale=7. ! scale hight
257  vert_dzmin=0.02 ! width of first layer
258  vert_dzlow=1. ! dz in the low atmosphere
259  vert_z0low=8. ! height at which resolution recches dzlow
260  vert_dzmid=3. ! dz in the mid atmsophere
261  vert_z0mid=70. ! height at which resolution recches dzmid
262  vert_h_mid=20. ! width of the transition
263  vert_dzhig=11. ! dz in the high atmsophere
264  vert_z0hig=80. ! height at which resolution recches dz
265  vert_h_hig=20. ! width of the transition
266 !===================================================================
267 
268  call getin('vert_scale',vert_scale)
269  call getin('vert_dzmin',vert_dzmin)
270  call getin('vert_dzlow',vert_dzlow)
271  call getin('vert_z0low',vert_z0low)
272  CALL getin('vert_dzmid',vert_dzmid)
273  CALL getin('vert_z0mid',vert_z0mid)
274  call getin('vert_h_mid',vert_h_mid)
275  call getin('vert_dzhig',vert_dzhig)
276  call getin('vert_z0hig',vert_z0hig)
277  call getin('vert_h_hig',vert_h_hig)
278 
279  scaleheight=vert_scale ! for consistency with further computations
280  ! Build geometrical distribution
281  zz(1)=0.
282  DO l=1,llm
283  z=zz(l)
284  zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*tanh(z/vert_z0low)+ &
285 & (vert_dzmid-vert_dzlow)* &
286 & (tanh((z-vert_z0mid)/vert_h_mid)-tanh((-vert_z0mid)/vert_h_mid)) &
287 & +(vert_dzhig-vert_dzmid-vert_dzlow)* &
288 & (tanh((z-vert_z0hig)/vert_h_hig)-tanh((-vert_z0hig)/vert_h_hig))
289  ENDDO
290 
291 
292 !===================================================================
293 ! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
294 ! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
295 ! sig0 is p/p0
296 ! sig is an intermediate distribution introduce to estimate bp
297 ! 1. sig0=exp(-z/H)
298 ! 2. inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
299 ! 3. bp=exp(1-1/sig**2)
300 ! 4. ap deduced from the combination of 2 and 3 so that sig0=ap/p0+bp
301 !===================================================================
302 
303  sig0=exp(-zz(:)/vert_scale)
304  sig0(llm+1)=0.
305  sig=racinesig(sig0)
306  bp(1)=1.
307  bp(2:llm)=exp(1.-1./sig(2:llm)**2)
308  bp(llm+1)=0.
309  ap=pa*(sig-bp)
310 
311  case("read")
312  ! Read "ap" and "bp". First line is skipped (title line). "ap"
313  ! should be in Pa. First couple of values should correspond to
314  ! the surface, that is : "bp" should be in descending order.
315  call new_unit(unit)
316  open(unit, file="hybrid.txt", status="old", action="read", &
317  position="rewind")
318  read(unit, fmt=*) ! skip title line
319  do l = 1, llm + 1
320  read(unit, fmt=*) ap(l), bp(l)
321  end do
322  close(unit)
323  call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
324  bp(llm + 1) == 0., "disvert: bad ap or bp values")
325  case default
326  call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
327  END select
328 
329  DO l=1, llm
330  nivsigs(l) = REAL(l)
331  ENDDO
332 
333  DO l=1, llmp1
334  nivsig(l)= REAL(l)
335  ENDDO
336 
337  write(lunout, *) trim(modname),': BP '
338  write(lunout, *) bp
339  write(lunout, *) trim(modname),': AP '
340  write(lunout, *) ap
341 
342  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
343  write(lunout, *)'couches calcules pour une pression de surface =', preff
344  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
345  write(lunout, *) scaleheight,' km'
346  DO l = 1, llm
347  dpres(l) = bp(l) - bp(l+1)
348  presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
349  write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
350  log(preff/presnivs(l))*scaleheight &
351  , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
352  max(ap(l+1)+bp(l+1)*preff, 1.e-10))
353  ENDDO
354 
355  write(lunout, *) trim(modname),': PRESNIVS '
356  write(lunout, *) presnivs
357 
358 CONTAINS
359 
360 !-------------------------------------------------------------------------------
361 !
362 FUNCTION ridders(sig) RESULT(sg)
363 !
364 !-------------------------------------------------------------------------------
365  IMPLICIT NONE
366 !-------------------------------------------------------------------------------
367 ! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
368 ! Notes: Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
369 ! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
370 ! Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
371 !-------------------------------------------------------------------------------
372 ! Arguments:
373  REAL, INTENT(IN) :: sig(:)
374  REAL :: sg(size(sig))
375 !-------------------------------------------------------------------------------
376 ! Local variables:
377  INTEGER :: it, ns, maxit
378  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
379 !-------------------------------------------------------------------------------
380  ns=SIZE(sig); maxit=9999
381  c1=pa/preff; c2=1.-c1
382  DO l=1,ns
383  xx=huge(1.)
384  x1=0.0; f1=distrib(x1,c1,c2,sig(l))
385  x2=1.0; f2=distrib(x2,c1,c2,sig(l))
386  DO it=1,maxit
387  x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
388  s=sqrt(f3**2-f1*f2); IF(s==0.) EXIT
389  x4=x3+(x3-x1)*(sign(1.,f1-f2)*f3/s); IF(abs(10.*log(x4-xx))<=1e-5) EXIT
390  xx=x4; f4=distrib(x4,c1,c2,sig(l)); IF(f4==0.) EXIT
391  IF(sign(f3,f4)/=f3) then; x1=x3; f1=f3; x2=xx; f2=f4
392  ELSE IF(sign(f1,f4)/=f1) then; x2=xx; f2=f4
393  ELSE IF(sign(f2,f4)/=f2) then; x1=xx; f1=f4
394  else; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
395  END IF
396  IF(abs(10.*log(abs(x2-x1)))<=1e-5) EXIT !--- ERROR ON SIG <= 0.01m
397  END DO
398  IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
399  &e for level ',l
400  sg(l)=xx
401  END DO
402  sg(1)=1.; sg(ns)=0.
403 
404 END FUNCTION ridders
405 
406 FUNCTION racinesig(sig) RESULT(sg)
407 !
408 !-------------------------------------------------------------------------------
409  IMPLICIT NONE
410 !-------------------------------------------------------------------------------
411 ! Fredho 2014/05/18
412 ! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
413 ! Notes: Uses Newton Raphson search
414 !-------------------------------------------------------------------------------
415 ! Arguments:
416  REAL, INTENT(IN) :: sig(:)
417  REAL :: sg(size(sig))
418 !-------------------------------------------------------------------------------
419 ! Local variables:
420  INTEGER :: it, ns, maxit
421  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
422 !-------------------------------------------------------------------------------
423  ns=SIZE(sig); maxit=100
424  c1=pa/preff; c2=1.-c1
425  DO l=2,ns-1
426  sg(l)=sig(l)
427  DO it=1,maxit
428  f1=exp(1-1./sg(l)**2)*(1.-c1)
429  sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
430  ENDDO
431 ! print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
432  ENDDO
433  sg(1)=1.; sg(ns)=0.
434 
435 END FUNCTION racinesig
436 
437 
438 
439 
440 END SUBROUTINE disvert
441 
442 !-------------------------------------------------------------------------------
443 
444 FUNCTION distrib(x,c1,c2,x0) RESULT(res)
445 !
446 !-------------------------------------------------------------------------------
447 ! Arguments:
448  REAL, INTENT(IN) :: x, c1, c2, x0
449  REAL :: res
450 !-------------------------------------------------------------------------------
451  res=c1*x+c2*exp(1-1/(x**2))-x0
452 
453 END FUNCTION distrib
454 
455 
ridders
real function, dimension(size(sig)) ridders(sig)
Definition: disvert.F90:363
distrib
real function distrib(x, c1, c2, x0)
Definition: disvert.F90:445
preff
!$Id preff
Definition: comvert.h:8
llmp1
!$Header llmp1
Definition: paramet.h:14
bp
!$Id bp(llm+1)
abort_gcm
subroutine abort_gcm(modname, message, ierr)
Definition: abort_gcm.F:7
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
presnivs
!$Id presnivs(llm)
assert_m::assert
Definition: assert_m.F90:6
ok_strato
!$Id ysinus ok_strato
Definition: logic.h:10
new_unit_m
Definition: new_unit_m.F90:1
nivsigs
!$Id nivsigs(llm)
ioipsl_getincom
Definition: ioipsl_getincom.F90:6
dpres
!$Id dpres(llm)
ioipsl_getincom::getin
Definition: ioipsl_getincom.F90:22
assert_m
Definition: assert_m.F90:2
disvert
subroutine disvert()
Definition: disvert.F90:4
pa
!$Id && pa
Definition: comvert.h:8
scaleheight
!$Id scaleheight
Definition: comvert.h:9
new_unit_m::new_unit
subroutine new_unit(unit)
Definition: new_unit_m.F90:8
racinesig
real function, dimension(size(sig)) racinesig(sig)
Definition: disvert.F90:407
nivsig
!$Id nivsig(llm+1)
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7