Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

nonoro_gwd_mix_mod.F90 @ 3935

Last change on this file since 3935 was 3935, checked in by jliu, 2 months ago
#Non-orographic gravity waves and induced turbulence #are added to the model. Further tunning are still needed.
File size: 41.9 KB

Line
1	MODULE nonoro_gwd_mix_mod
2
3	IMPLICIT NONE
4
5	REAL,allocatable,save :: du_eddymix_gwd(:,:) ! Zonal wind tendency due to GWD
6	REAL,allocatable,save :: dv_eddymix_gwd(:,:) ! Meridional wind tendency due to GWD
7	REAL,allocatable,save :: dh_eddymix_gwd(:,:) ! PT tendency due to GWD
8	REAL,allocatable,save :: dq_eddymix_gwd(:,:,:) ! tracers tendency due to GWD
9	REAL,allocatable,save :: de_eddymix_rto(:,:) ! Meridional wind tendency due to GWD
10	REAL,allocatable,save :: df_eddymix_flx(:,:) ! Meridional wind tendency due to GWD
11	!REAL,ALLOCATABLE,SAVE :: east_gwstress(:,:) ! Profile of eastward stress
12	!REAL,ALLOCATABLE,SAVE :: west_gwstress(:,:) ! Profile of westward stress
13	LOGICAL, save :: calljliu_gwimix ! flag for using non-orographic GW-induced mixing
14
15	!$OMP THREADPRIVATE(du_eddymix_gwd,dv_eddymix_gwd,dh_eddymix_gwd,dq_eddymix_gwd,de_eddymix_rto,df_eddymix_flx,calljliu_gwimix)
16	!,east_gwstress,west_gwstress)
17
18	CONTAINS
19
20	SUBROUTINE NONORO_GWD_MIX(ngrid,nlayer,DTIME,nq,cpnew, rnew, pp, &
21	zmax_therm, pt, pu, pv, pq, pht, pdt, pdu, pdv, pdq, pdht, &
22	d_pq, d_t, d_u, d_v)
23
24	!--------------------------------------------------------------------------------
25	! Parametrization of the eddy diffusion coefficient due to a discrete
26	! number of gravity waves.
27	! J.LIU
28	! Version 01, Gaussian distribution of the source
29	! LMDz model online version
30	!
31	!---------------------------------------------------------------------------------
32
33	use comcstfi_mod, only: g, pi, r,rcp
34	USE ioipsl_getin_p_mod, ONLY : getin_p
35	use vertical_layers_mod, only : presnivs
36	use geometry_mod, only: cell_area
37	use write_output_mod, only: write_output
38
39	implicit none
40
41	CHARACTER (LEN=20) :: modname='NONORO_GWD_MIX'
42	CHARACTER (LEN=80) :: abort_message
43
44
45	! 0. DECLARATIONS:
46
47	! 0.1 INPUTS
48	INTEGER, intent(in):: ngrid ! number of atmospheric columns
49	INTEGER, intent(in):: nlayer ! number of atmospheric layers
50	INTEGER, INTENT(IN) :: nq ! number of tracer species in traceurs.def
51	! integer, parameter:: nesp =42 ! number of traceurs in chemistry
52	REAL, intent(in):: DTIME ! Time step of the Physics(s)
53	REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m)
54	! REAL, intent(in):: loss(nesp) ! Chemical reaction loss rate
55	REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere
56	REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere
57	REAL, intent(in):: pp(ngrid,nlayer) ! Pressure at full levels(Pa)
58	REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels(K)
59	REAL, intent(in):: pu(ngrid,nlayer) ! Zonal wind at full levels(m/s)
60	REAL, intent(in):: pv(ngrid,nlayer) ! Meridional winds at full levels(m/s)
61	REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
62	REAL, INTENT(IN) :: pht(ngrid,nlayer) ! advected field of potential temperature
63	REAL, INTENT(IN) :: pdht(ngrid,nlayer) ! tendancy of potential temperature
64	REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq)! tendancy field pq
65	REAL,INTENT(in) :: pdt(ngrid,nlayer) ! Tendency on temperature (K/s)
66	REAL,INTENT(in) :: pdu(ngrid,nlayer) ! Tendency on zonal wind (m/s/s)
67	REAL,INTENT(in) :: pdv(ngrid,nlayer) ! Tendency on meridional wind (m/s/s)
68
69	! 0.2 OUTPUTS
70	! REAL, intent(out):: zustr(ngrid) ! Zonal surface stress
71	! REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress
72	REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature (K/s) due to gravity waves (not used set to zero)
73	REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind (m/s/s) due to gravity waves
74	REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves
75	REAL, INTENT(out) :: d_pq(ngrid,nlayer,nq)! tendancy field pq
76	REAL :: d_h(ngrid, nlayer) ! Tendency on PT (T/s/s) due to gravity waves mixing
77	! 0.3 INTERNAL ARRAYS
78	REAL :: TT(ngrid, nlayer) ! Temperature at full levels
79	REAL :: RHO(ngrid, nlayer) ! Mass density at full levels
80	REAL :: UU(ngrid, nlayer) ! Zonal wind at full levels
81	REAL :: VV(ngrid, nlayer) ! Meridional winds at full levels
82	REAL :: HH(ngrid, nlayer) ! potential temperature at full levels
83	REAL :: BVLOW(ngrid) ! initial N at each grid (not used)
84
85	INTEGER II, JJ, LL, QQ
86
87	! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
88	REAL, parameter:: DELTAT = 24. * 3600.
89
90	! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
91	INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers
92	INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed
93	INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed
94	INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed
95	INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves
96	INTEGER JK, JP, JO, JW ! Loop index for NK,NP,NO, and NW
97
98	REAL, save :: kmax ! Max horizontal wavenumber (lambda min,lambda=2*pi/kmax),kmax=N/u, u(=30~50) zonal wind velocity at launch altitude
99	!$OMP THREADPRIVATE(kmax)
100	REAL, save :: kmin ! Min horizontal wavenumber (lambda max = 314 km,lambda=2*pi/kmin)
101	!$OMP THREADPRIVATE(kmin)
102	REAL kstar ! Min horizontal wavenumber constrained by horizontal resolution
103
104	REAL :: max_k(ngrid) ! max_k=max(kstar,kmin)
105	REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity (not used)
106	REAL CPHA(ngrid) ! absolute PHASE VELOCITY frequency
107	REAL ZK(NW, ngrid) ! Horizontal wavenumber amplitude
108	REAL ZP(NW, ngrid) ! Horizontal wavenumber angle
109	REAL ZO(NW, ngrid) ! Absolute frequency
110
111	REAL maxp(NW,ngrid) ! wave saturation index
112	REAL maxs(NW,ngrid) ! wave saturation index
113	integer LLSATURATION(NW,ngrid) ! layer where the gravity waves break
114	integer LLZCRITICAL(NW,ngrid) ! layer where the gravity waves depleting
115	integer ZHSATURATION(NW,ngrid) ! altitude of the layer where the gravity waves break
116	INTEGER ll_zb(ngrid) ! layer where the gravity waves break
117	INTEGER ll_zc(ngrid) ! layer where the gravity waves break
118	integer ll_zb_ii,ll_zc_ii
119	integer ll_zb_max, ll_zb_max_r
120	REAL d_eddy_mix_p(NW,ngrid) ! Diffusion coefficients where ll> = ll_zb
121	REAL d_eddy_mix_m(NW,ngrid) ! Diffusion coefficients where ll< ll_zb
122	REAL d_eddy_mix_s(NW,ngrid) ! Diffusion coefficients where ll< ll_zb
123	REAL lambda_img(NW,ngrid)
124	REAL d_eddy_mix_p_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll> = ll_zb
125	REAL d_eddy_mix_m_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll< ll_zb
126	REAL d_eddy_mix_tot_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll> = ll_zb
127	REAL d_eddy_mix_tot(ngrid, nlayer+1)
128	REAL d_eddy_mix(NW,ngrid) ! Comprehensive Diffusion coefficients
129	REAL d_wave(NW,ngrid) ! coefficients consider the tracers' gradients
130	REAL u_eddy_mix_p(NW, ngrid) ! Zonal Diffusion coefficients
131	REAL v_eddy_mix_p(NW, ngrid) ! Meridional Diffusion coefficients
132	REAL h_eddy_mix_p(NW, ngrid) ! potential temperature DC
133	Real u_eddy_mix_tot(ngrid, nlayer+1) ! Total zonal D
134	Real v_eddy_mix_tot(ngrid, nlayer+1) ! Total meridional D
135	Real h_eddy_mix_tot(ngrid, nlayer+1) ! Total PT D
136	REAL U_shear(ngrid,nlayer)
137	Real wwp_vertical_tot(nlayer+1, NW, ngrid) ! Total meridional D
138	Real wwp_vertical_ll(nlayer+1)
139	real eddy_mix_ll(nlayer)
140	real eddy_mix_tot_ll(nlayer)
141	REAL pq_eddy_mix_p(NW,ngrid,nq)
142	REAL pq_eddy_mix_tot(ngrid, nlayer+1,nq)
143	REAL zq(ngrid,nlayer,nq) ! advected field nq
144	REAL zq_var(ngrid,nlayer,nq)
145	REAL dzq_var(ngrid,nlayer,nq)
146	REAL mdzq_var(nlayer,nq)
147	REAL zq_ave(nlayer,nq)
148	REAL dzq_ave(nlayer,nq)
149	REAL zq_ratio(nlayer,nq)
150	REAL, save:: eff
151	!$OMP THREADPRIVATE(eff)
152	REAL, save:: eff1
153	!$OMP THREADPRIVATE(eff1)
154	REAL, save:: vdl
155	!$OMP THREADPRIVATE(vdl)
156
157
158	REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
159	REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
160	REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level
161	REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level
162	REAL wwpsat(nw,ngrid) ! Wave EP-fluxes of saturation
163	REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
164	REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
165	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)
166	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)
167	REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0)
168	REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable)
169	!$OMP THREADPRIVATE(epflux_max)
170	INTEGER LAUNCH ! Launching altitude
171	REAL, save :: xlaunch ! Control the launching altitude by pressure
172	!$OMP THREADPRIVATE(xlaunch)
173	REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx. not used)
174	REAL cmax(ngrid,nlayer) ! Max horizontal absolute phase velocity (at the maginitide of zonal wind u at the launch altitude)
175
176	! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
177	REAL, save :: sat ! saturation parameter(tunable)
178	!$OMP THREADPRIVATE(sat)
179	REAL, save :: rdiss ! dissipation coefficient (tunable)
180	!$OMP THREADPRIVATE(rdiss)
181	REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic frequency
182
183	! 0.3.3 Background flow at 1/2 levels and vertical coordinate
184	REAL H0bis(ngrid, nlayer) ! Characteristic Height of the atmosphere (specific locations)
185	REAL, save:: H0 ! Characteristic Height of the atmosphere (constant)
186	!$OMP THREADPRIVATE(H0)
187	REAL, parameter:: pr = 250 ! Reference pressure [Pa]
188	REAL, parameter:: tr = 220. ! Reference temperature [K]
189	REAL ZH(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0)
190	REAL ZHbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H0bis)
191	REAL UH(ngrid, nlayer + 1) ! zonal wind at 1/2 levels
192	REAL VH(ngrid, nlayer + 1) ! meridional wind at 1/2 levels
193	REAL PH(ngrid, nlayer + 1) ! Pressure at 1/2 levels
194	REAL, parameter:: psec = 1.e-20 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure)
195	REAL BV(ngrid, nlayer + 1) ! Brunt Vaisala freq. (N^2) at 1/2 levels
196	REAL, parameter:: bvsec = 1.e-3 ! Security to avoid negative BV
197	REAL HREF(nlayer + 1) ! Reference pressure for launching altitude
198
199
200	REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
201	INTEGER first_breaking_flag(ngrid)
202	INTEGER first_satuatio_flag(ngrid)
203
204	LOGICAL,SAVE :: firstcall = .true.
205	!$OMP THREADPRIVATE(firstcall)
206
207	!-----------------------------------------------------------------------------------------------------------------------
208	! 1. INITIALISATIONS
209	!-----------------------------------------------------------------------------------------------------------------------
210	IF (firstcall) THEN
211	write(,) "nonoro_gwd_mix: non-oro GW-induced mixing is active!"
212	epflux_max = 5.E-4 ! Mars' value !!
213	call getin_p("nonoro_gwd_epflux_max", epflux_max)
214	write(,) "NONORO_GWD_MIX: epflux_max=", epflux_max
215	sat = 1.5 ! default gravity waves saturation value !!
216	call getin_p("nonoro_gwd_sat", sat)
217	write(,) "NONORO_GWD_MIX: sat=", sat
218	! cmax = 50. ! default gravity waves phase velocity value !!
219	! call getin_p("nonoro_gwd_cmax", cmax)
220	! write(,) "NONORO_GWD_MIX: cmax=", cmax
221	rdiss=0.07 ! default coefficient of dissipation !!
222	call getin_p("nonoro_gwd_rdiss", rdiss)
223	write(,) "NONORO_GWD_MIX: rdiss=", rdiss
224	kmax=1.E-4 ! default Max horizontal wavenumber !!
225	call getin_p("nonoro_gwd_kmax", kmax)
226	write(,) "NONORO_GWD_MIX: kmax=", kmax
227	kmin=7.E-6 ! default Max horizontal wavenumber !!
228	call getin_p("nonoro_gwd_kmin", kmin)
229	write(,) "NONORO_GWD_MIX: kmin=", kmin
230	xlaunch=0.6 ! default GW luanch altitude controller !!
231	call getin_p("nonoro_gwd_xlaunch", xlaunch)
232	write(,) "NONORO_GWD_MIX: xlaunch=", xlaunch
233	eff=0.1 ! Diffusion effective factor !!
234	call getin_p("nonoro_gwimixing_eff", eff)
235	write(,) "NONORO_GWD_MIX: eff=", eff
236	eff1=0.1 ! Diffusion effective factor !!
237	call getin_p("nonoro_gwimixing_eff1", eff1)
238	write(,) "NONORO_GWD_MIX: eff1=", eff1
239	vdl=1.5 ! Diffusion effective factor !!
240	call getin_p("nonoro_gwimixing_vdl", vdl)
241	write(,) "NONORO_GWD_MIX: vdl=", vdl
242	! Characteristic vertical scale height
243	H0 = r * tr / g
244	! Control
245	if (deltat .LT. dtime) THEN
246	call abort_physic("NONORO_GWD_MIX","gwd random: deltat lower than dtime!",1)
247	endif
248	if (nlayer .LT. nw) THEN
249	call abort_physic("NONORO_GWD_MIX","gwd random: nlayer lower than nw!",1)
250	endif
251	firstcall = .false.
252	ENDIF
253
254	! Compute current values of temperature and winds
255	tt(:,:)=pt(:,:)+dtime*pdt(:,:)
256	uu(:,:)=pu(:,:)+dtime*pdu(:,:)
257	vv(:,:)=pv(:,:)+dtime*pdv(:,:)
258	zq(:,:,:)=pq(:,:,:)+dtime*pdq(:,:,:)
259	hh(:,:)=pht(:,:)+dtime*pdht(:,:)
260
261	! tracer average and gradients for mixing
262	zq_ave(:,:)=0.
263	zq_ave(:,:) = SUM(zq(:,:,:), dim=1) / real(ngrid)
264
265	zq_var(:,:,:)=0.
266	DO ii=1,ngrid
267	zq_var(ii,:,:)=zq(ii,:,:)-zq_ave(:,:)
268	endDO
269
270	dzq_var(:,:,:)=0.
271	mdzq_var(:,:) =0.
272	dzq_ave(:,:) =0.
273	zq_ratio(:,:) =0.
274	DO LL=1,nlayer-1
275	dzq_var(:, LL, :) = (zq_var(:, LL+1, :) - zq_var(:, LL, :))**2.
276	mdzq_var(LL, :) = SUM(dzq_var(:, LL, :), dim=1)/real(ngrid)
277	dzq_ave(LL, :) = (zq_ave(LL+1, :) - zq_ave(LL, :))**2.
278	zq_ratio(LL,:) = MAX(0., (mdzq_var(LL+1,:) + mdzq_var(LL,:)))&
279	/MAX(1E-15, (dzq_ave(LL+1,:) + dzq_ave(LL,:)))
280	! do qq=1,nq
281	! print*, 'ratio=', zq_ratio(LL,QQ)
282	! enddo
283	endDO
284	dzq_var(:,nlayer,:) = dzq_var(:,nlayer-1, :)
285	mdzq_var(nlayer,:) = mdzq_var(nlayer-1, :)
286	dzq_ave(nlayer,:) = dzq_ave(nlayer-1, :)
287	! Compute the real mass density by rho=p/R(T)T
288	DO ll=1,nlayer
289	DO ii=1,ngrid
290	rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll))
291	ENDDO
292	ENDDO
293	! print*,'epflux_max just after firstcall:', epflux_max
294
295	!-----------------------------------------------------------------------------------------------------------------------
296	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
297	!-----------------------------------------------------------------------------------------------------------------------
298	! Pressure and inverse of pressure at 1/2 level
299	DO LL = 2, nlayer
300	PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
301	end DO
302	PH(:, nlayer + 1) = 0.
303	PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
304
305	! call write_output('nonoro_pp','nonoro_pp', 'm',PP(:,:))
306	! call write_output('nonoro_ph','nonoro_ph', 'm',PH(:,:))
307
308	! Launching level for reproductible case
309	!Pour revenir a la version non reproductible en changeant le nombre de
310	!process
311	! Reprend la formule qui calcule PH en fonction de PP=play
312	DO LL = 2, nlayer
313	HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
314	end DO
315	HREF(nlayer + 1) = 0.
316	HREF(1) = 2. * presnivs(1) - HREF(2)
317
318	LAUNCH=0
319	DO LL = 1, nlayer
320	IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
321	ENDDO
322
323	! Log pressure vert. coordinate
324	ZH(:,1) = 0.
325	DO LL = 2, nlayer + 1
326	!ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
327	H0bis(:, LL-1) = r * TT(:, LL-1) / g
328	ZH(:, LL) = ZH(:, LL-1) &
329	+ H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
330	end DO
331	ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC))
332
333	! call write_output('nonoro_zh','nonoro_zh', 'm',ZH(:,2:nlayer+1))
334
335	! Winds and Brunt Vaisala frequency
336	DO LL = 2, nlayer
337	UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
338	VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
339	! Brunt Vaisala frequency (=g/T*[dT/dz + g/cp] )
340	BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1))) &
341	*((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+ &
342	G / cpnew(:,LL))
343
344	BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive
345	BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small
346	end DO
347	BV(:, 1) = BV(:, 2)
348	UH(:, 1) = 0.
349	VH(:, 1) = 0.
350	BV(:, nlayer + 1) = BV(:, nlayer)
351	UH(:, nlayer + 1) = UU(:, nlayer)
352	VH(:, nlayer + 1) = VV(:, nlayer)
353	cmax(:,launch)=UU(:, launch)
354	DO ii=1,ngrid
355	KSTAR = PI/SQRT(cell_area(II))
356	MAX_K(II)=MAX(kmin,kstar)
357	ENDDO
358	call write_output('nonoro_bv','Brunt Vaisala frequency in nonoro', 'Hz',BV(:,2:nlayer+1))
359
360	!-----------------------------------------------------------------------------------------------------------------------
361	! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
362	!-----------------------------------------------------------------------------------------------------------------------
363	! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way
364	DO JW = 1, NW
365	! Angle
366	DO II = 1, ngrid
367	! Angle (0 or PI so far)
368	RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
369	RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
370	ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.)* PI / 2.
371
372	! Horizontal wavenumber amplitude
373	! From Venus model: TN+GG April/June 2020 (rev2381)
374	! "Individual waves are not supposed to occupy the entire domain, but only a fraction of it" (Lott+2012)
375	! ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
376	KSTAR = PI/SQRT(cell_area(II))
377	ZK(JW, II) = MAX_K(II) + (KMAX - MAX_K(II)) *RAN_NUM_2
378
379	! Horizontal phase speed
380	! this computation allows to get a gaussian distribution centered on 0 (right ?)
381	! then cmin is not useful, and it favors small values.
382	CPHA(:) = 0.
383	DO JJ = 1, NA
384	RAN_NUM_3=MOD(TT(II, JW+3JJ)*2, 1.)
385	CPHA(ii) = CPHA(ii) + 2.CMAX(ii,launch) &
386	(RAN_NUM_3 -0.5)* &
387	SQRT(3.)/SQRT(NA*1.)
388	END DO
389	IF (CPHA(ii).LT.0.) THEN
390	CPHA(ii) = -1.*CPHA(ii)
391	ZP(JW,II) = ZP(JW,II) + PI
392	ENDIF
393	! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
394	! ran_num_3 = mod(tt(ii, jw)**2, 1.)
395	! cpha = cmin + (cmax - cmin) * ran_num_3
396
397	! Intrinsic frequency
398	ZO(JW, II) = CPHA(II) * ZK(JW, II)
399	! Intrinsic frequency is imposed
400	ZO(JW, II) = ZO(JW, II) &
401	+ ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
402	+ ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
403
404	! Momentum flux at launch level
405	epflux_0(JW, II) = epflux_max &
406	* MOD(100. * (UU(II, JW)2 + VV(II, JW)2), 1.)
407	ENDDO
408	end DO
409	! print*,'epflux_0 just after waves charac. ramdon:', epflux_0
410
411	!-----------------------------------------------------------------------------------------------------------------------
412	! 4. COMPUTE THE FLUXES
413	!-----------------------------------------------------------------------------------------------------------------------
414	! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes.
415	!------------------------------------------------------
416	DO JW = 1, NW
417	! Evaluate intrinsic frequency at launching altitude:
418	intr_freq_p(JW, :) = ZO(JW, :) &
419	- ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
420	- ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
421	end DO
422
423	! VERSION WITHOUT CONVECTIVE SOURCE
424	! Vertical velocity at launch level, value to ensure the imposed
425	! mom flux:
426	DO JW = 1, NW
427	! WW is directly a flux, here, not vertical velocity anymore
428	WWP(JW, :) = epflux_0(JW,:)
429	u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
430	v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
431	end DO
432
433	! 4.2 Initial flux at launching altitude
434	!------------------------------------------------------
435	u_epflux_tot(:, LAUNCH) = 0
436	v_epflux_tot(:, LAUNCH) = 0
437	DO JW = 1, NW
438	u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :)
439	v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :)
440	end DO
441
442	! 4.3 Loop over altitudes, with passage from one level to the next done by:
443	!----------------------------------------------------------------------------
444	! i) conserving the EP flux,
445	! ii) dissipating a little,
446	! iii) testing critical levels,
447	! iv) testing the breaking.
448	!----------------------------------------------------------------------------
449
450	wwp_vertical_tot(:, :,:) =0.
451	DO LL = LAUNCH, nlayer - 1
452	! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL TO THE NEXT)
453	DO JW = 1, NW
454	intr_freq_m(JW, :) = intr_freq_p(JW, :)
455	WWM(JW, :) = WWP(JW, :)
456	! Intrinsic Frequency
457	intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) &
458	- ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1)
459	! Vertical velocity in flux formulation
460
461	wwpsat(JW,:) = ABS(intr_freq_p(JW, :))*3 / (2.BV(:, LL+1)) &
462	* rho(:,launch)*exp(-zh(:, ll + 1) / H0) &
463	* SAT*2 KMIN2 / ZK(JW, :)4
464	WWP(JW, :) = MIN( &
465	! No breaking (Eq.6)
466	WWM(JW, :) &
467	! Dissipation (Eq. 8)(real rho used here rather than pressure
468	! parametration because the original code has a bug if the density of
469	! the planet at the launch altitude not approximate 1): !
470	* EXP(- RDISS*2./rho(:, LL + 1) &
471	* ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
472	/ MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 &
473	* ZK(JW, :)*3 (ZH(:, LL + 1) - ZH(:, LL))), &
474	! Critical levels (forced to zero if intrinsic frequency
475	! changes sign)
476	MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) &
477	! Saturation (Eq. 12) (rho at launch altitude is imposed by J.Liu.
478	! Same reason with Eq. 8)
479	* WWPSAT(JW,:))
480	end DO
481	! END OF W(KB)ARNING
482
483	! first_breaking_flag(:)=0 !mixing start at first breaking
484	! first_satuatio_flag(:)=0
485
486	DO JW=1,NW
487	wwp_vertical_tot(ll, JW, :) = WWP(JW,:)/rho(:,ll)
488	ENDDO
489	end DO ! DO LL = LAUNCH, nlayer - 1
490	! print*, "this line is for tunning"
491
492	DO II=1, ngrid
493	DO JW=1,NW
494	wwp_vertical_ll(:)=wwp_vertical_tot(:, JW, II)
495	ll_zb_max = MAXloc(wwp_vertical_ll(:),1)
496	!if (LLSATURATION(JW,II).ne.-1000) then
497	LLSATURATION(JW,II)=ll_zb_max
498	!endif
499	! print*, "this line is for tunning"
500	if (ll_zb_max .gt. 1) then
501	!ll_zb_max_r =MAXloc(wwp_vertical_ll(nlayer:1:-1),1)
502	!if (ll_zb_max_r .gt. ll_zb_max) LLSATURATION(JW,II)=ll_zb_max_r
503	DO ll = nlayer,ll_zb_max,-1
504	if (wwp_vertical_ll(ll)-wwp_vertical_ll(ll-1).gt. 0) then
505	LLSATURATION(JW,II)=ll
506	goto 119
507	endif
508	ENDDO
509	119 continue
510	endif
511	! if (ll_zb_max .gt. launch) then
512	! DO LL = (nlayer + 1), ll_zb_max,-1
513	! if (abs(wwp_vertical_ll(ll)).le. 1.e-9) LLZCRITICAL(JW,II)=LL
514	! ENDDO
515	! else
516	! LLZCRITICAL(JW,II)=1
517	! endif
518	!print*, "this line is for tunning"
519	ENDDO
520	ENDDO
521	!print*, "this line is for tunning"
522
523
524	DO JW = 1, NW
525	! Evaluate intrinsic frequency at launching altitude:
526	intr_freq_p(JW, :) = ZO(JW, :) &
527	- ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
528	- ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
529	end DO
530	d_eddy_mix_p_ll(:,:,:)=0.
531	d_eddy_mix_m_ll(:,:,:)=0.
532	d_eddy_mix_p(:,:)=0.
533	d_eddy_mix_m(:,:)=0.
534	d_eddy_mix_s(:,:)=0.
535	lambda_img(:,:) = 0.
536	u_shear(:,:)= 0.
537	DO LL = LAUNCH, nlayer - 1
538	U_shear(:,ll)=(UU(:, LL+1)-UU(:, LL)) /(ZH(:, LL+1) - ZH(:, LL))
539	DO II=1,ngrid
540	! all the eddy diffusion parameters are culculated at here
541	DO JW = 1, NW
542	intr_freq_m(JW, II) = intr_freq_p(JW, II)
543	! Intrinsic Frequency
544	intr_freq_p(JW, II) = ZO(JW, II) - ZK(JW, II) * COS(ZP(JW,II)) * UH(II, LL + 1) &
545	- ZK(JW, II) * SIN(ZP(JW,II)) * VH(II, LL + 1)
546	ll_zb(II)= LLSATURATION(JW,II)
547	ll_zb_ii = ll_zb(II)
548	! ll_zc(II)= LLZCRITICAL(JW,II)
549	! ll_zc_ii = ll_zc(II)
550	! If (ll_zb_ii.le. launch .or. ll .gt. ll_zc_ii) then
551	If (ll .lt. ll_zb_ii .or. ll_zb_ii.lt. launch) then
552	d_eddy_mix_p(JW,II)=0.
553	d_eddy_mix_m(JW,II)=0.
554	d_eddy_mix_p_ll(ll,JW,ii)=0.
555	d_eddy_mix_m_ll(ll,JW,ii)=0.
556	endif
557	IF (LL.GE.ll_zb_ii .and. ll_zb_ii.ge. launch) THEN
558	lambda_img(JW,II) = 0.5 / H0bis(II,LL) &
559	- 1.5 /((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.) &
560	*(intr_freq_p(JW, II) - intr_freq_m(JW, II)) &
561	/(ZH(II, LL+1) - ZH(II, LL)) &
562	+ 1.5/((BV(II,LL)+BV(II, LL+1))/2.) &
563	*(BV(II, LL+1)-BV(II, LL)) /(ZH(II, LL+1) - ZH(II, LL))
564	!lambda_img(JW,II) = max(0.5 / H0bis(II,LL), abs(lambda_img(JW,II)))
565	if (lambda_img(JW,II).lt. 0) lambda_img(JW,II) = 0.
566	!if (lambda_img(JW,II).gt. 0.5/H0bis(II,LL)) lambda_img(JW,II) = 0.5/H0bis(II,LL)
567	d_eddy_mix_s(JW,II) = eff*MAX(ABS(intr_freq_p(JW, II) + intr_freq_m(JW, II)) / 2., &
568	ZOISEC)4 / (ZK(JW, II)3 * ((BV(II,LL)+BV(II, LL+1))/2.)**3)&
569	*lambda_img(JW,II)
570	! abs(0.5 / H0bis(II,LL) - 1.5ZK(JW, II) &
571	! /abs((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.) &
572	! *(UU(II, LL+1)-UU(II, LL)) /(ZH(II, LL+1) - ZH(II, LL)) &
573	! + 1.5/((BV(II,LL)+BV(II, LL+1))/2.) &
574	! *(BV(II, LL+1)-BV(II, LL)) /(ZH(II, LL+1) - ZH(II, LL)))
575	! *MAX(0., sign(1., lambda_img(JW,II)))
576
577	d_eddy_mix_p(JW,II) = Min( d_eddy_mix_s(JW,II) , &
578	-((wwp_vertical_tot(ll+1, JW, II)-wwp_vertical_tot(ll, JW, II)) &
579	/(ZH(II, LL+1) - ZH(II, LL))) &
580	((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.)eff1 &
581	/(((BV(II, LL+1) + BV(II,LL)) / 2.)*2 ZK(JW, II)))
582
583	if (d_eddy_mix_p(jw,ii) .lt. 0.) then
584	print*, "this line is for tunning"
585	endif
586	ENDIF
587
588
589	end DO !JW = 1, NW
590	end DO !II=1,ngrid
591	! print*, "this line is for tunning"
592	d_eddy_mix_p_ll(ll,:,:)=d_eddy_mix_p(:,:)
593	! print*, "this line is for tunning"
594	! if (ll.ge.55) then
595	! print*, "this line is for tunning"
596	! endif
597	end DO ! DO LL = LAUNCH, nlayer - 1
598
599	call write_output('zonal_shear','u shear', 's-1',u_shear(:,:))
600
601	DO II=1,ngrid
602	DO JW = 1, NW
603	ll_zb(II)= LLSATURATION(JW,II)
604	ll_zb_ii = ll_zb(II)
605
606	do LL=launch, nlayer-1
607	IF (LL.eq.ll_zb_ii .and. ll_zb_ii .gt.launch) THEN
608	d_eddy_mix_m(JW,II) = d_eddy_mix_p_ll(ll,JW,ii)
609	! if (ii.eq. 848) then
610	! print*, "this line is for tunning"
611	! endif
612	ENDIF
613	enddo
614
615	ENDDO
616	ENDDO
617
618	DO II=1,ngrid
619	DO JW = 1, NW
620	ll_zb(II)= LLSATURATION(JW,II)
621	ll_zb_ii = ll_zb(II)
622	DO LL= launch, (ll_zb_ii - 1)
623	if (ll_zb_ii .gt. launch) then
624	d_eddy_mix_m_ll(ll,JW,II) = d_eddy_mix_m(JW,II) &
625	* exp(vdl*(ZH(II, LL)-ZH(II, ll_zb_ii) ) / H0)
626	else
627	d_eddy_mix_m_ll(ll,JW,II) = 0.
628	endif
629	endDO
630	ENDDO
631	ENDDO
632	! print*, "this line is for tunning"
633
634
635	DO II=1, ngrid
636	DO JW=1,NW
637	eddy_mix_ll(:) = d_eddy_mix_p_ll(:,JW,II)+ d_eddy_mix_m_ll(:,JW,II)
638	! print*, "this line is for tunning"
639	enddo
640	ENDDO
641
642
643	DO JW = 1, NW
644	! Evaluate intrinsic frequency at launching altitude:
645	intr_freq_p(JW, :) = ZO(JW, :) &
646	- ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
647	- ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
648	end DO
649	u_eddy_mix_tot(:, :) = 0.
650	v_eddy_mix_tot(:, :) = 0.
651	h_eddy_mix_tot(:, :) = 0.
652	u_eddy_mix_p(:, :)=0.
653	v_eddy_mix_p(:, :)=0.
654	h_eddy_mix_p(:, :)=0.
655	d_eddy_mix_tot_ll(:,:,:)=0.
656	pq_eddy_mix_p(:,:,:)=0.
657	pq_eddy_mix_tot(:, :,:)=0.
658	d_eddy_mix(:,:)=0.
659	d_wave(:, :) =0.
660	d_eddy_mix_p_ll(nlayer,:,:)=d_eddy_mix_p_ll(nlayer-1,:,:)
661	d_eddy_mix_tot(:, :) =0.
662	DO LL = LAUNCH, nlayer - 1
663	d_eddy_mix(:,:) = d_eddy_mix_m_ll(ll,:,:) + d_eddy_mix_p_ll(ll,:,:)
664	DO JW = 1, NW
665	u_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(UU(:, LL + 1) - UU(:, LL)) &
666	/(ZH(:, LL + 1) - ZH(:, LL)) &
667	SIGN(1.,intr_freq_p(JW, :)) COS(ZP(JW, :))
668	v_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(VV(:, LL + 1) - VV(:, LL)) &
669	/(ZH(:, LL + 1) - ZH(:, LL)) &
670	SIGN(1.,intr_freq_p(JW, :)) SIN(ZP(JW, :))
671	h_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(HH(:, LL + 1) - HH(:, LL)) &
672	/(ZH(:, LL + 1) - ZH(:, LL)) &
673	SIGN(1.,intr_freq_p(JW, :)) SIN(ZP(JW, :))
674
675	ENDDO
676	u_eddy_mix_tot(:, LL+1)=0.
677	v_eddy_mix_tot(:, LL+1)=0.
678	h_eddy_mix_tot(:, LL+1)=0.
679	DO JW=1,NW
680	u_eddy_mix_tot(:, LL+1) = u_eddy_mix_tot(:, LL+1) + u_eddy_mix_p(JW, :)
681	v_eddy_mix_tot(:, LL+1) = v_eddy_mix_tot(:, LL+1) + v_eddy_mix_p(JW, :)
682	h_eddy_mix_tot(:, LL+1) = h_eddy_mix_tot(:, LL+1) + h_eddy_mix_p(JW, :)
683	ENDDO
684	DO JW=1,NW
685	! d_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL+1) + d_eddy_mix(JW,:)
686	d_eddy_mix_tot(:, LL+1) = d_eddy_mix_tot(:, LL+1) + d_eddy_mix(JW,:)
687	! u_eddy_mix_tot(:, LL+1) = u_eddy_mix_tot(:, LL+1)+ u_eddy_mix_p(JW, :)
688	! v_eddy_mix_tot(:, LL+1) = v_eddy_mix_tot(:, LL+1)+ v_eddy_mix_p(JW, :)
689	ENDDO
690	! if (ll.ge.55) then
691	! print*, "this line is for tunning"
692	!endif
693	! d_wave(JW,:) = d_eddy_mix(JW, :)(mdzq_var(LL,QQ)/dzq_ave(LL,QQ))*2.
694	DO QQ=1,NQ
695	DO JW=1,NW
696	d_wave(JW, :) = d_eddy_mix(JW, :)*zq_ratio(LL,QQ)
697	pq_eddy_mix_p(JW, :, QQ) = d_wave(JW, :)* (zq(:, LL + 1,QQ)- zq(:, LL, QQ)) &
698	/(ZH(:, LL + 1)- ZH(:, LL)) &
699	SIGN(1.,intr_freq_p(JW, :)) SIN(ZP(JW, :))
700	! pq_eddy_mix_tot(:, LL+1,QQ) = pq_eddy_mix_tot(:, LL+1,QQ) &
701	! + pq_eddy_mix_p(JW, :, QQ)
702	ENDDO
703	ENDDO
704	pq_eddy_mix_tot(:, LL+1,:)=0.
705	DO JW=1,NW
706	DO QQ=1, NQ
707	pq_eddy_mix_tot(:, LL+1,QQ) = pq_eddy_mix_tot(:, LL+1,QQ) + pq_eddy_mix_p(JW, :, QQ)
708	ENDDO
709	endDO
710	! print*, "this line is for tunning"
711	ENDDO !LL = LAUNCH, nlayer - 1
712
713	d_eddy_mix_tot(:, LAUNCH) = d_eddy_mix_tot(:, LAUNCH+1)
714	d_eddy_mix_tot(:, nlayer + 1) = d_eddy_mix_tot(:, nlayer)
715	d_eddy_mix_tot(:,:) = DTIME/DELTAT/REAL(NW) * d_eddy_mix_tot(:,:) &
716	+ (1.-DTIME/DELTAT) * de_eddymix_rto(:,:)
717	call write_output('nonoro_d_mixing_tot','Total EP Flux along U in nonoro', 'm2s-1',d_eddy_mix_tot(:,2:nlayer+1))
718	! u_eddy_mix_tot(:, :) = 0.
719	! v_eddy_mix_tot(:, :) = 0.
720	! pq_eddy_mix_tot(:, :, :) = 0.
721	! DO LL = 1, nlayer-1
722	!u_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL)*(UU(:, LL + 1) - UU(:, LL)) &
723	! /(ZH(:, LL + 1) - ZH(:, LL)) !*sign(1.,u_shear(:, LL))
724	!v_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL)*(VV(:, LL + 1) - VV(:, LL)) &
725	!
726	! /(ZH(:, LL + 1) - ZH(:, LL)) !*sign(1.,u_shear(:, LL))
727
728
729
730	! DO QQ=1, NQ
731	! pq_eddy_mix_tot(:, LL,QQ) = d_eddy_mix_tot(:, LL) &
732	! * (zq(:, LL + 1,QQ)- zq(:, LL, QQ)) &
733	! /(ZH(:, LL + 1)- ZH(:, LL))
734	! ENDDO
735
736
737
738
739	!ENDDO !LL = LAUNCH, nlayer - 1
740
741	de_eddymix_rto(:,:) = d_eddy_mix_tot(:,:)
742	!-----------------------------------------------------------------------------------------------------------------------
743	! 5. TENDENCY CALCULATIONS
744	!-----------------------------------------------------------------------------------------------------------------------
745
746	! 5.1 Flow rectification at the top and in the low layers
747	! --------------------------------------------------------
748	! Warning, this is the total on all GW
749	u_eddy_mix_tot(:, nlayer + 1) = 0.
750	v_eddy_mix_tot(:, nlayer + 1) = 0.
751	h_eddy_mix_tot(:, nlayer + 1) = 0.
752	pq_eddy_mix_tot(:, nlayer + 1,:)=0.
753	! Here, big change compared to FLott version:
754	! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
755	! over the layers max(1,LAUNCH-3) to LAUNCH-1
756	DO LL = 1, max(1,LAUNCH-3)
757	u_eddy_mix_tot(:, LL) = 0.
758	v_eddy_mix_tot(:, LL) = 0.
759	h_eddy_mix_tot(:, LL) = 0.
760	end DO
761
762	DO QQ=1,NQ
763	DO LL = 1, max(1,LAUNCH-3)
764	pq_eddy_mix_tot(:, LL,QQ) = 0.
765	end DO
766	ENDDO
767
768	DO LL = max(2,LAUNCH-2), LAUNCH-1
769	u_eddy_mix_tot(:, LL) = u_eddy_mix_tot(:, LL - 1) + u_eddy_mix_tot(:, LAUNCH) &
770	* (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
771	v_eddy_mix_tot(:, LL) = v_eddy_mix_tot(:, LL - 1) + v_eddy_mix_tot(:, LAUNCH) &
772	* (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
773	h_eddy_mix_tot(:, LL) = h_eddy_mix_tot(:, LL - 1) + h_eddy_mix_tot(:, LAUNCH) &
774	* (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
775	end DO
776
777
778	DO QQ=1,NQ
779	DO LL = max(2,LAUNCH-2), LAUNCH-1
780	pq_eddy_mix_tot(:, LL,QQ) = pq_eddy_mix_tot(:, LL-1,QQ)+pq_eddy_mix_tot(:, LAUNCH,QQ)&
781	* (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
782	end DO
783	ENDDO
784	!u_eddy_mix_tot(:,:) = DTIME/DELTAT/REAL(NW) * u_eddy_mix_tot(:,:) &
785	! + (1.-DTIME/DELTAT) * df_eddymix_flx(:,:)
786
787	!do ii=1,ngrid
788	! if (u_eddy_mix_tot(ii,58).lt. -8000.) then
789	! print*, ii
790	! endif
791	!enddo
792
793	! This way, the total flux from GW is zero, but there is a net transport
794	! (upward) that should be compensated by circulation
795	! and induce additional friction at the surface
796	call write_output('nonoro_u_mixing_tot','Total EP Flux along U in nonoro', '',u_eddy_mix_tot(:,2:nlayer+1))
797	call write_output('nonoro_v_mixing_tot','Total EP Flux along V in nonoro', '',v_eddy_mix_tot(:,2:nlayer+1))
798
799	! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
800	!---------------------------------------------
801	DO LL = 1, nlayer
802	!Notice here the d_u and d_v are tendency (For Mars) but not increment(Venus).
803	d_u(:, LL) = (u_eddy_mix_tot(:, LL + 1) - u_eddy_mix_tot(:, LL)) &
804	/ (ZH(:, LL + 1) - ZH(:, LL))
805	!d_v(:, LL) = (v_eddy_mix_tot(:, LL + 1) - v_eddy_mix_tot(:, LL)) &
806	! / (ZH(:, LL + 1) - ZH(:, LL))
807	d_h(:, LL) = (h_eddy_mix_tot(:, LL + 1) - h_eddy_mix_tot(:, LL)) &
808	/ (ZH(:, LL + 1) - ZH(:, LL))
809	ENDDO
810
811	DO QQ=1,NQ
812	DO ll = 1, nlayer
813	d_pq(:, ll, QQ) = (pq_eddy_mix_tot(:, LL + 1,QQ) - pq_eddy_mix_tot(:, LL, QQ)) &
814	/ (ZH(:, LL + 1) - ZH(:, LL))
815	end DO
816	ENDDO
817	!df_eddymix_flx(:,:) = u_eddy_mix_tot(:,:)
818	!d_pq(:, :, :)=0.
819	!d_t(:,:) = 0.
820	!d_v(:,:) = 0.
821	!zustr(:) = 0.
822	!zvstr(:) = 0.
823	! call write_output('nonoro_d_u','nonoro_d_u', '',d_u(:,:))
824	! call write_output('nonoro_d_v','nonoro_d_v', '',d_v(:,:))
825
826	! 5.3 Update tendency of wind with the previous (and saved) values
827	!-----------------------------------------------------------------
828	d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) &
829	+ (1.-DTIME/DELTAT) * du_eddymix_gwd(:,:)
830	!d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) &
831	! + (1.-DTIME/DELTAT) * dv_eddymix_gwd(:,:)
832	du_eddymix_gwd(:,:) = d_u(:,:)
833	!dv_eddymix_gwd(:,:) = d_v(:,:)
834	d_v(:,:)=0.
835	call write_output('du_eddymix_gwd','Tendency on U due to nonoro GW', 'm.s-2',du_eddymix_gwd(:,:))
836	!call write_output('dv_eddymix_gwd','Tendency on V due to nonoro GW', 'm.s-2',dv_eddymix_gwd(:,:))
837	d_h(:,:) = DTIME/DELTAT/REAL(NW) * d_h(:,:) &
838	+ (1.-DTIME/DELTAT) * dh_eddymix_gwd(:,:)
839	do ii=1,ngrid
840	d_t(ii,:) = d_h(ii,:) * (PP(ii,:) / PH(ii,1))**rcp
841	enddo
842
843	dh_eddymix_gwd(:,:)=d_h(:,:)
844
845	DO QQ=1,NQ
846	d_pq(:, :, QQ) =DTIME/DELTAT/REAL(NW) * d_pq(:, :, QQ) &
847	+ (1.-DTIME/DELTAT) * dq_eddymix_gwd(:, :, QQ)
848	endDO
849
850	do QQ=1,NQ
851	dq_eddymix_gwd(:, :, QQ)=d_pq(:, :, QQ)
852	ENDdo
853
854	END SUBROUTINE NONORO_GWD_MIX
855
856
857
858	! ========================================================
859	! Subroutines used to allocate/deallocate module variables
860	! ========================================================
861	SUBROUTINE ini_nonoro_gwd_mix(ngrid,nlayer,nq)
862
863	IMPLICIT NONE
864
865	INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns
866	INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
867	INTEGER, INTENT (in) :: nq ! number of atmospheric tracers
868
869	allocate(du_eddymix_gwd(ngrid,nlayer))
870	allocate(dv_eddymix_gwd(ngrid,nlayer))
871	allocate(dh_eddymix_gwd(ngrid,nlayer))
872	allocate(dq_eddymix_gwd(ngrid,nlayer,nq))
873	allocate(de_eddymix_rto(ngrid,nlayer+1))
874	allocate(df_eddymix_flx(ngrid,nlayer+1))
875
876	!du_eddymix_gwd(:,:)=0
877	!dv_eddymix_gwd(:,:)=0
878	! allocate(east_gwstress(ngrid,nlayer))
879	! east_gwstress(:,:)=0
880	! allocate(west_gwstress(ngrid,nlayer))
881	! west_gwstress(:,:)=0
882
883	END SUBROUTINE ini_nonoro_gwd_mix
884	! ----------------------------------
885	SUBROUTINE end_nonoro_gwd_mix
886
887	IMPLICIT NONE
888
889	if (allocated(du_eddymix_gwd)) deallocate(du_eddymix_gwd)
890	if (allocated(dv_eddymix_gwd)) deallocate(dv_eddymix_gwd)
891	if (allocated(dh_eddymix_gwd)) deallocate(dh_eddymix_gwd)
892	if (allocated(dq_eddymix_gwd)) deallocate(dq_eddymix_gwd)
893	if (allocated(de_eddymix_rto)) deallocate(de_eddymix_rto)
894	if (allocated(df_eddymix_flx)) deallocate(df_eddymix_flx)
895	! if (allocated(east_gwstress)) deallocate(east_gwstress)
896	! if (allocated(west_gwstress)) deallocate(west_gwstress)
897
898	END SUBROUTINE end_nonoro_gwd_mix
899	!-----------------------------------
900	! FUNCTION MAXLOCATION(gwd_normal,gwd_satura)
901
902	! implicit none
903
904	! INTEGER MAXLOCATION
905	! REAL gwd_normal,gwd_satura
906
907	! IF (gwd_normal .GT. gwd_satura ) THEN
908	! MAXLOCATION=1
909	! ELSEIF (gwd_normal .LT.gwd_satura) THEN
910	! MAXLOCATION=2
911	! ELSE
912	! MAXLOCATION=1
913	! ENDIF
914
915	! return
916
917	! END FUNCTION
918
919
920	END MODULE nonoro_gwd_mix_mod

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

Download in other formats:

Original Format