Context Navigation

nonoro_gwd_ran_mod.F90 @ 2640

Last change on this file since 2640 was 2640, checked in by emillour, 3 years ago
Mars GCM: Update for non-orographic Gravity Waves; some key parameters can now be set via callphys.def JL
File size: 25.6 KB

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

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

Download in other formats:

Original Format