source: LMDZ5/branches/testing/libf/dyn3d_common/disvert.F90 @ 3209

Last change on this file since 3209 was 2839, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2785:2838 into testing branch
  • 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: 15.0 KB
Line 
1! $Id: disvert.F90 2839 2017-03-30 14:16:38Z jyg $
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  USE comvert_mod, ONLY: ap, bp, aps, bps, nivsigs, nivsig, dpres, presnivs, &
13                         pseudoalt, pa, preff, scaleheight
14  USE logic_mod, ONLY: ok_strato
15
16  IMPLICIT NONE
17
18  include "dimensions.h"
19  include "paramet.h"
20  include "iniprint.h"
21
22!-------------------------------------------------------------------------------
23! Purpose: Vertical distribution functions for LMDZ.
24!          Triggered by the levels number llm.
25!-------------------------------------------------------------------------------
26! Read    in "comvert_mod":
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_mod":
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     aps(l) =  0.5 *( ap(l) +ap(l+1))
349     bps(l) =  0.5 *( bp(l) +bp(l+1))
350     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
351     pseudoalt(l) = log(preff/presnivs(l))*scaleheight
352     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
353          pseudoalt(l) &
354          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
355          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
356  ENDDO
357
358  write(lunout, *) trim(modname),': PRESNIVS '
359  write(lunout, *) presnivs
360
361CONTAINS
362
363!-------------------------------------------------------------------------------
364!
365FUNCTION ridders(sig) RESULT(sg)
366!
367!-------------------------------------------------------------------------------
368  IMPLICIT NONE
369!-------------------------------------------------------------------------------
370! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
371! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
372! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
373!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
374!-------------------------------------------------------------------------------
375! Arguments:
376  REAL, INTENT(IN)  :: sig(:)
377  REAL              :: sg(SIZE(sig))
378!-------------------------------------------------------------------------------
379! Local variables:
380  INTEGER :: it, ns, maxit
381  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
382!-------------------------------------------------------------------------------
383  ns=SIZE(sig); maxit=9999
384  c1=Pa/Preff; c2=1.-c1
385  DO l=1,ns
386    xx=HUGE(1.)
387    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
388    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
389    DO it=1,maxit
390      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
391      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
392      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
393      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
394      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
395      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
396      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
397      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
398      END IF
399      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m           
400    END DO
401    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
402     &e for level ',l
403    sg(l)=xx
404  END DO
405  sg(1)=1.; sg(ns)=0.
406
407END FUNCTION ridders
408
409FUNCTION racinesig(sig) RESULT(sg)
410!
411!-------------------------------------------------------------------------------
412  IMPLICIT NONE
413!-------------------------------------------------------------------------------
414! Fredho 2014/05/18
415! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
416! Notes:   Uses Newton Raphson search
417!-------------------------------------------------------------------------------
418! Arguments:
419  REAL, INTENT(IN)  :: sig(:)
420  REAL              :: sg(SIZE(sig))
421!-------------------------------------------------------------------------------
422! Local variables:
423  INTEGER :: it, ns, maxit
424  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
425!-------------------------------------------------------------------------------
426  ns=SIZE(sig); maxit=100
427  c1=Pa/Preff; c2=1.-c1
428  DO l=2,ns-1
429    sg(l)=sig(l)
430    DO it=1,maxit
431       f1=exp(1-1./sg(l)**2)*(1.-c1)
432       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
433    ENDDO
434!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
435  ENDDO
436  sg(1)=1.; sg(ns)=0.
437
438END FUNCTION racinesig
439
440
441
442
443END SUBROUTINE disvert
444
445!-------------------------------------------------------------------------------
446
447FUNCTION distrib(x,c1,c2,x0) RESULT(res)
448!
449!-------------------------------------------------------------------------------
450! Arguments:
451  REAL, INTENT(IN) :: x, c1, c2, x0
452  REAL             :: res
453!-------------------------------------------------------------------------------
454  res=c1*x+c2*EXP(1-1/(x**2))-x0
455
456END FUNCTION distrib
457
458
Note: See TracBrowser for help on using the repository browser.