Context Navigation

disvert.F90 @ 5227

Last change on this file since 5227 was 2786, checked in by Ehouarn Millour, 8 years ago

Further convergence with planetary models:

move vertical_layers_mod.F90 to phy_common and call init_vertical_layers in inigeomphy_mod
add aps and bps (mid-layer coordinates) to vertical_layers_mod.F90
same adaptations for the 1D case

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.0 KB

Line
1	! $Id: disvert.F90 2786 2017-01-30 10:24:47Z abarral $
2
3	SUBROUTINE 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
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	! 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))*(-alphak0/scaleheight) &
95	exp(-alpha/scaleheighttanh((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=alphatanh((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 = 2asin(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 = 2asin(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 = 2asin(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 = 2asin(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+z1.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.88TANH(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.80TANH(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 ~ ', scaleheightlog((ap(l)+bp(l)preff)/ &
355	max(ap(l+1)+bp(l+1)*preff, 1.e-10))
356	ENDDO
357
358	write(lunout, *) trim(modname),': PRESNIVS '
359	write(lunout, *) presnivs
360
361	CONTAINS
362
363	!-------------------------------------------------------------------------------
364	!
365	FUNCTION ridders(sig) RESULT(sg)
366	!
367	!-------------------------------------------------------------------------------
368	IMPLICIT NONE
369	!-------------------------------------------------------------------------------
370	! Purpose: Search for s solving (Pa/Preff)s+(1-Pa/Preff)EXP(1-1./s**2)=sg
371	! Notes: Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
372	! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
373	! Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
374	!-------------------------------------------------------------------------------
375	! Arguments:
376	REAL, INTENT(IN) :: sig(:)
377	REAL :: sg(SIZE(sig))
378	!-------------------------------------------------------------------------------
379	! Local variables:
380	INTEGER :: it, ns, maxit
381	REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
382	!-------------------------------------------------------------------------------
383	ns=SIZE(sig); maxit=9999
384	c1=Pa/Preff; c2=1.-c1
385	DO l=1,ns
386	xx=HUGE(1.)
387	x1=0.0; f1=distrib(x1,c1,c2,sig(l))
388	x2=1.0; f2=distrib(x2,c1,c2,sig(l))
389	DO it=1,maxit
390	x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
391	s=SQRT(f3*2-f1f2); IF(s==0.) EXIT
392	x4=x3+(x3-x1)(SIGN(1.,f1-f2)f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
393	xx=x4; f4=distrib(x4,c1,c2,sig(l)); IF(f4==0.) EXIT
394	IF(SIGN(f3,f4)/=f3) THEN; x1=x3; f1=f3; x2=xx; f2=f4
395	ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
396	ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
397	ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
398	END IF
399	IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT !--- ERROR ON SIG <= 0.01m
400	END DO
401	IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
402	&e for level ',l
403	sg(l)=xx
404	END DO
405	sg(1)=1.; sg(ns)=0.
406
407	END FUNCTION ridders
408
409	FUNCTION racinesig(sig) RESULT(sg)
410	!
411	!-------------------------------------------------------------------------------
412	IMPLICIT NONE
413	!-------------------------------------------------------------------------------
414	! Fredho 2014/05/18
415	! Purpose: Search for s solving (Pa/Preff)sg+(1-Pa/Preff)EXP(1-1./sg**2)=s
416	! Notes: Uses Newton Raphson search
417	!-------------------------------------------------------------------------------
418	! Arguments:
419	REAL, INTENT(IN) :: sig(:)
420	REAL :: sg(SIZE(sig))
421	!-------------------------------------------------------------------------------
422	! Local variables:
423	INTEGER :: it, ns, maxit
424	REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
425	!-------------------------------------------------------------------------------
426	ns=SIZE(sig); maxit=100
427	c1=Pa/Preff; c2=1.-c1
428	DO l=2,ns-1
429	sg(l)=sig(l)
430	DO it=1,maxit
431	f1=exp(1-1./sg(l)*2)(1.-c1)
432	sg(l)=sg(l)-(c1sg(l)+f1-sig(l))/(c1+2f1sg(l)*(-3))
433	ENDDO
434	! print,'SSSSIG ',sig(l),sg(l),c1sg(l)+exp(1-1./sg(l)*2)(1.-c1)
435	ENDDO
436	sg(1)=1.; sg(ns)=0.
437
438	END FUNCTION racinesig
439
440
441
442
443	END SUBROUTINE disvert
444
445	!-------------------------------------------------------------------------------
446
447	FUNCTION distrib(x,c1,c2,x0) RESULT(res)
448	!
449	!-------------------------------------------------------------------------------
450	! Arguments:
451	REAL, INTENT(IN) :: x, c1, c2, x0
452	REAL :: res
453	!-------------------------------------------------------------------------------
454	res=c1x+c2EXP(1-1/(x**2))-x0
455
456	END FUNCTION distrib
457
458

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format