source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90 @ 5139

Last change on this file since 5139 was 5134, checked in by abarral, 4 months ago

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