source: trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90 @ 3553

Last change on this file since 3553 was 2426, checked in by emillour, 4 years ago

Common dynamics:
Fix a couple of harmless typos spoted by YJ.
EM
File size: 14.9 KB
Line 
1! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
2
3SUBROUTINE disvert()
4
5  ! Auteur : P. Le Van
6
7#ifdef CPP_IOIPSL
8  use ioipsl, only: getin
9#else
10  USE ioipsl_getincom, only: getin
11#endif
12  use new_unit_m, only: new_unit
13  use assert_m, only: assert
14  USE comvert_mod, ONLY: ap,bp,nivsigs,nivsig,preff,pa,presnivs,dpres,scaleheight
15  USE logic_mod, ONLY: ok_strato
16
17  IMPLICIT NONE
18
19  include "dimensions.h"
20  include "paramet.h"
21  include "iniprint.h"
22
23!-------------------------------------------------------------------------------
24! Purpose: Vertical distribution functions for LMDZ.
25!          Triggered by the levels number llm.
26!-------------------------------------------------------------------------------
27! Read    in "comvert_mod":
28
29! pa !--- vertical coordinate is close to a PRESSURE COORDINATE FOR P
30! < 0.3 * pa (relative variation of p on a model level is < 0.1 %)
31
32! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
33! Written in "comvert_mod":
34! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
35! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
36! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
37! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
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     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
350     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
351          log(preff/presnivs(l))*scaleheight &
352          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
353          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
354  ENDDO
355
356  write(lunout, *) trim(modname),': PRESNIVS '
357  write(lunout, *) presnivs
358
359CONTAINS
360
361!-------------------------------------------------------------------------------
362!
363FUNCTION ridders(sig) RESULT(sg)
364!
365!-------------------------------------------------------------------------------
366  IMPLICIT NONE
367!-------------------------------------------------------------------------------
368! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
369! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
370! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
371!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
372!-------------------------------------------------------------------------------
373! Arguments:
374  REAL, INTENT(IN)  :: sig(:)
375  REAL              :: sg(SIZE(sig))
376!-------------------------------------------------------------------------------
377! Local variables:
378  INTEGER :: it, ns, maxit
379  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
380!-------------------------------------------------------------------------------
381  ns=SIZE(sig); maxit=9999
382  c1=Pa/Preff; c2=1.-c1
383  DO l=1,ns
384    xx=HUGE(1.)
385    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
386    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
387    DO it=1,maxit
388      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
389      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
390      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
391      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
392      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
393      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
394      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
395      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...',1)
396      END IF
397      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m           
398    END DO
399    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
400     &e for level ',l
401    sg(l)=xx
402  END DO
403  sg(1)=1.; sg(ns)=0.
404
405END FUNCTION ridders
406
407FUNCTION racinesig(sig) RESULT(sg)
408!
409!-------------------------------------------------------------------------------
410  IMPLICIT NONE
411!-------------------------------------------------------------------------------
412! Fredho 2014/05/18
413! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
414! Notes:   Uses Newton Raphson search
415!-------------------------------------------------------------------------------
416! Arguments:
417  REAL, INTENT(IN)  :: sig(:)
418  REAL              :: sg(SIZE(sig))
419!-------------------------------------------------------------------------------
420! Local variables:
421  INTEGER :: it, ns, maxit
422  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
423!-------------------------------------------------------------------------------
424  ns=SIZE(sig); maxit=100
425  c1=Pa/Preff; c2=1.-c1
426  DO l=2,ns-1
427    sg(l)=sig(l)
428    DO it=1,maxit
429       f1=exp(1-1./sg(l)**2)*(1.-c1)
430       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
431    ENDDO
432!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
433  ENDDO
434  sg(1)=1.; sg(ns)=0.
435
436END FUNCTION racinesig
437
438
439
440
441END SUBROUTINE disvert
442
443!-------------------------------------------------------------------------------
444
445FUNCTION distrib(x,c1,c2,x0) RESULT(res)
446!
447!-------------------------------------------------------------------------------
448! Arguments:
449  REAL, INTENT(IN) :: x, c1, c2, x0
450  REAL             :: res
451!-------------------------------------------------------------------------------
452  res=c1*x+c2*EXP(1-1/(x**2))-x0
453
454END FUNCTION distrib
455
456
457
Note: See TracBrowser for help on using the repository browser.