source: LMDZ5/branches/AI-cosp/libf/dyn3d_common/disvert.F90

Last change on this file was 2153, checked in by lguez, 10 years ago

Corrected a comment. See discussion on POIHL mailing list, June 30th
2014.

  • 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: 14.7 KB
Line 
1! $Id: disvert.F90 2153 2014-11-24 15:18:11Z aclsce $
2
3SUBROUTINE 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
358CONTAINS
359
360!-------------------------------------------------------------------------------
361!
362FUNCTION 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
404END FUNCTION ridders
405
406FUNCTION 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
435END FUNCTION racinesig
436
437
438
439
440END SUBROUTINE disvert
441
442!-------------------------------------------------------------------------------
443
444FUNCTION 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
453END FUNCTION distrib
454
455
Note: See TracBrowser for help on using the repository browser.