Context Navigation

nonoro_gwd_ran_mod.F90 @ 2220

Last change on this file since 2220 was 2220, checked in by dbardet, 5 years ago

Update the nonorographic gravity waves drag parametrization (from the r3599 of Earth s model LMDZ6): 1) add east_ and west_gwdstress variables, 2) delete aleas function: random waves are producted by the MOD function, 3) reproductibility concerning the number of procs is now validated, 4) changing name of some variables like RUW, RVW, ZOP, ZOM,..., 5) tendency of winds due to GW drag are now module variables and written in restart files.

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

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

Download in other formats:

Original Format