Context Navigation

inidissip.F90 @ 1793

Last change on this file since 1793 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 11.5 KB

Line
1	!
2	! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
3	!
4	SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , &
5	tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
6	!=======================================================================
7	! Initialization for horizontal (lateral) dissipation
8	! - in all cases, there is a multiplicative coefficient which increases
9	! the dissipation in the higher levels of the atmosphere, but there
10	! are different ways of seting the vertical profile of this coefficient
11	! (see code below).
12	! - the call to dissipation, every 'dissip_period' dynamical time step,
13	! can be imposed via 'dissip_period=...' in run.def (or via the
14	! 'idissip=...' flag, but this flag should become obsolete, and is
15	! overridden by the 'dissip_period' flag). Note that setting
16	! dissip_period=0 has the special meaning of requesting an "optimal"
17	! value for "dissip_period" (then taken as the largest possible value)
18	! - the three characteristic time scales (relative to the smallest
19	! longitudinal grid mesh size), which are privided in run.def, which
20	! are used for the dissipations steps are:
21	! tetagdiv : time scale for the gradient of the divergence of velocities
22	! tetagrot : time scale for the curl of the curl of velocities
23	! tetatemp : time scale for the laplacian of the potential temperature
24	!=======================================================================
25
26	USE control_mod, only : dissip_period,iperiod,planet_type
27	USE comvert_mod, ONLY: preff,presnivs,scaleheight,pseudoalt
28	USE comconst_mod, ONLY: dtvr,dtdiss,rad,pi,dissip_zref,dissip_deltaz, &
29	dissip_factz,dissip_fac_mid,dissip_fac_up,dissip_pupstart, &
30	dissip_hdelta
31	USE logic_mod, ONLY: ok_strato
32
33	IMPLICIT NONE
34	include "dimensions.h"
35	include "paramet.h"
36	include "comdissipn.h"
37	include "iniprint.h"
38
39	LOGICAL,INTENT(in) :: lstardis
40	INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
41	REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
42
43	integer, INTENT(in):: vert_prof_dissip ! Vertical profile of horizontal dissipation
44	! For the Earth model:
45	! Allowed values:
46	! 0: rational fraction, function of pressure
47	! 1: tanh of altitude
48	! For planets:
49	! 0: use fac_mid (read from run.def)
50	! 1: use fac_mid, fac_up, startalt, delta (hard coded in inidissip.F90)
51	! Local variables:
52	REAL fact,zvert(llm),zz
53	REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
54	real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
55	REAL ullm,vllm,umin,vmin,zhmin,zhmax
56	REAL zllm
57
58	INTEGER l,ij,idum,ii
59	REAL tetamin
60	REAL pseudoz
61	REAL Pup
62	character (len=80) :: abort_message
63
64	REAL ran1
65	logical,save :: firstcall=.true.
66	real,save :: fac_mid,fac_up,startalt,delta,middle
67
68	if (firstcall) then
69	firstcall=.false.
70	if ((planet_type.ne."earth").and.(vert_prof_dissip.eq.1)) then
71	! initialize values for dissipation variation along the vertical (Mars)
72	fac_mid=3 ! coefficient for lower/middle atmosphere
73	fac_up=30 ! coefficient for upper atmosphere
74	startalt=70. ! altitude (in km) for the transition from middle to upper atm.
75	delta=30.! Size (in km) of the transition region between middle/upper atm.
76
77	middle=startalt+delta/2
78	write(lunout,*)"inidissip: multiplicative factors in altitude:", &
79	" fac_mid=",fac_mid," fac_up=",fac_up
80	write(lunout,*)" transition mid/up : startalt (km) =",startalt, &
81	" delta (km) =",delta
82	endif
83	endif !of if firstcall
84
85	!-----------------------------------------------------------------------
86	!
87	! calcul des valeurs propres des operateurs par methode iterrative:
88	! -----------------------------------------------------------------
89
90	crot = -1.
91	cdivu = -1.
92	cdivh = -1.
93
94	! calcul de la valeur propre de divgrad:
95	! --------------------------------------
96	idum = 0
97	DO l = 1, llm
98	DO ij = 1, ip1jmp1
99	deltap(ij,l) = 1.
100	ENDDO
101	ENDDO
102
103	idum = -1
104	zh(1) = RAN1(idum)-.5
105	idum = 0
106	DO ij = 2, ip1jmp1
107	zh(ij) = RAN1(idum) -.5
108	ENDDO
109
110	CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
111
112	CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
113
114	IF ( zhmin .GE. zhmax ) THEN
115	write(lunout,*)' Inidissip zh min max ',zhmin,zhmax
116	abort_message='probleme generateur alleatoire dans inidissip'
117	call abort_gcm('inidissip',abort_message,1)
118	ENDIF
119
120	zllm = ABS( zhmax )
121	DO l = 1,50
122	IF(lstardis) THEN
123	CALL divgrad2(1,zh,deltap,niterh,divgra)
124	ELSE
125	CALL divgrad (1,zh,niterh,divgra)
126	ENDIF
127
128	zllm = ABS(maxval(divgra))
129	zh = divgra / zllm
130	ENDDO
131
132	IF(lstardis) THEN
133	cdivh = 1./ zllm
134	ELSE
135	cdivh = zllm ** ( -1./niterh )
136	ENDIF
137
138	! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2)
139	! -----------------------------------------------------------------
140	write(lunout,*)'inidissip: calcul des valeurs propres'
141
142	DO ii = 1, 2
143	!
144	DO ij = 1, ip1jmp1
145	zu(ij) = RAN1(idum) -.5
146	ENDDO
147	CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
148	DO ij = 1, ip1jm
149	zv(ij) = RAN1(idum) -.5
150	ENDDO
151	CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
152
153	CALL minmax(iip1*jjp1,zu,umin,ullm )
154	CALL minmax(iip1*jjm, zv,vmin,vllm )
155
156	ullm = ABS ( ullm )
157	vllm = ABS ( vllm )
158
159	DO l = 1, 50
160	IF(ii.EQ.1) THEN
161	!cccc CALL covcont( 1,zu,zv,zu,zv )
162	IF(lstardis) THEN
163	CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
164	ELSE
165	CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
166	ENDIF
167	ELSE
168	IF(lstardis) THEN
169	CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
170	ELSE
171	CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
172	ENDIF
173	ENDIF
174
175	zllm = max(abs(maxval(gx)), abs(maxval(gy)))
176	zu = gx / zllm
177	zv = gy / zllm
178	end DO
179
180	IF ( ii.EQ.1 ) THEN
181	IF(lstardis) THEN
182	cdivu = 1./zllm
183	ELSE
184	cdivu = zllm **( -1./nitergdiv )
185	ENDIF
186	ELSE
187	IF(lstardis) THEN
188	crot = 1./ zllm
189	ELSE
190	crot = zllm **( -1./nitergrot )
191	ENDIF
192	ENDIF
193
194	end DO
195
196	! petit test pour les operateurs non star:
197	! ----------------------------------------
198
199	! IF(.NOT.lstardis) THEN
200	fact = rad*24./REAL(jjm)
201	fact = fact*fact
202	write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
203	write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
204	write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
205	! ENDIF
206
207	!-----------------------------------------------------------------------
208	! variation verticale du coefficient de dissipation:
209	! --------------------------------------------------
210
211	if (planet_type.eq."earth") then
212
213	if (vert_prof_dissip == 1) then
214	do l=1,llm
215	pseudoz=8.*log(preff/presnivs(l))
216	zvert(l)=1+ &
217	(tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
218	*(dissip_factz-1.)
219	enddo
220	else
221	DO l=1,llm
222	zvert(l)=1.
223	ENDDO
224	fact=2.
225	DO l = 1, llm
226	zz = 1. - preff/presnivs(l)
227	zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
228	ENDDO
229	endif ! of if (vert_prof_dissip == 1)
230
231	else ! other planets
232
233	if (vert_prof_dissip == 0) then
234	! First step: going from 1 to dissip_fac_mid (in gcm.def)
235	!============
236	DO l=1,llm
237	zz = 1. - preff/presnivs(l)
238	zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
239	ENDDO
240
241	write(lunout,*) 'Dissipation : '
242	write(lunout,*) 'Multiplication de la dissipation en altitude :'
243	write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid
244
245	! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def)
246	!==========================
247	! Utilisation de la fonction tangente hyperbolique pour augmenter
248	! arbitrairement la dissipation et donc la stabilite du modele a
249	! partir d'une certaine altitude.
250
251	! Le facteur multiplicatif de basse atmosphere etant deja pris
252	! en compte, il faut diviser le facteur multiplicatif de haute
253	! atmosphere par celui-ci.
254
255	if (ok_strato) then
256
257	Pup = dissip_pupstartexp(-0.5dissip_deltaz/dissip_hdelta)
258	do l=1,llm
259	zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
260	(1-(0.5(1+tanh(-6./dissip_deltaz* &
261	(-dissip_hdelta*log(presnivs(l)/Pup)) )))) ))
262	enddo
263
264	write(,) ' dissip_fac_up =', dissip_fac_up
265	write(,) 'Transition mid /up: Pupstart,delta =', &
266	dissip_pupstart,'Pa', dissip_deltaz , '(km)'
267
268	endif ! of if (ok_strato)
269	elseif (vert_prof_dissip==1) then
270	DO l=1,llm
271	zz = 1. - preff/presnivs(l)
272	! zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
273	zvert(l)= fac_mid -( fac_mid-1.)/( 1.+zz*zz )
274
275	zvert(l)= zvert(l)(1.0+((fac_up/fac_mid-1) &
276	(1-(0.5*(1+tanh(-6./ &
277	delta(scaleheight(-log(presnivs(l)/preff))-middle))))) &
278	))
279	ENDDO
280	write(lunout,*) "inidissip: vert_prof_disip=1, scaleheight=",scaleheight
281	write(lunout,*) " fac_mid=",fac_mid,", fac_up=",fac_up
282
283	else
284	write(lunout,*) 'wrong value for vert_prof_dissip:',vert_prof_dissip
285	abort_message='wrong value for vert_prof_dissip'
286	call abort_gcm('inidissip',abort_message,1)
287	endif ! of (vert_prof_dissip == 0)
288	endif ! of if (planet_type.eq."earth")
289
290
291	write(lunout,*)'inidissip: Time constants for lateral dissipation'
292
293	tetamin = 1.e+6
294
295	DO l=1,llm
296	tetaudiv(l) = zvert(l)/tetagdiv
297	tetaurot(l) = zvert(l)/tetagrot
298	tetah(l) = zvert(l)/tetatemp
299
300	IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
301	IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
302	IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l)
303	ENDDO
304
305	! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
306	IF (dissip_period == 0) THEN
307	dissip_period = INT( tetamin/( 2.dtvriperiod) ) * iperiod
308	write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
309	dissip_period = MAX(iperiod,dissip_period)
310	END IF
311
312	dtdiss = dissip_period * dtvr
313	write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
314
315	write(lunout,*)'pseudoZ(km) zvert dt(tetagdiv) dt(tetagrot) dt(divgrad)'
316	DO l = 1,llm
317	write(lunout,'(f6.1,x,4(1pe14.7,x))') &
318	pseudoalt(l),zvert(l),dtdisstetaudiv(l),dtdisstetaurot(l),dtdiss*tetah(l)
319	! test if disipation is not too strong (for Explicit Euler time marching)
320	if (dtdiss*tetaudiv(l).gt.1.9) then
321	write(lunout,*)"STOP : lateral dissipation is too intense and will"
322	write(lunout,*)" generate instabilities in the model !"
323	write(lunout,*)" You must increase tetagdiv (or increase dissip_period"
324	write(lunout,*)" or increase day_stap)"
325	endif
326	if (dtdiss*tetaurot(l).gt.1.9) then
327	write(lunout,*)"STOP : lateral dissipation is too intense and will"
328	write(lunout,*)" generate instabilities in the model !"
329	write(lunout,*)" You must increase tetaurot (or increase dissip_period"
330	write(lunout,*)" or increase day_stap)"
331	endif
332	if (dtdiss*tetah(l).gt.1.9) then
333	write(lunout,*)"STOP : lateral dissipation is too intense and will"
334	write(lunout,*)" generate instabilities in the model !"
335	write(lunout,*)" You must increase tetah (or increase dissip_period"
336	write(lunout,*)" or increase day_stap)"
337	endif
338	ENDDO ! of DO l=1,llm
339
340	END SUBROUTINE inidissip

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

Download in other formats:

Original Format