source: LMDZ6/branches/cirrus/libf/dyn3d_common/disvert.F90 @ 5220

Last change on this file since 5220 was 4257, checked in by lguez, 2 years ago

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