Context Navigation

nonoro_gwd_ran_mod.F90 @ 2257

Last change on this file since 2257 was 2225, checked in by emillour, 6 years ago
Mars GCM: Further cleanup of nonorographic gravity waves drag parametrization east_gwstress and west_gwstress need be initialized and are better as module variables. EM
File size: 21.3 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	CONTAINS
11
12	SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, pp, &
13	zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
14	zustr,zvstr,d_t, d_u, d_v)
15
16	!--------------------------------------------------------------------------------
17	! Parametrization of the momentum flux deposition due to a discrete
18	! number of gravity waves.
19	! F. Lott
20	! Version 14, Gaussian distribution of the source
21	! LMDz model online version
22	! ADAPTED FOR VENUS / F. LOTT + S. LEBONNOIS
23	! Version adapted on 03/04/2013:
24	! - input flux compensated in the deepest layers
25	!
26	! ADAPTED FOR MARS G.GILLI 02/2016
27	! Revision with F.Forget 06/2016 Variable EP-flux according to
28	! PBL variation (max velocity thermals)
29	! UPDATED D.BARDET 01/2020 - reproductibility of the
30	! launching altitude calculation
31	! - wave characteristic
32	! calculation using MOD
33	! - adding east_gwstress and
34	! west_gwstress variables
35	!---------------------------------------------------------------------------------
36
37	use comcstfi_h, only: g, pi, cpp, r
38	USE ioipsl_getin_p_mod, ONLY : getin_p
39	use assert_m, only : assert
40	use vertical_layers_mod, only : presnivs
41	implicit none
42
43	include "dimensions.h"
44	include "paramet.h"
45	include "yoegwd.h"
46	include "callkeys.h"
47
48	CHARACTER (LEN=20) :: modname='flott_gwd_rando'
49	CHARACTER (LEN=80) :: abort_message
50
51
52	! 0. DECLARATIONS:
53
54	! 0.1 INPUTS
55	INTEGER, intent(in):: ngrid, nlayer
56	REAL, intent(in):: DTIME ! Time step of the Physics
57	REAL, intent(in):: zmax_therm(ngrid) ! altitude of max velocity thermals (m)
58	REAL, intent(in):: pp(ngrid,nlayer) ! Pressure at full levels
59	REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels
60	REAL, intent(in):: pu(ngrid,nlayer),pv(ngrid,nlayer) ! Hor winds at full levels
61	REAL,INTENT(in) :: pdt(ngrid,nlayer) ! tendency on temperature
62	REAL,INTENT(in) :: pdu(ngrid,nlayer) ! tendency on zonal wind
63	REAL,INTENT(in) :: pdv(ngrid,nlayer) ! tendency on meridional wind
64
65	! 0.2 OUTPUTS
66	REAL, intent(out):: zustr(ngrid), zvstr(ngrid) ! Surface Stresses
67	REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on Temp.
68	REAL, intent(out):: d_u(ngrid, nlayer), d_v(ngrid, nlayer) ! tendency on winds
69
70	! O.3 INTERNAL ARRAYS
71	REAL :: TT(ngrid, nlayer) ! Temp at full levels
72	REAL :: UU(ngrid, nlayer) , VV(ngrid, nlayer) ! Hor winds at full levels
73	REAL :: BVLOW(ngrid)
74	REAL :: DZ
75
76	INTEGER II, JJ, LL
77
78	! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
79
80	REAL, parameter:: DELTAT = 24. * 3600.
81
82	! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
83
84	INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers
85	INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed
86	INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed
87	INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves
88	INTEGER JK, JP, JO, JW
89	INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed
90	REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber
91	REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber
92	REAL, parameter:: cmax = 30. ! Max horizontal absolute phase velocity
93	REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity
94	REAL CPHA ! absolute PHASE VELOCITY frequency
95	REAL ZK(NW, ngrid) ! Horizontal wavenumber amplitude
96	REAL ZP(NW, ngrid) ! Horizontal wavenumber angle
97	REAL ZO(NW, ngrid) ! Absolute frequency
98
99	REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
100	REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
101	REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level
102	REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level
103	REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
104	REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
105	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)
106	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)
107	REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0)
108	REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX)
109	INTEGER LAUNCH ! Launching altitude
110	REAL, parameter:: xlaunch = 0.4 ! Control the launching altitude
111	REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.)
112
113
114	REAL PREC(ngrid)
115	REAL PRMAX ! Maximum value of PREC, and for which our linear formula
116
117
118	! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
119	REAL, parameter:: sat = 1. ! saturation parameter
120	REAL, parameter:: rdiss = 1. ! coefficient of dissipation
121	REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency
122
123	! 0.3.3 Background flow at 1/2 levels and vertical coordinate
124	REAL H0bis(ngrid, nlayer) ! Characteristic Height of the atmosphere
125	REAL, save:: H0 ! Characteristic Height of the atmosphere
126	REAL, parameter:: pr = 250 ! Reference pressure [Pa]
127	REAL, parameter:: tr = 220. ! Reference temperature [K]
128	REAL ZH(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0)
129	REAL ZHbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H)
130	REAL UH(ngrid, nlayer + 1) ! zonal wind at 1/2 levels
131	REAL VH(ngrid, nlayer + 1) ! meridional wind at 1/2 levels
132	REAL PH(ngrid, nlayer + 1) ! Pressure at 1/2 levels
133	REAL, parameter:: psec = 1.e-6 ! Security to avoid division by 0 pressure
134	REAL BV(ngrid, nlayer + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
135	REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV
136	REAL HREF(nlayer + 1) ! Reference altitude for launching alt.
137
138
139	! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
140	logical,save :: output=.false.
141	! CAUTION ! IF output is .true. THEN change NO to 10 at least !
142	character*14 outform
143	character*2 str2
144	integer ieq
145
146	REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
147
148
149	LOGICAL,SAVE :: firstcall = .true.
150
151
152	!-----------------------------------------------------------------
153	! 1. INITIALISATIONS
154
155	IF (firstcall) THEN
156	write(,) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
157	epflux_max = 7.E-7 ! Mars' value !!
158	call getin_p("nonoro_gwd_epflux_max", epflux_max)
159	write(,) "nonoro_gwd_ran: epflux_max=", epflux_max
160	! Characteristic vertical scale height
161	H0 = r * tr / g
162	! Control
163	if (deltat .LT. dtime) THEN
164	call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
165	endif
166	if (nlayer .LT. nw) THEN
167	call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
168	endif
169	firstcall = .false.
170	ENDIF
171
172	gwd_convective_source=.false.
173
174	! Compute current values of temperature and winds
175	tt(:,:)=pt(:,:)+dtime*pdt(:,:)
176	uu(:,:)=pu(:,:)+dtime*pdu(:,:)
177	vv(:,:)=pv(:,:)+dtime*pdv(:,:)
178
179
180
181	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
182	!-------------------------------------------------------------
183
184	!Online output
185	if (output) OPEN(11,file="impact-gwno.dat")
186
187	! Pressure and Inv of pressure, Temperature / at 1/2 level
188	DO LL = 2, nlayer
189	PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
190	end DO
191
192	PH(:, nlayer + 1) = 0.
193	PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
194
195	! Launching altitude
196
197	!Pour revenir a la version non reproductible en changeant le nombre de
198	!process
199	! Reprend la formule qui calcule PH en fonction de PP=play
200	DO LL = 2, nlayer
201	HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
202	end DO
203	HREF(nlayer + 1) = 0.
204	HREF(1) = 2. * presnivs(1) - HREF(2)
205
206	LAUNCH=0
207	DO LL = 1, nlayer
208	IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
209	ENDDO
210
211	if (output) print*, " WE ARE IN FLOTT GW SCHEME "
212
213	! Log pressure vert. coordinate
214	DO LL = 1, nlayer + 1
215	ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
216	end DO
217
218	if (output) then
219	! altitude above surface
220	ZHbis(:,1) = 0.
221	DO LL = 2, nlayer + 1
222	H0bis(:, LL-1) = r * TT(:, LL-1) / g
223	ZHbis(:, LL) = ZHbis(:, LL-1) &
224	+ H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
225	end DO
226	endif
227
228	! Winds and BV frequency
229	DO LL = 2, nlayer
230	UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
231	VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
232	! GG test
233	!print*, 'TT, UH, VH, ZH at launch', TT(ngrid/2,LAUNCH), UH(ngrid/2,LAUNCH),VH(ngrid/2, LAUNCH), ZH(ngrid/2,LAUNCH)
234	! BVSEC: BV Frequency
235	BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
236	* r2 / cpp / H02 + (TT(:, LL) &
237	- TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * r / H0
238	BV(:,LL) =SQRT(MAX(BVSEC,BV(:,LL)))
239	end DO
240	!GG test
241	!print*, 'BV freq in flott_gwnoro:',LAUNCH, BV(ngrid/2, LAUNCH)
242
243	BV(:, 1) = BV(:, 2)
244	UH(:, 1) = 0.
245	VH(:, 1) = 0.
246	BV(:, nlayer + 1) = BV(:, nlayer)
247	UH(:, nlayer + 1) = UU(:, nlayer)
248	VH(:, nlayer + 1) = VV(:, nlayer)
249
250
251	! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
252	!-------------------------------------------
253
254	! The mod function of here a weird arguments
255	! are used to produce the waves characteristics
256	! in a stochastic way
257
258	DO JW = 1, NW
259	! Angle
260	DO II = 1, ngrid
261	! Angle (0 or PI so far)
262	RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
263	RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
264	ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
265	* PI / 2.
266	! Horizontal wavenumber amplitude
267	ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
268	! Horizontal phase speed
269	CPHA = 0.
270	DO JJ = 1, NA
271	RAN_NUM_3=MOD(TT(II, JW+3JJ)*2, 1.)
272	CPHA = CPHA + &
273	CMAX2.(RAN_NUM_3 -0.5)SQRT(3.)/SQRT(NA1.)
274	END DO
275	IF (CPHA.LT.0.) THEN
276	CPHA = -1.*CPHA
277	ZP(JW,II) = ZP(JW,II) + PI
278	ENDIF
279	!Online output: linear
280	if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1)
281	! Intrinsic frequency
282	ZO(JW, II) = CPHA * ZK(JW, II)
283	! Intrinsic frequency is imposed
284	ZO(JW, II) = ZO(JW, II) &
285	+ ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
286	+ ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
287	! Momentum flux at launch lev
288	! epflux_0(JW, II) = epflux_max / REAL(NW) &
289	epflux_0(JW, II) = epflux_max &
290	* MOD(100. * (UU(II, JW)2 + VV(II, JW)2), 1.)
291	!Online output: fixed to max
292	if (output) epflux_0(JW, II) = epflux_max
293	ENDDO
294	end DO
295
296	! 4. COMPUTE THE FLUXES
297	!--------------------------
298
299	! 4.1 Vertical velocity at launching altitude to ensure
300	! the correct value to the imposed fluxes.
301	!
302	DO JW = 1, NW
303	! Evaluate intrinsic frequency at launching altitude:
304	intr_freq_p(JW, :) = ZO(JW, :) &
305	- ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
306	- ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
307	end DO
308
309	IF (gwd_convective_source) THEN
310	DO JW = 1, NW
311	! VERSION WITH CONVECTIVE SOURCE (designed for Earth)
312
313	! Vertical velocity at launch level, value to ensure the
314	! imposed mmt flux factor related to the convective forcing:
315	! precipitations.
316
317	! tanh limitation to values above prmax:
318	! WWP(JW, :) = epflux_0(JW, :) &
319	! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
320	! Here, we neglected the kinetic energy providing of the thermodynamic
321	! phase change
322
323	!
324
325	! Factor related to the characteristics of the waves:
326	WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:) &
327	/ MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**3
328
329	! Moderation by the depth of the source (dz here):
330	WWP(JW, :) = WWP(JW, :) &
331	* EXP(- BVLOW(:)2 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)2 &
332	* ZK(JW, :)*2 DZ**2)
333
334	! Put the stress in the right direction:
335	u_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
336	* BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
337	v_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
338	* BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
339	end DO
340	ELSE ! VERSION WITHOUT CONVECTIVE SOURCE
341	! Vertical velocity at launch level, value to ensure the imposed
342	! mom flux:
343	DO JW = 1, NW
344	! WW is directly a flux, here, not vertical velocity anymore
345	WWP(JW, :) = epflux_0(JW,:)
346	u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
347	v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
348
349	end DO
350	ENDIF
351	! 4.2 Initial flux at launching altitude
352
353	u_epflux_tot(:, LAUNCH) = 0
354	v_epflux_tot(:, LAUNCH) = 0
355	DO JW = 1, NW
356	u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :)
357	v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :)
358	end DO
359
360	! 4.3 Loop over altitudes, with passage from one level to the
361	! next done by i) conserving the EP flux, ii) dissipating
362	! a little, iii) testing critical levels, and vi) testing
363	! the breaking.
364
365	!Online output
366	if (output) then
367	ieq=ngrid/2+1
368	write(str2,'(i2)') NW+2
369	outform="("//str2//"(E12.4,1X))"
370	WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
371	(ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW, IEQ)), JW = 1, NW)
372	endif
373
374	DO LL = LAUNCH, nlayer - 1
375
376
377	! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL
378	! TO THE NEXT)
379	DO JW = 1, NW
380	intr_freq_m(JW, :) = intr_freq_p(JW, :)
381	WWM(JW, :) = WWP(JW, :)
382	! Intrinsic Frequency
383	intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) &
384	- ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1)
385
386	WWP(JW, :) = MIN( &
387	! No breaking (Eq.6)
388	WWM(JW, :) &
389	! Dissipation (Eq. 8):
390	* EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
391	* ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
392	/ MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 &
393	* ZK(JW, :)*3 (ZH(:, LL + 1) - ZH(:, LL))), &
394	! Critical levels (forced to zero if intrinsic frequency changes sign)
395	MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) &
396	! Saturation (Eq. 12)
397	* ABS(intr_freq_p(JW, :))**3 /BV(:, LL+1) &
398	* EXP(-ZH(:, LL + 1)/H0) * SAT*2KMIN2/ZK(JW, :)4)
399	end DO
400
401	! END OF W(KB)ARNING
402	! Evaluate EP-flux from Eq. 7 and
403	! Give the right orientation to the stress
404	DO JW = 1, NW
405	u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) * WWP(JW, :)
406	v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) * WWP(JW, :)
407	end DO
408
409	u_epflux_tot(:, LL + 1) = 0.
410	v_epflux_tot(:, LL + 1) = 0.
411
412	DO JW = 1, NW
413	u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :)
414	v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :)
415	EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW)
416	WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW)
417	end DO
418	!Online output
419	if (output) then
420	do JW=1,NW
421	if(u_epflux_p(JW, IEQ).gt.0.) then
422	u_epflux_p(JW, IEQ) = max(u_epflux_p(JW, IEQ), 1.e-99)
423	else
424	u_epflux_p(JW, IEQ) = min(u_epflux_p(JW, IEQ), -1.e-99)
425	endif
426	enddo
427	WRITE(11,outform) ZH(IEQ, LL+1) / 1000., &
428	ZHbis(IEQ, LL+1) / 1000., &
429	(u_epflux_p(JW, IEQ), JW = 1, NW)
430	endif
431
432	end DO
433
434	! 5 CALCUL DES TENDANCES:
435	!------------------------
436
437	! 5.1 Rectification des flux au sommet et dans les basses couches:
438
439	! Attention, ici c'est le total sur toutes les ondes...
440
441	u_epflux_tot(:, nlayer + 1) = 0.
442	v_epflux_tot(:, nlayer + 1) = 0.
443
444	! Here, big change compared to FLott version:
445	! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
446	! over the layers max(1,LAUNCH-3) to LAUNCH-1
447	DO LL = 1, max(1,LAUNCH-3)
448	u_epflux_tot(:, LL) = 0.
449	v_epflux_tot(:, LL) = 0.
450	end DO
451	DO LL = max(2,LAUNCH-2), LAUNCH-1
452	u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * &
453	(PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
454	v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * &
455	(PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
456	EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + &
457	EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
458	(PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
459	WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + &
460	WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
461	(PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
462	end DO
463	! This way, the total flux from GW is zero, but there is a net transport
464	! (upward) that should be compensated by circulation
465	! and induce additional friction at the surface
466
467	!Online output
468	if (output) then
469	DO LL = 1, nlayer - 1
470	WRITE(11,*) ZHbis(IEQ, LL)/1000.,u_epflux_tot(IEQ,LL)
471	end DO
472	CLOSE(11)
473	stop
474	endif
475
476
477	! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
478	!---------------------------------------------
479	DO LL = 1, nlayer
480	d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
481	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
482	d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
483	/ (PH(:, LL + 1) - PH(:, LL)) * DTIME
484	ENDDO
485	d_t(:,:) = 0.
486
487	! 5.3 Update tendency of wind with the previous (and saved) values
488	!-----------------------------------------------------------------
489	d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) &
490	+ (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
491	d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) &
492	+ (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
493	du_nonoro_gwd(:,:) = d_u(:,:)
494	dv_nonoro_gwd(:,:) = d_v(:,:)
495
496	! Cosmetic: evaluation of the cumulated stress
497
498	ZUSTR(:) = 0.
499	ZVSTR(:) = 0.
500	DO LL = 1, nlayer
501	ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
502	ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
503	ENDDO
504
505	END SUBROUTINE NONORO_GWD_RAN
506
507	! ========================================================
508	! Subroutines used to allocate/deallocate module variables
509	! ========================================================
510
511	SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
512
513	IMPLICIT NONE
514
515	INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns
516	INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
517
518	allocate(du_nonoro_gwd(ngrid,nlayer))
519	allocate(dv_nonoro_gwd(ngrid,nlayer))
520	allocate(east_gwstress(ngrid,nlayer))
521	east_gwstress(:,:)=0
522	allocate(west_gwstress(ngrid,nlayer))
523	west_gwstress(:,:)=0
524
525	END SUBROUTINE ini_nonoro_gwd_ran
526	! ----------------------------------
527	SUBROUTINE end_nonoro_gwd_ran
528
529	IMPLICIT NONE
530
531	if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
532	if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
533	if (allocated(east_gwstress)) deallocate(east_gwstress)
534	if (allocated(west_gwstress)) deallocate(west_gwstress)
535
536	END SUBROUTINE end_nonoro_gwd_ran
537
538	END MODULE nonoro_gwd_ran_mod

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

Download in other formats:

Original Format