Context Navigation

flott_gwd_ran.F90 @ 780

Last change on this file since 780 was 778, checked in by slebonnois, 12 years ago
SL: bug correction for a file specific to Venus and Titan, + Venus gravity waves routine update (work still in progress)
File size: 12.2 KB

Line
1	SUBROUTINE FLOTT_GWD_RAN(NLON,NLEV,DTIME, pp, pn2, &
2	tt,uu,vv,zustr,zvstr,d_t, d_u, d_v)
3
4	!----------------------------------------------------------------------
5	! Parametrization of the momentum flux deposition due to a discrete
6	! number of gravity waves.
7	! F. Lott (version 9: 16 February, 2012), reproduce v3 but with only
8	! two waves present at each time step
9	! LMDz model online version
10	! ADAPTED FOR VENUS
11	!---------------------------------------------------------------------
12
13	use dimphy
14	implicit none
15
16	#include "dimensions.h"
17	#include "paramet.h"
18
19	#include "YOEGWD.h"
20	#include "YOMCST.h"
21
22	! 0. DECLARATIONS:
23
24	! 0.1 INPUTS
25	INTEGER, intent(in):: NLON, NLEV
26	REAL, intent(in):: DTIME ! Time step of the Physics
27	REAL, intent(in):: pp(NLON, NLEV) ! Pressure at full levels
28	! VENUS ATTENTION: CP VARIABLE PN2 CALCULE EN AMONT DES PARAMETRISATIONS
29	REAL, intent(in):: pn2(NLON,NLEV) ! N2 (BV^2) at 1/2 levels
30	REAL, intent(in):: TT(NLON, NLEV) ! Temp at full levels
31
32	REAL, intent(in):: UU(NLON, NLEV) , VV(NLON, NLEV)
33	! Hor winds at full levels
34
35	! 0.2 OUTPUTS
36	REAL, intent(out):: zustr(NLON), zvstr(NLON) ! Surface Stresses
37	REAL, intent(inout):: d_t(NLON, NLEV) ! Tendency on Temp.
38
39	REAL, intent(inout):: d_u(NLON, NLEV), d_v(NLON, NLEV)
40	! Tendencies on winds
41
42	! O.3 INTERNAL ARRAYS
43
44	INTEGER II, LL, IEQ
45
46	! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
47
48	REAL DELTAT
49
50	! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
51
52	!VENUS INTEGER, PARAMETER:: NK = 4, NP = 4, NO = 4, NW = NK * NP * NO
53	!Online output: change NO
54	INTEGER, PARAMETER:: NK = 1, NP = 2, NO = 10, NW = NK * NP * NO
55	INTEGER JK, JP, JO, JW
56	REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
57	REAL CMIN, CMAX ! Min and Max absolute ph. vel.
58	REAL CPHA ! absolute PHASE VELOCITY frequency
59	REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
60	REAL ZP(NW) ! Horizontal wavenumber angle
61	REAL ZO(NW, KLON) ! Absolute frequency !
62
63	! Waves Intr. freq. at the 1/2 lev surrounding the full level
64	REAL ZOM(NW, KLON), ZOP(NW, KLON)
65
66	! Wave vertical velocities at the 2 1/2 lev surrounding the full level
67	REAL WWM(NW, KLON), WWP(NW, KLON)
68
69	REAL RUW0(NW, KLON) ! Fluxes at launching level
70
71	REAL RUWP(NW, KLON), RVWP(NW, KLON)
72	! Fluxes X and Y for each waves at 1/2 Levels
73
74	INTEGER LAUNCH ! Launching altitude
75
76	REAL RUWMAX,SAT ! saturation parameter
77	REAL XLAUNCH ! Controle the launching altitude
78	REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
79	REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
80
81	! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
82
83	REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
84
85	! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
86
87	REAL H0(KLON, KLEV) ! Characteristic Height of the atmosphere
88	REAL PR ! Reference Pressure
89
90	REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
91
92	REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
93	REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
94	REAL PSEC ! Security to avoid division by 0 pressure
95	REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
96	REAL BVSEC ! Security to avoid negative BVF
97
98	! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
99	logical output
100	data output/.false./
101	! CAUTION ! IF output is .true. THEN change NO to 10 at least !
102	character*14 outform
103	character*2 str2
104
105	! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
106	real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
107	LOGICAL firstcall
108	SAVE firstcall
109	DATA firstcall/.true./
110
111	!-----------------------------------------------------------------
112	! 1. INITIALISATIONS
113
114	IF (firstcall) THEN
115	allocate(d_u_sav(NLON,NLEV),d_v_sav(NLON,NLEV))
116	firstcall=.false.
117	ENDIF
118
119	! 1.1 Basic parameter
120
121	! PARAMETERS CORRESPONDING TO V3:
122	RUWMAX = 0.005 ! Max EP-Flux at Launch altitude
123	SAT = 0.85 ! Saturation parameter: Sc in (12)
124	RDISS = 10. ! Diffusion parameter
125
126	DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)
127
128	KMIN = 1.E-6 ! Min horizontal wavenumber
129	KMAX = 2.E-5 ! Max horizontal wavenumber
130	!Online output: one value only
131	if (output) then
132	KMIN = 1.3E-5
133	KMAX = 1.3E-5
134	endif
135	CMIN = 1. ! Min phase velocity
136	CMAX = 60. ! Max phase speed velocity
137	XLAUNCH=0.6 ! Parameter that control launching altitude
138
139	PR = 9.2e6 ! Reference pressure ! VENUS!!
140
141	BVSEC = 1.E-5 ! Security to avoid negative BVF
142	PSEC = 1.E-6 ! Security to avoid division by 0 pressure
143	ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
144
145	IF(DELTAT.LT.DTIME)THEN
146	PRINT *,'GWD RANDO: DELTAT LT DTIME!'
147	STOP
148	ENDIF
149
150
151	IF (NLEV < NW) THEN
152	PRINT *, 'YOU WILL HAVE PROBLEM WITH RANDOM NUMBERS'
153	PRINT *, 'FLOTT GWD STOP'
154	STOP 1
155	ENDIF
156
157	! 1.2 WAVES CHARACTERISTICS CHOSEN RANDOMLY
158	!-------------------------------------------
159
160	! The mod function of here a weird arguments
161	! are used to produce the waves characteristics
162	! in a stochastic way
163
164	JW = 0
165	DO JP = 1, NP
166	DO JK = 1, NK
167	DO JO = 1, NO
168	JW = JW + 1
169	! Angle
170	ZP(JW) = 2. * RPI * REAL(JP - 1) / REAL(NP)
171	DO II = 1, KLON
172	! Horizontal wavenumber amplitude
173	ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
174	! Horizontal phase speed
175	CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
176	!Online output: linear
177	if (output) CPHA = CMIN + (CMAX - CMIN) * JO/NO
178	! Intrinsic frequency
179	ZO(JW, II) = CPHA * ZK(JW, II)
180	! Momentum flux at launch lev
181	! RUW0(JW, II) = RUWMAX / REAL(NW) &
182	RUW0(JW, II) = RUWMAX &
183	* MOD(100. * (UU(II, JW)2 + VV(II, JW)2), 1.)
184	ENDDO
185	end DO
186	end DO
187	end DO
188
189	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
190	!-------------------------------------------------------------
191
192	IEQ = KLON / 2
193	!Online output
194	if (output) OPEN(11,file="impact-gwno.dat")
195
196	! Pressure and Inv of pressure, Temperature / at 1/2 level
197	DO LL = 2, KLEV
198	PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
199	end DO
200
201	PH(:, KLEV + 1) = 0.
202	PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
203
204	! Launching altitude
205
206	DO LL = 1, NLEV
207	IF (PH(IEQ, LL) / PH(IEQ, 1) > XLAUNCH) LAUNCH = LL
208	ENDDO
209
210	! Log pressure vert. coordinate (altitude above surface)
211	ZH(:,1) = 0.
212	DO LL = 2, KLEV + 1
213	H0(:, LL-1) = RD * TT(:, LL-1) / RG
214	ZH(:, LL) = ZH(:, LL-1) + H0(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
215	end DO
216
217	! Winds and BV frequency
218	DO LL = 2, KLEV
219	UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
220	VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
221	! BVSEC: BV Frequency
222	! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
223	BV(:, LL) = SQRT(MAX(BVSEC,pn2(:,LL)))
224	end DO
225	BV(:, 1) = BV(:, 2)
226	UH(:, 1) = 0.
227	VH(:, 1) = 0.
228	BV(:, KLEV + 1) = BV(:, KLEV)
229	UH(:, KLEV + 1) = UU(:, KLEV)
230	VH(:, KLEV + 1) = VV(:, KLEV)
231
232
233	! 3. COMPUTE THE FLUXES
234	!--------------------------
235
236	! 3.1 Vertical velocity at launching altitude to ensure
237	! the correct value to the imposed fluxes.
238	!
239	DO JW = 1, NW
240
241	! Evaluate intrinsic frequency at launching altitude:
242	ZOP(JW, :) = ZO(JW, :) &
243	- ZK(JW, :) * COS(ZP(JW)) * UH(:, LAUNCH) &
244	- ZK(JW, :) * SIN(ZP(JW)) * VH(:, LAUNCH)
245	! Vertical velocity at launch level, value to ensure the imposed
246	! mom flux:
247	WWP(JW, :) = SQRT(ABS(ZOP(JW, :)) / MAX(BV(:, LAUNCH),BVSEC) &
248	* RUW0(JW,:))
249	RUWP(JW, :) = COS(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
250	RVWP(JW, :) = SIN(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
251
252	end DO
253
254	! 3.2 Uniform values below the launching altitude
255
256	DO LL = 1, LAUNCH
257	RUW(:, LL) = 0
258	RVW(:, LL) = 0
259	DO JW = 1, NW
260	RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
261	RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
262	end DO
263	end DO
264
265	! 3.3 Loop over altitudes, with passage from one level to the
266	! next done by i) conserving the EP flux, ii) dissipating
267	! a little, iii) testing critical levels, and vi) testing
268	! the breaking.
269
270	!Online output
271	write(str2,'(i2)') NW+1
272	outform="("//str2//"(E12.4,1X))"
273	if (output) WRITE(11,outform) ZH(IEQ, 1) / 1000., (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW)), JW = 1, NW)
274
275	DO LL = LAUNCH, KLEV - 1
276
277
278	! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL
279	! TO THE NEXT)
280	DO JW = 1, NW
281	ZOM(JW, :) = ZOP(JW, :)
282	WWM(JW, :) = WWP(JW, :)
283	! Intrinsic Frequency
284	ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW)) * UH(:, LL + 1) &
285	- ZK(JW, :) * SIN(ZP(JW)) * VH(:, LL + 1)
286
287	WWP(JW, :) = MIN( &
288	! No breaking (Eq.6)
289	WWM(JW, :) &
290	* SQRT(BV(:, LL ) / BV(:, LL+1) &
291	* ABS(ZOP(JW, :)) / MAX(ABS(ZOM(JW, :)), ZOISEC)) &
292	! Dissipation (Eq. 8):
293	* EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
294	* ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
295	/ MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
296	* ZK(JW, :)*3 (ZH(:, LL + 1) - ZH(:, LL))), &
297	! Critical levels (forced to zero if intrinsic
298	! frequency changes sign)
299	MAX(0., SIGN(1., ZOP(JW, :) * ZOM(JW, :))) &
300	! Saturation (Eq. 12)
301	* ZOP(JW, :)**2 / ZK(JW, :)/BV(:, LL+1) &
302	* EXP(-ZH(:, LL + 1)/2./H0(:,LL)) * SAT)
303	end DO
304
305	! END OF W(KB)ARNING
306	! Evaluate EP-flux from Eq. 7 and
307	! Give the right orientation to the stress
308
309	DO JW = 1, NW
310	RUWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
311	*BV(:,LL+1)&
312	* COS(ZP(JW)) * MAX(WWP(JW, :),1e-30)**2
313	RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
314	*BV(:,LL+1)&
315	* SIN(ZP(JW)) * MAX(WWP(JW, :),1e-30)**2
316	end DO
317	!
318	RUW(:, LL + 1) = 0.
319	RVW(:, LL + 1) = 0.
320
321	DO JW = 1, NW
322	RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
323	RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
324	end DO
325	!Online output
326	if (output) WRITE(11,outform) ZH(IEQ, LL + 1) / 1000., (RUWP(JW, IEQ), JW = 1, NW)
327
328	end DO
329
330	!Online output
331	if (output) then
332	CLOSE(11)
333	stop
334	endif
335
336	! 4 CALCUL DES TENDANCES:
337	!------------------------
338
339	! 4.1 Rectification des flux au sommet et dans les basses couches:
340
341	RUW(:, KLEV + 1) = 0.
342	RVW(:, KLEV + 1) = 0.
343	RUW(:, 1) = RUW(:, LAUNCH)
344	RVW(:, 1) = RVW(:, LAUNCH)
345	DO LL = 2, LAUNCH
346	RUW(:, LL) = RUW(:, LL - 1) + (RUW(:, LAUNCH + 1) - RUW(:, 1)) * &
347	(PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
348	RVW(:, LL) = RVW(:, LL - 1) + (RVW(:, LAUNCH + 1) - RVW(:, 1)) * &
349	(PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
350	end DO
351
352	! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
353	DO LL = 1, KLEV
354	d_u(:, LL) = RG * (RUW(:, LL + 1) - RUW(:, LL)) &
355	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
356	d_v(:, LL) = RG * (RVW(:, LL + 1) - RVW(:, LL)) &
357	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
358	ENDDO
359	! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
360	d_u = DTIME/DELTAT/REAL(NW) * d_u + (1.-DTIME/DELTAT) * d_u_sav
361	d_v = DTIME/DELTAT/REAL(NW) * d_v + (1.-DTIME/DELTAT) * d_v_sav
362	d_u_sav = d_u
363	d_v_sav = d_v
364
365	! Cosmetic: evaluation of the cumulated stress
366
367	ZUSTR(:) = 0.
368	ZVSTR(:) = 0.
369	DO LL = 1, KLEV
370	ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
371	ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
372	ENDDO
373
374	END SUBROUTINE FLOTT_GWD_RAN

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

Download in other formats:

Original Format