Context Navigation

nonoro_gwd_mix_mod.F90 @ 3393

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

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

Download in other formats:

Original Format