Context Navigation

nonoro_gwd_ran_mod.F90 @ 3523

Last change on this file since 3523 was 2595, checked in by emillour, 3 years ago

Generic GCM:
Fixes and improvements in the Non-orographic GW scheme, namely:

increments are not tendencies
missing rho factor in EP flux computation
missing rho at launch altitude
changed inputs, because R and Cp are needed to compute rho and BV

File size: 26.2 KB

Line
1	MODULE nonoro_gwd_ran_mod
2
3	IMPLICIT NONE
4
5	REAL, allocatable, save :: du_nonoro_gwd(:, :) ! Zonal wind tendency due to GWD
6	REAL, allocatable, save :: dv_nonoro_gwd(:, :) ! Meridional wind tendency due to GWD
7	REAL, allocatable, save :: east_gwstress(:, :) ! Eastward stress profile
8	REAL, allocatable, save :: west_gwstress(:, :) ! Westward stress profile
9	!$OMP THREADPRIVATE(du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
10	CONTAINS
11
12	SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, cpnew, rnew, pp, &
13	zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
14	zustr, zvstr, d_t, d_u, d_v)
15
16	!--------------------------------------------------------------------------------
17	! Parametrization of the momentum flux deposition due to a discrete
18	! number of gravity waves.
19	! F. Lott
20	! Version 14, Gaussian distribution of the source
21	! LMDz model online version
22	! ADAPTED FOR VENUS / F. LOTT + S. LEBONNOIS
23	! Version adapted on 03/04/2013:
24	! - input flux compensated in the deepest layers
25	!
26	! ADAPTED FOR MARS G.GILLI 02/2016
27	! Revision with F.Forget 06/2016 Variable EP-flux
28	! according to
29	! PBL variation (max
30	! velocity thermals)
31	! D.BARDET 01/2020 - reproductibility of
32	! the launching altitude
33	! calculation
34	! - wave characteristic
35	! calculation using MOD
36	! - adding east_gwstress
37	! and west_gwstress variables
38	!
39	! ADAPTED FOR GENERIC D.BARDET 01/2020
40	! UPDATED J.LIU 12/2021 The rho (density) at the specific
41	! locations is introduced. The equation
42	! of EP-flux is corrected by adding the
43	! term of density at launch (source)
44	! altitude(level). Bugs in BV,d_u,d_v fixed.
45	!---------------------------------------------------------------------------------
46
47
48
49	use comcstfi_mod, only: g, pi, r
50	USE ioipsl_getin_p_mod, ONLY : getin_p
51	use vertical_layers_mod, only : presnivs
52	use geometry_mod, only: cell_area
53	#ifdef CPP_XIOS
54	use xios_output_mod, only: send_xios_field
55	#endif
56	implicit none
57
58	! 0. DECLARATIONS:
59
60	! 0.1 INPUTS
61
62	INTEGER, intent(in):: ngrid ! number of atmospheric columns
63	INTEGER, intent(in):: nlayer ! number of atmospheric columns
64	REAL, intent(in):: dtime ! Time step of the Physics
65	REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m)
66	REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere
67	REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere
68	REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels(Pa)
69	REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels(K)
70	REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels(m/s)
71	REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels(m/s)
72	REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature(K/s)
73	REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind(m/s/s)
74	REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind(m/s/s)
75
76	! 0.2 OUTPUTS
77
78	REAL, intent(out):: zustr(ngrid) ! Zonal surface stress
79	REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress
80	REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature (K/s) due to gravity waves (not used set to zero)
81	REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind (m/s/s) due to gravity waves
82	REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves
83
84	! 0.3 INTERNAL ARRAYS
85
86	REAL :: tt(ngrid, nlayer) ! Temperature at full levels
87	REAL :: rho(ngrid, nlayer) ! Mass density at full levels
88	REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels
89	REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels
90	REAL :: bvlow(ngrid) ! initial N at each grid (not used)
91	REAL :: dz ! depth of the GW source if gwd_convective_source=.true.
92
93	INTEGER ii, jj, ll
94
95	! 0.3.0 Time scale of the like cycle of the waves parametrized
96	REAL, parameter:: deltat = 24. * 3600.
97
98	! 0.3.1 Gravity waves specifications
99	INTEGER, parameter:: nk = 2 ! number of horizontal wavenumbers
100	INTEGER, parameter:: np = 2 ! directions (eastward and westward) phase speed
101	INTEGER, parameter:: no = 2 ! absolute values of phase speed
102	INTEGER, parameter:: na = 5 ! Number of realizations to get the phase speed
103	INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves
104	INTEGER jk, jp, jo, jw
105
106	REAL kstar ! Control value to constrain the min horizontal wavenumber by the horizontal resolution
107	REAL, save :: kmax ! Max horizontal wavenumber=N/u, u(=30) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax=62.8 km
108	!$OMP THREADPRIVATE(kmax)
109	!REAL, parameter:: kmax = 4.e-5 ! Max horizontal wavenumber=N/u, u(=70) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax
110	REAL, parameter:: kmin = 2.e-6 ! Min horizontal wavenumber=1/sqrt(DxDy) Dx and Dy horizontal grid,lambda max,lambda=2*pi/kmax=314 km
111	REAL, save :: cmax ! Max horizontal absolute phase velocity=zonal wind at launch
112	!$OMP THREADPRIVATE(cmax)
113	!REAL, parameter:: cmax = 100. ! Test for Saturn: Max horizontal absolute phase velocity
114	!REAL, parameter:: cmax = 70. ! Test for Saturn: Max horizontal absolute phase velocity
115	REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity(not used)
116	REAL cpha ! absolute phase velocity frequency
117	REAL zk(nw, ngrid) ! horizontal wavenumber amplitude
118	REAL zp(nw, ngrid) ! horizontal wavenumber angle
119	REAL zo(nw, ngrid) ! absolute frequency
120
121	REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
122	REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
123	REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level
124	REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level
125	REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
126	REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
127	REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW)
128	REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW)
129	REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0)
130	REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable)
131	!$OMP THREADPRIVATE(epflux_max)
132	INTEGER launch ! Launching altitude
133	REAL, parameter:: xlaunch = 0.2 ! Control the launching altitude by pressure
134	REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer in m (approx.)
135
136	REAL prec(ngrid) ! precipitations
137	REAL prmax ! Max value of precipitation
138
139	! 0.3.2 Parameters of waves dissipations
140	REAL, save :: sat ! saturation parameter(tunable)
141	!$OMP THREADPRIVATE(sat)
142	REAL, save :: rdiss ! coefficient of dissipation(tunable)
143	!$OMP THREADPRIVATE(rdiss)
144	REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic freguency
145
146	! 0.3.3 Background flow at 1/2 levels and vertical coordinate
147	REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere
148	REAL, save:: H0 ! characteristic height of the atmosphere
149	!$OMP THREADPRIVATE(H0)
150	!REAL, parameter:: pr = 250 ! Reference pressure [Pa]
151	!REAL, parameter:: tr = 220. ! Reference temperature [K]
152	! TEST FOR SATURN STRATOSPHERE
153	REAL, parameter:: pr = 110 ! Reference pressure [Pa]
154	REAL, parameter:: tr = 130. ! Reference temperature [K]
155	REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0)
156	REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H)
157	REAL uh(ngrid, nlayer + 1) ! Zonal wind at 1/2 levels
158	REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels
159	REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels
160	REAL, parameter:: psec = 1.e-12 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure)
161	REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels
162	REAL, parameter:: bvsec = 1.e-6 ! Security to avoid negative BV (!!IMPORTANT: tunable, depends on which Planet)
163	REAL href(nlayer + 1) ! Reference altitude for launching altitude
164
165
166
167	REAL ran_num_1, ran_num_2, ran_num_3
168
169
170	LOGICAL, save :: firstcall = .true.
171	!$OMP THREADPRIVATE(firstcall)
172	LOGICAL, save :: gwd_convective_source = .false.
173	!$OMP THREADPRIVATE(gwd_convective_source)
174
175	!--------------------------------------------------------------------------------------------------
176	! 1. INITIALISATION
177	!--------------------------------------------------------------------------------------------------
178	IF (firstcall) THEN
179	write(,) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
180	epflux_max = 0.0 ! Default value !!
181	call getin_p("nonoro_gwd_epflux_max", epflux_max)
182	write(,) "nonoro_gwd_ran: epflux_max=", epflux_max
183	sat = 1. ! default gravity waves saturation value !!
184	call getin_p("nonoro_gwd_sat", sat)
185	write(,) "nonoro_gwd_ran: sat=", sat
186	cmax = 30. ! default gravity waves phase velocity value !!
187	call getin_p("nonoro_gwd_cmax", cmax)
188	write(,) "nonoro_gwd_ran: cmax=", cmax
189	rdiss=1 ! default coefficient of dissipation !!
190	call getin_p("nonoro_gwd_rdiss", rdiss)
191	write(,) "nonoro_gwd_ran: rdiss=", rdiss
192	kmax = 1.e-4 ! default Max horizontal wavenumber !!
193	call getin_p("nonoro_gwd_kmax", kmax)
194	write(,) "nonoro_gwd_ran: kmax=", kmax
195
196	! Characteristic vertical scale height
197	H0 = r * tr / g
198	! Control
199	if (deltat .LT. dtime) THEN
200	call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
201	endif
202	if (nlayer .LT. nw) THEN
203	call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
204	endif
205	firstcall = .false.
206	ENDIF
207	gwd_convective_source = .false.
208
209	! Compute subroutine's current values of temperature and winds
210	tt(:,:) = pt(:,:) + dtime * pdt(:,:)
211	uu(:,:) = pu(:,:) + dtime * pdu(:,:)
212	vv(:,:) = pv(:,:) + dtime * pdv(:,:)
213	! Compute the real mass density by rho=p/R(T)T
214	DO ll=1,nlayer
215	DO ii=1,ngrid
216	rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll))
217	ENDDO
218	ENDDO
219	! print*,'epflux_max just after firstcall:', epflux_max
220
221	!--------------------------------------------------------------------------------------------------
222	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
223	!--------------------------------------------------------------------------------------------------
224	! Pressure and Inv of pressure, temperature at 1/2 level
225	DO ll = 2, nlayer
226	ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.)
227	end DO
228
229	ph(:, nlayer + 1) = 0.
230	ph(:, 1) = 2. * pp(:, 1) - ph(:, 2)
231
232	! Launching altitude for reproductible case
233	DO ll = 2, nlayer
234	href(ll) = exp((log(presnivs(ll)) + log(presnivs(ll - 1))) / 2.)
235	end DO
236	href(nlayer + 1) = 0.
237	href(1) = 2. * presnivs(1) - href(2)
238
239	launch = 0.
240	DO ll =1, nlayer
241	IF (href (ll) / href(1) > xlaunch) launch = ll
242	end DO
243
244	! Log-pressure vertical coordinate
245	!DO ll = 1, nlayer + 1
246	! zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec))
247	!end DO
248	ZH(:,1) = 0.
249	DO LL = 2, nlayer + 1
250	!ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
251	H0bis(:, LL-1) = r * TT(:, LL-1) / g
252	ZH(:, LL) = ZH(:, LL-1) &
253	+ H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
254	end DO
255	ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC))
256
257	! Winds and Brunt Vaisala frequency
258	DO ll = 2, nlayer
259	uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind
260	vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind
261
262	BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1))) &
263	*((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+ &
264	G / cpnew(:,LL))
265
266	BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive
267	BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small
268	end DO
269
270	bv(:, 1) = bv(:, 2)
271	uh(:, 1) = 0.
272	vh(:, 1) = 0.
273	bv(:, nlayer + 1) = bv(:, nlayer)
274	uh(:, nlayer + 1) = uu(:, nlayer)
275	vh(:, nlayer + 1) = vv(:, nlayer)
276
277	!-----------------------------------------------------------------------------------------------------------------------
278	! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
279	!-----------------------------------------------------------------------------------------------------------------------
280	! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way
281	DO jw = 1, nw
282	!Angle
283	DO ii = 1, ngrid
284
285	ran_num_1 = mod(tt(ii, jw) * 10., 1.)
286	ran_num_2 = mod(tt(ii, jw) * 100., 1.)
287
288	! Angle (0 or pi so far)
289	zp(jw, ii) = (sign(1., 0.5 - ran_num_1) + 1.) * pi / 2.
290	! Horizontal wavenumber amplitude
291	! TN+GG April/June 2020 (rev2381)
292	! "Individual waves are not supposed to occupy the entire domain, but only a faction of it" Lott+2012
293	!zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
294	kstar = pi / sqrt(cell_area(ii))
295	zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2
296
297	! Horizontal phase speed
298	! this computation allows to get a gaussian distribution centered on 0 (right ?)
299	! then cmin is not useful, and it favors small values.
300	cpha = 0.
301	DO jj = 1, na
302	ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.)
303	cpha = cpha + 2. * cmax * &
304	(ran_num_3 - 0.5) * &
305	sqrt(3.) / sqrt(na * 1.)
306	end DO
307	IF (cpha < 0.) THEN
308	cpha = - 1. * cpha
309	zp (jw, ii) = zp(jw, ii) + pi
310	ENDIF
311	! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
312	! ran_num_3 = mod(tt(ii, jw)**2, 1.)
313	! cpha = cmin + (cmax - cmin) * ran_num_3
314
315	! Intrinsic frequency
316	zo(jw, ii) = cpha * zk(jw, ii)
317	! Intrinsic frequency is imposed
318	zo(jw, ii) = zo(jw, ii) &
319	+ zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) &
320	+ zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch)
321	! Momentum flux at launch level
322	epflux_0(jw, ii) = epflux_max &
323	* mod(100. * (uu(ii, jw)2 + vv(ii, jw)2), 1.)
324	ENDDO
325	end DO !DO jw = 1, nw
326
327	! print*,'epflux_max just after waves charac. ramdon:', epflux_max
328	! print*,'epflux_0 just after waves charac. ramdon:', epflux_0
329
330	!-----------------------------------------------------------------------------------------------------------------------
331	! 4. COMPUTE THE FLUXES
332	!-----------------------------------------------------------------------------------------------------------------------
333	! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes.
334	!------------------------------------------------------
335	DO jw = 1, nw
336	! Evaluate intrinsic frequency at launching altutide:
337	intr_freq_p(jw, :) = zo(jw, :) &
338	- zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) &
339	- zk(jw, :) * sin(zp(jw, :)) * vh(:, launch)
340	end DO
341
342	IF (gwd_convective_source) THEN
343	DO jw = 1, nw
344	! VERSION WITH CONVECTIVE SOURCE ! designed for Earth
345
346	! Vertical velocity at launch level, value to ensure the
347	! imposed mmt flux factor related to the convective forcing:
348	! precipitations.
349
350	! tanh limitation to values above prmax:
351	! WWP(JW, :) = epflux_0(JW, :) &
352	! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
353	! Here, we neglected the kinetic energy providing of the thermodynamic
354	! phase change
355
356	!
357
358	! Factor related to the characteristics of the waves:
359	wwp(jw, :) = wwp(jw, :) * zk(jw, :)**3 / kmin / bvlow(:) &
360	/ MAX(ABS(intr_freq_p(jw, :)), zoisec)**3
361
362	! Moderation by the depth of the source (dz here):
363	wwp(jw, :) = wwp(jw, :) &
364	* exp(- bvlow(:)2 / max(abs(intr_freq_p(jw, :)), zoisec)2 &
365	* zk(jw, :)*2 dz**2)
366
367	! Put the stress in the right direction:
368	u_epflux_p(jw, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 &
369	* bv(:, launch) * cos(zp(jw, :)) * wwp(jw, :)**2
370	v_epflux_p(JW, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 &
371	* bv(:, launch) * sin(zp(jw, :)) * wwp(jw, :)**2
372	end DO
373	ELSE ! VERSION WITHOUT CONVECTIVE SOURCE
374	! Vertical velocity at launch level, value to ensure the imposed
375	! mom flux:
376	DO jw = 1, nw
377	! WW is directly a flux, here, not vertical velocity anymore
378	wwp(jw, :) = epflux_0(JW,:)
379	u_epflux_p(jw, :) = cos(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :)
380	v_epflux_p(jw, :) = sin(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :)
381
382	end DO
383	ENDIF
384
385	! 4.2 Initial flux at launching altitude
386	!---------------------------------------
387	u_epflux_tot(:, launch) = 0.
388	v_epflux_tot(:, launch) = 0.
389	DO jw = 1, nw
390	u_epflux_tot(:, launch) = u_epflux_tot(:, launch) + u_epflux_p(jw, :)
391	v_epflux_tot(:, launch) = v_epflux_tot(:, launch) + v_epflux_p(jw, :)
392	end DO
393
394	! 4.3 Loop over altitudes, with passage from one level to the next done by:
395	!--------------------------------------------------------------------------
396	! i) conserving the EP flux,
397	! ii) dissipating a little,
398	! iii) testing critical levels,
399	! iv) testing the breaking.
400	!--------------------------------------------------------------------------
401	DO ll = launch, nlayer - 1
402	! Warning! all the physics is here (passage from one level to the next
403	DO jw = 1, nw
404	intr_freq_m(jw, :) = intr_freq_p(jw, :)
405	wwm(jw, :) = wwp(jw, :)
406	! Intrinsic frequency
407	intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) &
408	- zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
409	! Vertical velocity in flux formulation
410	wwp(jw, :) = min( &
411	! No breaking (Eq. 6):
412	wwm(jw, :) &
413	! Dissipation (Eq. 8)! (real rho used here rather than pressure
414	! parametration because the original code has a bug if the density of
415	! the planet at the launch altitude not approximate 1):
416	* exp(-rdiss * 2./rho(:, LL + 1) &
417	* ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 &
418	/ max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 &
419	* zk(jw, :)*3 (zh(:, ll + 1) - zh(:, ll))) , &
420	! Critical levels (forced to zero if intrinsic frequency
421	! changes sign)
422	max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) &
423	! Saturation (Eq. 12)(rho at launch altitude is imposed by J.Liu.
424	! Same reason with Eq. 8)
425	* abs(intr_freq_p(jw, :))*3 /(2.0 bv(:, ll + 1)) &
426	* rho(:,launch)*exp(-zh(:, ll + 1) / H0) &
427	* sat*2 kmin2 / zk(jw, :)4)
428	end DO
429	! END OF W(KB)ARNING
430
431	! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress
432	DO jw = 1, nw
433	u_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * cos(zp(jw, :)) * wwp(jw, :)
434	v_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * sin(zp(jw, :)) * wwp(jw, :)
435	end DO
436
437	u_epflux_tot(:, ll + 1) = 0.
438	v_epflux_tot(:, ll + 1) = 0.
439	DO jw = 1, nw
440	u_epflux_tot(:, ll + 1) = u_epflux_tot(:, ll + 1) + u_epflux_p(jw, :)
441	v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :)
442	east_gwstress(:, ll) = east_gwstress(:, ll) &
443	+ max(0., u_epflux_p(jw, ll)) / float(nw)
444	west_gwstress(:, ll) = west_gwstress(:, ll) &
445	+ min(0., u_epflux_p(jw, ll)) / float(nw)
446	end DO
447	end DO ! DO LL = LAUNCH, nlayer - 1
448
449	! 5. TENDENCY CALCULATIONS
450	!-------------------------
451	! 5.1 Flow rectification at the top and in the low layers
452	! --------------------------------------------------------
453	! Warning, this is the total on all GW...
454
455	u_epflux_tot(:, nlayer + 1) = 0.
456	v_epflux_tot(:, nlayer + 1) = 0.
457
458	! Here, big change compared to FLott version:
459	! We compensate (u_epflux_(:, LAUNCH), ie total emitted upward flux
460	! over the layers max(1,LAUNCH-3) to LAUNCH-1
461	DO LL = 1, max(1,LAUNCH-3)
462	u_epflux_tot(:, LL) = 0.
463	v_epflux_tot(:, LL) = 0.
464	end DO
465	DO LL = max(2,LAUNCH-2), LAUNCH-1
466	u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) &
467	+ u_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
468	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
469	v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) &
470	+ v_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
471	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
472	east_gwstress(:,LL) = east_gwstress(:, LL - 1) &
473	+ east_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
474	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
475	west_gwstress(:,LL) = west_gwstress(:, LL - 1) &
476	+ west_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
477	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
478	end DO
479	! This way, the total flux from GW is zero, but there is a net transport
480	! (upward) that should be compensated by circulation
481	! and induce additional friction at the surface
482
483	! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
484	!---------------------------------------------
485	DO LL = 1, nlayer
486	!Notice here the d_u and d_v are tendency (i.e. -1/rho*dEP/dz For Mars) but not increment(Venus).
487	d_u(:, LL) = - (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
488	/ (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL)))
489	d_v(:, LL) = - (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
490	/ (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL)))
491	ENDDO
492	d_t(:,:) = 0.
493
494	! 5.3 Update tendency of wind with the previous (and saved) values
495	!-----------------------------------------------------------------
496	d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) &
497	+ (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
498	d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) &
499	+ (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
500	du_nonoro_gwd(:,:) = d_u(:,:)
501	dv_nonoro_gwd(:,:) = d_v(:,:)
502
503	! Cosmetic: evaluation of the cumulated stress
504	ZUSTR(:) = 0.
505	ZVSTR(:) = 0.
506	DO LL = 1, nlayer
507	ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
508	ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
509	ENDDO
510
511	#ifdef CPP_XIOS
512	call send_xios_field("du_nonoro", d_u)
513	call send_xios_field("dv_nonoro", d_v)
514	call send_xios_field("east_gwstress", east_gwstress)
515	call send_xios_field("west_gwstress", west_gwstress)
516	#endif
517
518	END SUBROUTINE nonoro_gwd_ran
519
520	! ===================================================================
521	! Subroutines used to write variables of memory in start files
522	! ===================================================================
523
524	SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
525
526	IMPLICIT NONE
527
528	INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns
529	INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
530
531	allocate(du_nonoro_gwd(ngrid, nlayer))
532	allocate(dv_nonoro_gwd(ngrid, nlayer))
533	allocate(east_gwstress(ngrid, nlayer))
534	allocate(west_gwstress(ngrid, nlayer))
535
536	END SUBROUTINE ini_nonoro_gwd_ran
537	! ----------------------------------
538	SUBROUTINE end_nonoro_gwd_ran
539
540	IMPLICIT NONE
541
542	if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
543	if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
544	if (allocated(east_gwstress)) deallocate(east_gwstress)
545	if (allocated(west_gwstress)) deallocate(west_gwstress)
546
547	END SUBROUTINE end_nonoro_gwd_ran
548
549	END MODULE nonoro_gwd_ran_mod
550

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

Download in other formats:

Original Format