source: LMDZ6/trunk/libf/dyn3d_common/disvert.f90 @ 5408

Last change on this file since 5408 was 5285, checked in by abarral, 7 weeks ago

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