Context Navigation

nonoro_gwd_ran_mod.F90 @ 3556

Last change on this file since 3556 was 3035, checked in by amartinez, 16 months ago

======= VENUS PCM ===============

COMMIT BY ANTOINE MARTINEZ

Implementation of vertical ambipolar diffusion in physics

=================================

NEW KEYWORD OF PHYSIQ.DEF

=================================

NEW VERSION OF "physiq.def"

deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def

ok_iondiff: keyword supposed to activate ion ambipolar diffusion

nbapp_chem: replaces nbapp_chim, in order to characterize the Number of calls to the chemistry routines (per Venusian day)

================

phyvenus

================

Iondiff_red.F

Calculation of the Ion Ambipolar Diffusion for 13 ions on 14:

O2+, O+, C+, N+, CO2+,
NO+, CO+, H+, H2O+, H3O+,

HCO+, N2+, OH+

The ion temperature is fixed as the half of the electron temperature
calculated in ion_h.F for the stability of the program and because the
ion temperature is lower than electron temperature.

plasmaphys_venus_mod.F90

parameters of the ambipolar diffusion scheme:

parameter (Pdiffion = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
parameter (naltvert = 168) ! number of level on the altimetric subgrid
parameter (nsubvert = 24000) ! ptimephysiq/nsubvert - minimum sub-timestep allowed
parameter ( nsubmin = 40) ! ptimephysiq/nsubmin - maximum sub-timestep allowed

physic_mod.F

nbapp_chem is not fixed anymore here but deftank/in physiq.def
Ambipolar diffusion activated if (ok_iondiff) is true

conf_phys.f90

add ok_iondiff as parameters to read from physiq.def file set to .false. by default
add nbapp_chem as parameters to characterize the Number of calls to the chemistry routines (per Venusian day) instead of be fixed in physic_mod.F

to read from physiq.def file set to 24000 by default

cleshphys.h

added ok_iondiff & nbapp_chem in COMMON/clesphys_l/
removed nbapp_chim

phytrac_chemistry.F

Added security in the calculation of the sza_local in order to avoid the rare case where the |range| is above 1

File size: 23.6 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, pp, &
13	pn2, presnivs, pt, pu, pv, &
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	! READAPTED FOR VENUS S.LEBONNOIS 08/2021
41	!---------------------------------------------------------------------------------
42
43	USE ioipsl_getin_p_mod, ONLY : getin_p
44	USE geometry_mod, ONLY: cell_area, latitude_deg
45	!#ifdef CPP_XIOS
46	! use xios_output_mod, only: send_xios_field
47	!#endif
48	implicit none
49	#include "YOMCST.h"
50
51	! 0. DECLARATIONS:
52
53	! 0.1 INPUTS
54
55	INTEGER, intent(in):: ngrid, nlayer
56	REAL, intent(in):: dtime ! Time step of the Physics
57	REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels
58	REAL, intent(in):: pn2(ngrid, nlayer) ! N2 (BV^2) at 1/2 levels
59	REAL, intent(in):: presnivs(nlayer) ! Approximate pressure for reproductible launching altitude
60	REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels
61	REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels
62	REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels
63
64	! 0.2 OUTPUTS
65
66	REAL, intent(out):: zustr(ngrid) ! Zonal surface stress
67	REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress
68	REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature
69	REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind
70	REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind
71
72	! 0.3 INTERNAL ARRAYS
73
74	REAL :: tt(ngrid, nlayer) ! Temperature at full levels
75	REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels
76	REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels
77
78	INTEGER ii, jj, ll
79
80	! 0.3.0 Time scale of the like cycle of the waves parametrized
81	REAL, parameter:: deltat = 24. * 3600.
82
83	! 0.3.1 Gravity waves specifications
84	INTEGER, parameter:: nk = 2 ! number of horizontal wavenumbers
85	INTEGER, parameter:: np = 2 ! directions (eastward and westward) phase speed
86	INTEGER, parameter:: no = 2 ! absolute values of phase speed
87	INTEGER, parameter:: na = 5 ! Number of realizations to get the phase speed
88	INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves
89	INTEGER jk, jp, jo, jw
90	REAL kstar ! Control value to constrain the min horizontal wavenumber by the horizontal resolution
91
92	!!! TEST Kobs and BestFit Diogo
93	REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber (lambda min 50 km)
94	! generic : kmax=N/u, u(=30) zonal wind velocity at launch
95	REAL, parameter:: kmin = 1.e-5 ! Min horizontal wavenumber (lambda max 628 km)
96	! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid
97
98	!----------------------------------------
99	! VCD 1.1 tuning
100	! REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity
101	!----------------------------------------
102	! VCD 2.0 tuning
103	REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity
104	!----------------------------------------
105
106	! generic : cmax=zonal wind at launch
107	REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity
108
109	!----------------------------------------
110	! Nominal as Gilli2017
111	! REAL, parameter:: epflux_max = 0.005 ! Max EP flux value at launching altitude (previous name: RUWMAX)
112	!----------------------------------------
113	! VCD 2.1 tuning
114	REAL, parameter:: epflux_max = 0.001 ! Max EP flux value at launching altitude (previous name: RUWMAX)
115	!----------------------------------------
116
117	REAL cpha ! absolute phase velocity frequency
118	REAL zk(nw, ngrid) ! horizontal wavenumber amplitude
119	REAL zp(nw, ngrid) ! horizontal wavenumber angle
120	REAL zo(nw, ngrid) ! absolute frequency
121
122	REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
123	REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
124	REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level
125	REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level
126	REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
127	REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
128	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)
129	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)
130	REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0)
131	INTEGER launch ! Launching altitude
132	REAL, parameter:: xlaunch = 5e-3 ! Value for top of cloud convective region
133
134	! 0.3.2 Parameters of waves dissipations
135	!----------------------------------------
136	! VCD 1.1 tuning
137	! REAL, parameter:: sat = 0.85 ! saturation parameter
138	! REAL, parameter:: rdiss = 0.1 ! coefficient of dissipation
139	!----------------------------------------
140	! VCD 2.0 tuning
141	REAL, parameter:: sat = 0.6 ! saturation parameter
142	REAL, parameter:: rdiss = 8.e-4 ! coefficient of dissipation
143	!----------------------------------------
144	REAL, parameter:: zoisec = 1.e-8 ! 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 min_k(ngrid) ! min(kstar,kmin)
149	REAL, save:: H0 ! characteristic height of the atmosphere
150	REAL, save:: MDR ! characteristic mass density
151	!----------------------------------------
152	! VCD 1.1 tuning
153	! REAL, parameter:: pr = 5e5 ! Reference pressure [Pa] ! VENUS: cloud layer
154	!----------------------------------------
155	! VCD 2.0 tuning
156	REAL, parameter:: pr = 5e4 ! Reference pressure [Pa] ! VENUS: cloud layer
157	!----------------------------------------
158	REAL, parameter:: tr = 300. ! Reference temperature [K] ! VENUS: cloud layer
159	REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0)
160	REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H)
161	REAL uh(ngrid, nlayer + 1) ! Zonal wind at 1/2 levels
162	REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels
163	REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels
164	REAL, parameter:: psec = 1.e-9 ! Security to avoid division by 0 pressure
165	REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels
166	REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV
167	REAL href(nlayer + 1) ! Reference altitude for launching altitude
168
169	REAL ran_num_1, ran_num_2, ran_num_3
170
171	LOGICAL, save :: firstcall = .true.
172
173	!--------------------------------------------------------------------------------------------------
174	! 1. INITIALISATION
175	!------------------
176	IF (firstcall) THEN
177	write(,) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
178	! Characteristic vertical scale height
179	H0 = RD * tr / RG
180	! Characteristic mass density at launch altitude
181	MDR = pr / ( RD * tr )
182	! Control
183	if (deltat .LT. dtime) THEN
184	call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
185	endif
186	if (nlayer .LT. nw) THEN
187	call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
188	endif
189	firstcall = .false.
190	ENDIF
191
192	! For VENUS, we use *_seri variables, that already integrate the different previous tendencies
193	tt(:,:) = pt(:,:)
194	uu(:,:) = pu(:,:)
195	vv(:,:) = pv(:,:)
196
197	! print*,'epflux_max just after firstcall:', epflux_max
198
199	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
200	!----------------------------------------------------
201	! Pressure and Inv of pressure, temperature at 1/2 level
202	DO ll = 2, nlayer
203	ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.)
204	end DO
205
206	ph(:, nlayer + 1) = 0.
207	ph(:, 1) = 2. * pp(:, 1) - ph(:, 2)
208
209	! Launching altitude for reproductible case
210	DO ll = 2, nlayer
211	href(ll) = exp((log(presnivs(ll)) + log(presnivs(ll - 1))) / 2.)
212	end DO
213	href(nlayer + 1) = 0.
214	href(1) = 2. * presnivs(1) - href(2)
215
216	launch = 0.
217	DO ll =1, nlayer
218	IF (href (ll) / href(1) > xlaunch) launch = ll
219	end DO
220
221	! Log-pressure vertical coordinate
222	DO ll = 1, nlayer + 1
223	zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec))
224	end DO
225
226	! Winds and Brunt Vaisala frequency
227	DO ll = 2, nlayer
228	uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind
229	vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind
230	! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
231	bv(:, ll) = max(bvsec,sqrt(pn2(:,ll)))
232	end DO
233
234	bv(:, 1) = bv(:, 2)
235	uh(:, 1) = 0.
236	vh(:, 1) = 0.
237	bv(:, nlayer + 1) = bv(:, nlayer)
238	uh(:, nlayer + 1) = uu(:, nlayer)
239	vh(:, nlayer + 1) = vv(:, nlayer)
240
241	! TN+GG April/June 2020
242	! "Individual waves are not supposed
243	! to occupy the entire domain, but
244	! only a faction of it" Lott+2012
245	! minimum value of k between kmin and kstar
246	DO ii = 1, ngrid
247	kstar = RPI / sqrt(cell_area(ii))
248	min_k(ii) = max(kmin,kstar)
249	end DO
250
251	! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
252	!-----------------------------------------
253	! The mod function of here a weird arguments
254	! are used to produce the waves characteristics
255	! in a stochastic way
256	DO jw = 1, nw
257	DO ii = 1, ngrid
258
259	ran_num_1 = mod(tt(ii, jw) * 10., 1.)
260	ran_num_2 = mod(tt(ii, jw) * 100., 1.)
261
262	!! OPTIONS GENERIC DIFF VENUS !!
263	! angle (random) - reference
264	zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
265	+ 1.) * RPI / 2.
266	! zp(jw, ii) = atan2(4.*vh(ii, launch),uh(ii, launch)) &
267	! + (sign(1., 0.5 - ran_num_1) + 1.) * RPI / 2.
268	! --------- TRY 00 -----------------------
269	!IF((abs(latitude_deg(ii)) .le. 55.)) THEN
270	! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. &
271	! * (10.) ! +/- 10 deg par rapport equateur
272	!ELSE IF ((abs(latitude_deg(ii)) .le. 75.)) THEN
273	! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. &
274	! * (10. +(90.-10.) &
275	! * (latitude_deg(ii)-55.)/20.)
276	!ELSE
277	! zp(jw, ii) = ran_num_1 * RPI
278	!ENDIF
279	!zp(jw, ii) = mod(zp(jw, ii),2.*RPI)
280	! ------ TRY 01-------------------------------
281	!IF((abs(latitude_deg(ii)) .le. 55)) THEN
282	! zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
283	! + 1.) * RPI / 2.
284	!ELSE
285	! zp(jw, ii) = ran_num_1 * RPI
286	!ENDIF
287	! ---- angle (0 or RPI) -----
288	!!zp(jw, ii) = RPI*mod(jw,2)
289
290	! horizontal wavenumber amplitude
291	! zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
292	zk(jw, ii) = min_k(ii) + (kmax - min_k(ii)) * ran_num_2
293
294	! horizontal phase speed
295	! this computation allows to get a gaussian distribution centered on 0 (right ?)
296	! then cmin is not useful, and it favors small values.
297	! VCD 2.0 tuning
298	cpha = 0.
299	DO jj = 1, na
300	ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.)
301	cpha = cpha + 2. * cmax * &
302	(ran_num_3 - 0.5) * &
303	sqrt(3.) / sqrt(na * 1.)
304	end DO
305	!cpha = cpha - uh(ii, launch)
306	IF (cpha < 0.) THEN
307	cpha = - 1. * cpha
308	zp (jw, ii) = zp(jw, ii) + RPI
309	ENDIF
310	!IF (cpha < 1.) THEN
311	! cpha = 1.
312	!ENDIF
313	! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
314	! ran_num_3 = mod(tt(ii, jw)**2, 1.)
315	! cpha = cmin + (cmax - cmin) * ran_num_3
316
317	! Intrinsic frequency
318	zo(jw, ii) = cpha * zk(jw, ii)
319	! Intrinsic frequency is imposed
320	zo(jw, ii) = zo(jw, ii) &
321	+ zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) &
322	+ zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch)
323
324	! Momentum flux at launch level
325	epflux_0(jw, ii) = epflux_max &
326	* mod(100. * (uu(ii, jw)2 + vv(ii, jw)2), 1.)
327
328	end DO
329	end DO
330	! print*,'epflux_0 just after waves charac. ramdon:', epflux_0
331
332	! 4. COMPUTE THE FLUXES
333	!----------------------
334	! 4.1 Vertical velocity at launching altitude to ensure
335	! the correct value to the imposed fluxes.
336	!------------------------------------------------------
337	DO jw = 1, nw
338	! Evaluate intrinsic frequency at launching altutide:
339	intr_freq_p(jw, :) = zo(jw, :) &
340	- zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) &
341	- zk(jw, :) * sin(zp(jw, :)) * vh(:, launch)
342	end DO
343
344	! Vertical velocity at launch level, value to ensure the imposed
345	! mom flux:
346	DO jw = 1, nw
347	! WW is directly a flux, here, not vertical velocity anymore
348	wwp(jw, :) = epflux_0(jw,:)
349	u_epflux_p(jw, :) = cos(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :)
350	v_epflux_p(jw, :) = sin(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :)
351	end DO
352
353	! 4.2 Initial flux at launching altitude
354	!---------------------------------------
355	u_epflux_tot(:, launch) = 0.
356	v_epflux_tot(:, launch) = 0.
357	DO jw = 1, nw
358	u_epflux_tot(:, launch) = u_epflux_tot(:, launch) + u_epflux_p(jw, :)
359	v_epflux_tot(:, launch) = v_epflux_tot(:, launch) + v_epflux_p(jw, :)
360	end DO
361
362	! 4.3 Loop over altitudes, with passage from one level to the next done by:
363	!--------------------------------------------------------------------------
364	! i) conserving the EP flux,
365	! ii) dissipating a little,
366	! iii) testing critical levels,
367	! iv) testing the breaking.
368	!--------------------------------------------------------------------------
369	DO ll = launch, nlayer - 1
370	! Warning! all the physics is here (passage from one level to the next
371	DO jw = 1, nw
372	intr_freq_m(jw, :) = intr_freq_p(jw, :)
373	wwm(jw, :) = wwp(jw, :)
374	! Intrinsic frequency
375	intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) &
376	- zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
377	! Vertical velocity in flux formulation
378	wwp(jw, :) = min( &
379	! No breaking (Eq. 6):
380	wwm(jw, :) &
381	! Dissipation (Eq. 8):
382	* exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) &
383	* ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 &
384	/ max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 &
385	* zk(jw, :)*3 (zh(:, ll + 1) - zh(:, ll))) , &
386	! Critical levels (forced to zero if intrinsic frequency
387	! changes sign)
388	max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) &
389	! Saturation (Eq. 12)
390	* abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) &
391	* exp(-zh(:, ll + 1) / H0) &
392	!! Correction for VCD 2.0
393	! * sat*2 kmin2 / zk(jw, :)4)
394	* sat*2 min_k(:)*2 MDR / zk(jw, :)**4)
395	end DO
396
397	! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress
398	DO jw = 1, nw
399	u_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * cos(zp(jw, :)) * wwp(jw, :)
400	v_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * sin(zp(jw, :)) * wwp(jw, :)
401	end DO
402
403	u_epflux_tot(:, ll + 1) = 0.
404	v_epflux_tot(:, ll + 1) = 0.
405	DO jw = 1, nw
406	u_epflux_tot(:, ll + 1) = u_epflux_tot(:, ll + 1) + u_epflux_p(jw, :)
407	v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :)
408	east_gwstress(:, ll) = east_gwstress(:, ll) &
409	+ max(0., u_epflux_p(jw, :)) / float(nw)
410	west_gwstress(:, ll) = west_gwstress(:, ll) &
411	+ min(0., u_epflux_p(jw, ll)) / float(nw)
412	end DO
413	end DO ! DO LL = LAUNCH, nlayer - 1
414
415	! print*,'u_epflux_tot just after launching:', u_epflux_tot
416	! print*,'v_epflux_tot just after launching:', v_epflux_tot
417	! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p)
418	! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p)
419
420	! 5. TENDENCY CALCULATIONS
421	!-------------------------
422	! 5.1 Flow rectification at the top and in the low layers
423	! --------------------------------------------------------
424	! Warning, this is the total on all GW...
425
426	u_epflux_tot(:, nlayer + 1) = 0.
427	v_epflux_tot(:, nlayer + 1) = 0.
428
429	! Here, big change compared to FLott version:
430	! We compensate (u_epflux_(:, LAUNCH), ie total emitted upward flux
431	! over the layers max(1,LAUNCH-3) to LAUNCH-1
432	DO LL = 1, max(1,LAUNCH-3)
433	u_epflux_tot(:, LL) = 0.
434	v_epflux_tot(:, LL) = 0.
435	end DO
436	DO LL = max(2,LAUNCH-2), LAUNCH-1
437	u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) &
438	+ u_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
439	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
440	v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) &
441	+ v_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
442	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
443	east_gwstress(:,LL) = east_gwstress(:, LL - 1) &
444	+ east_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
445	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
446	west_gwstress(:,LL) = west_gwstress(:, LL - 1) &
447	+ west_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) &
448	/ (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
449	end DO
450	! This way, the total flux from GW is zero, but there is a net transport
451	! (upward) that should be compensated by circulation
452	! and induce additional friction at the surface
453
454	! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
455	!---------------------------------------------
456	DO LL = 1, nlayer
457	! d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))
458	d_u(:, LL) = RG * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
459	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
460	! d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))
461	d_v(:, LL) = RG * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
462	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
463	ENDDO
464	d_t(:,:) = 0.
465
466	! 5.3 Update tendency of wind with the previous (and saved) values
467	!-----------------------------------------------------------------
468	d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) &
469	+ (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
470	d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) &
471	+ (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
472	du_nonoro_gwd(:,:) = d_u(:,:)
473	dv_nonoro_gwd(:,:) = d_v(:,:)
474
475	! print*,'u_epflux_tot just after tendency:', u_epflux_tot
476	! print*,'v_epflux_tot just after tendency:', v_epflux_tot
477	! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u)
478	! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v)
479
480	! Cosmetic: evaluation of the cumulated stress
481	ZUSTR(:) = 0.
482	ZVSTR(:) = 0.
483	DO LL = 1, nlayer
484	ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
485	ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
486	ENDDO
487
488	!#ifdef CPP_XIOS
489	! call send_xios_field("du_nonoro", d_u)
490	! call send_xios_field("dv_nonoro", d_v)
491	! call send_xios_field("east_gwstress", east_gwstress)
492	! call send_xios_field("west_gwstress", west_gwstress)
493	!#endif
494
495	END SUBROUTINE nonoro_gwd_ran
496
497	! ===================================================================
498	! Subroutines used to write variables of memory in start files
499	! ===================================================================
500
501	SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
502
503	IMPLICIT NONE
504
505	INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns
506	INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
507
508	allocate(du_nonoro_gwd(ngrid, nlayer))
509	allocate(dv_nonoro_gwd(ngrid, nlayer))
510	allocate(east_gwstress(ngrid, nlayer))
511	allocate(west_gwstress(ngrid, nlayer))
512
513	END SUBROUTINE ini_nonoro_gwd_ran
514	! ----------------------------------
515	SUBROUTINE end_nonoro_gwd_ran
516
517	IMPLICIT NONE
518
519	if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
520	if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
521	if (allocated(east_gwstress)) deallocate(east_gwstress)
522	if (allocated(west_gwstress)) deallocate(west_gwstress)
523
524	END SUBROUTINE end_nonoro_gwd_ran
525
526	END MODULE nonoro_gwd_ran_mod
527
528

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

Download in other formats:

Original Format