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

Last change on this file since 5441 was 5159, checked in by abarral, 5 months ago

Put dimensions.h and paramet.h into modules
  • 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 5159 2024-08-02 19:58:25Z fhourdin $
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
14USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
15  USE lmdz_paramet
16  IMPLICIT NONE
17
18
19
20
21!-------------------------------------------------------------------------------
22! Purpose: Vertical distribution functions for LMDZ.
23!          Triggered by the levels number llm.
24!-------------------------------------------------------------------------------
25! Read    in "comvert_mod":
26
27! pa !--- vertical coordinate is close to a PRESSURE COORDINATE FOR P
28! < 0.3 * pa (relative variation of p on a model level is < 0.1 %)
29
30! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
31! Written in "comvert_mod":
32! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
33! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
34! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
35! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
36! presinter(llm+1)           !--- PRESSURE AT EACH INTERFACE
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  DO l=1, llmp1
358     presinter(l)= ( ap(l)+bp(l)*preff)
359     WRITE(lunout, *)'PRESINTER(', l, ')=', presinter(l)
360  ENDDO
361
362  WRITE(lunout, *) trim(modname),': PRESNIVS '
363  WRITE(lunout, *) presnivs
364
365CONTAINS
366
367!-------------------------------------------------------------------------------
368
369FUNCTION ridders(sig) RESULT(sg)
370
371!-------------------------------------------------------------------------------
372  IMPLICIT NONE
373!-------------------------------------------------------------------------------
374! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
375! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
376! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
377!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
378!-------------------------------------------------------------------------------
379! Arguments:
380  REAL, INTENT(IN)  :: sig(:)
381  REAL              :: sg(SIZE(sig))
382!-------------------------------------------------------------------------------
383! Local variables:
384  INTEGER :: it, ns, maxit
385  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
386!-------------------------------------------------------------------------------
387  ns=SIZE(sig); maxit=9999
388  c1=Pa/Preff; c2=1.-c1
389  DO l=1,ns
390    xx=HUGE(1.)
391    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
392    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
393    DO it=1,maxit
394      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
395      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
396      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
397      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
398      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
399      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
400      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
401      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...', 1)
402      END IF
403      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m           
404    END DO
405    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
406  e for level ',l
407    sg(l)=xx
408  END DO
409  sg(1)=1.; sg(ns)=0.
410
411END FUNCTION ridders
412
413FUNCTION racinesig(sig) RESULT(sg)
414
415!-------------------------------------------------------------------------------
416  IMPLICIT NONE
417!-------------------------------------------------------------------------------
418! Fredho 2014/05/18
419! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
420! Notes:   Uses Newton Raphson search
421!-------------------------------------------------------------------------------
422! Arguments:
423  REAL, INTENT(IN)  :: sig(:)
424  REAL              :: sg(SIZE(sig))
425!-------------------------------------------------------------------------------
426! Local variables:
427  INTEGER :: it, ns, maxit
428  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
429!-------------------------------------------------------------------------------
430  ns=SIZE(sig); maxit=100
431  c1=Pa/Preff; c2=1.-c1
432  DO l=2,ns-1
433    sg(l)=sig(l)
434    DO it=1,maxit
435       f1=exp(1-1./sg(l)**2)*(1.-c1)
436       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
437    ENDDO
438!   PRINT*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
439  ENDDO
440  sg(1)=1.; sg(ns)=0.
441
442END FUNCTION racinesig
443
444
445
446
447END SUBROUTINE disvert
448
449!-------------------------------------------------------------------------------
450
451FUNCTION distrib(x,c1,c2,x0) RESULT(res)
452
453!-------------------------------------------------------------------------------
454! Arguments:
455  REAL, INTENT(IN) :: x, c1, c2, x0
456  REAL             :: res
457!-------------------------------------------------------------------------------
458  res=c1*x+c2*EXP(1-1/(x**2))-x0
459
460END FUNCTION distrib
461
462
Note: See TracBrowser for help on using the repository browser.