Context Navigation

nonoro_gwd_mix_mod.F90 @ 3225

Last change on this file since 3225 was 3144, checked in by jliu, 15 months ago

Mars PCM
New scheme "nonoro_gwd_mix_mod" to simulate non-orographic gravity waves induced mixing. The scheme flag "calljliu_gwimix" has been put to "false" as default.
There are three tunable parameters that have been set up, as "nonoro_gwimixing_eff=0.1", "nonoro_gwimixing_eff1=0.1", "nonoro_gwimixing_vdl=1.5".
Remember to use full-layers model configuration (64x48x73) when you try to call this scheme because we need to evaluate the saturation altitude of each monochromatic wave.
More than ~75% of the waves saturated at altitude above 120 km.
Jliu

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

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

Download in other formats:

Original Format