Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lmdz_gwd_precip.f90

Last change on this file was 6064, checked in by evignon, 6 weeks ago
suppression de print restants
File size: 15.6 KB

Line
1	!
2	! $Id: flott_gwd_rando_m.f90 5767 2025-07-10 10:04:51Z rkazeroni $
3	!
4	!$gpum horizontal klon
5	module lmdz_gwd_precip
6
7	IMPLICIT NONE
8
9	contains
10
11	!=========================================================================================
12	SUBROUTINE FLOTT_GWD_rando(klon, klev, DTIME, PP, presnivs, tt, uu, vv, prec, zustr, zvstr, d_u, &
13	d_v, east_gwstress, west_gwstress)
14
15	! Parametrization of the momentum flux deposition due to a discrete
16	! number of gravity waves.
17	! Author: F. Lott
18	! July, 12th, 2012
19	! Gaussian distribution of the source, source is precipitation
20	! Reference: Lott (JGR, vol 118, page 8897, 2013)
21
22	!ONLINE:
23	USE lmdz_gwd_ini, ONLY: GWD_RANDO_RUWMAX, GWD_RANDO_SAT, RPI, RG, RD, RCPD
24	USE lmdz_gwd_ini, ONLY: RLVTT, NK, NP, NO, NA, NW
25	USE lmdz_gwd_ini, ONLY: gwd_reproductibilite_mpiomp
26	USE assert_m, ONLY: assert
27
28	! OFFLINE:
29	! include "dimensions_mod.f90"
30	! include "dimphy.h"
31	! END OF DIFFERENCE ONLINE-OFFLINE
32
33	IMPLICIT NONE
34
35	CHARACTER(LEN=20), PARAMETER :: modname = 'flott_gwd_rando'
36	CHARACTER(LEN=80) :: abort_message
37
38	! 0. DECLARATIONS:
39
40	! 0.1 INPUTS
41	INTEGER, intent(in):: klon, klev ! horizontal and vertical dimensions
42	REAL, intent(in):: DTIME ! Time step of the Physics [s]
43	REAL, intent(in):: pp(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels [Pa]
44	REAL, intent(in):: presnivs(KLEV) ! equivalent pressure of model levels [Pa]
45	REAL, intent(in):: prec(KLON) ! (klon) Precipitation [kg/m^2/s]
46	REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels [K]
47	REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels [m/s]
48	REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels [m/s]
49
50	! 0.2 OUTPUTS
51	REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses [Pa]
52
53	REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV) ! wind increments [m/s]
54	REAL, intent(inout):: east_gwstress(KLON, KLEV) ! Profile of eastward stress [Pa]
55	REAL, intent(inout):: west_gwstress(KLON, KLEV) ! Profile of westward stress [Pa]
56
57	! O.3 INTERNAL ARRAYS
58	REAL BVLOW(klon)
59	REAL DZ ! Characteristic depth of the Source
60
61	INTEGER II, JJ, LL
62
63	! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
64
65	REAL DELTAT
66
67	! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
68
69	INTEGER JK, JP, JO, JW
70	REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
71	REAL CMAX ! standard deviation of the phase speed distribution
72	REAL RUWMAX, SAT ! ONLINE SPECIFIED IN run.def
73	REAL CPHA ! absolute PHASE VELOCITY frequency
74	REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude
75	REAL ZP(KLON, NW) ! Horizontal wavenumber angle
76	REAL ZO(KLON, NW) ! Absolute frequency !
77
78	! Waves Intr. freq. at the 1/2 lev surrounding the full level
79	REAL ZOM(KLON, NW), ZOP(KLON, NW)
80
81	! Wave EP-fluxes at the 2 semi levels surrounding the full level
82	REAL WWM(KLON, NW), WWP(KLON, NW)
83
84	REAL RUW0(KLON, NW) ! Fluxes at launching level
85
86	REAL RUWP(KLON, NW), RVWP(KLON, NW)
87	! Fluxes X and Y for each waves at 1/2 Levels
88
89	INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
90
91	REAL XLAUNCH ! Controle the launching altitude
92	REAL XTROP ! SORT of Tropopause altitude
93	REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
94	REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
95
96	REAL PRMAX ! Maximum value of PREC, and for which our linear formula
97	! for GWs parameterisation apply
98
99	! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
100
101	REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
102
103	! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
104
105	REAL H0 ! Characteristic Height of the atmosphere
106	REAL PR, TR ! Reference Pressure and Temperature
107
108	REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
109
110	REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
111	REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
112	REAL PSEC ! Security to avoid division by 0 pressure
113	REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
114	REAL BVSEC ! Security to avoid negative BVF
115	REAL RAN_NUM_1, RAN_NUM_2, RAN_NUM_3
116
117	REAL, DIMENSION(klev + 1) ::HREF
118
119	!-----------------------------------------------------------------
120
121	! 1. INITIALISATIONS
122
123	! 1.1 Basic parameter
124
125	! Are provided from elsewhere (latent heat of vaporization, dry
126	! gaz constant for air, gravity constant, heat capacity of dry air
127	! at constant pressure, earth rotation rate, pi).
128
129	! 1.2 Tuning parameters of V14
130
131	RDISS = 0.5 ! Diffusion parameter
132	! ONLINE
133	RUWMAX = GWD_RANDO_RUWMAX
134	SAT = gwd_rando_sat
135	!END ONLINE
136	! OFFLINE
137	! RUWMAX= 1.75 ! Launched flux
138	! SAT=0.25 ! Saturation parameter
139	! END OFFLINE
140
141	PRMAX = 20./24./3600.
142	! maximum of rain for which our theory applies (in kg/m^2/s)
143
144	! Characteristic depth of the source
145	DZ = 1000.
146	XLAUNCH = 0.5 ! Parameter that control launching altitude
147	XTROP = 0.2 ! Parameter that control tropopause altitude
148	DELTAT = 24.*3600. ! Time scale of the waves (first introduced in 9b)
149	! OFFLINE
150	! DELTAT=DTIME
151	! END OFFLINE
152
153	KMIN = 2.E-5
154	! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
155
156	KMAX = 1.E-3 ! Max horizontal wavenumber
157	CMAX = 30. ! Max phase speed velocity
158
159	TR = 240. ! Reference Temperature
160	PR = 101300. ! Reference pressure
161	H0 = RD*TR/RG ! Characteristic vertical scale height
162
163	BVSEC = 5.E-3 ! Security to avoid negative BVF
164	PSEC = 1.E-6 ! Security to avoid division by 0 pressure
165	ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
166
167	IF (1 == 0) THEN
168	!ONLINE
169	call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
170	size(vv, 1), size(zustr), size(zvstr), size(d_u, 1), &
171	size(d_v, 1), &
172	size(east_gwstress, 1), size(west_gwstress, 1)/), &
173	"FLOTT_GWD_RANDO klon")
174	call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
175	size(vv, 2), size(d_u, 2), size(d_v, 2), &
176	size(east_gwstress, 2), size(west_gwstress, 2)/), &
177	"FLOTT_GWD_RANDO klev")
178	!END ONLINE
179	END IF
180
181	IF (DELTAT < DTIME) THEN
182	abort_message = 'flott_gwd_rando: deltat < dtime!'
183	CALL abort_physic(modname, abort_message, 1)
184	END IF
185
186	IF (KLEV < NW) THEN
187	abort_message = 'flott_gwd_rando: you will have problem with random numbers'
188	CALL abort_physic(modname, abort_message, 1)
189	END IF
190
191	! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
192
193	! Pressure and Inv of pressure
194	DO LL = 2, KLEV
195	PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1)))/2.)
196	end DO
197	PH(:, KLEV + 1) = 0.
198	PH(:, 1) = 2.*PP(:, 1) - PH(:, 2)
199
200	! Launching altitude
201
202	! To get back to non-reproducible version by changing the number of cores
203	IF (gwd_reproductibilite_mpiomp) THEN
204	! rework the formula that computes PH as a function of PP=play
205	DO LL = 2, KLEV
206	HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1)))/2.)
207	end DO
208	HREF(KLEV + 1) = 0.
209	HREF(1) = 2.*presnivs(1) - HREF(2)
210	ELSE
211	HREF(1:KLEV) = PH(KLON/2, 1:KLEV)
212	END IF
213
214	LAUNCH = 0
215	LTROP = 0
216	DO LL = 1, KLEV
217	IF (HREF(LL)/HREF(1) > XLAUNCH) LAUNCH = LL
218	END DO
219	DO LL = 1, KLEV
220	IF (HREF(LL)/HREF(1) > XTROP) LTROP = LL
221	END DO
222	!LAUNCH=22 ; LTROP=33
223
224	! Log pressure vert. coordinate
225	DO LL = 1, KLEV + 1
226	ZH(:, LL) = H0*LOG(PR/(PH(:, LL) + PSEC))
227	end DO
228
229	! BV frequency
230	DO LL = 2, KLEV
231	! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
232	UH(:, LL) = 0.5*(TT(:, LL) + TT(:, LL - 1)) &
233	RD2/RCPD/H0*2 + (TT(:, LL) &
234	- TT(:, LL - 1))/(ZH(:, LL) - ZH(:, LL - 1))*RD/H0
235	end DO
236	BVLOW(:) = 0.5*(TT(:, LTROP) + TT(:, LAUNCH)) &
237	RD2/RCPD/H0*2 + (TT(:, LTROP) &
238	- TT(:, LAUNCH))/(ZH(:, LTROP) - ZH(:, LAUNCH))*RD/H0
239
240	UH(:, 1) = UH(:, 2)
241	UH(:, KLEV + 1) = UH(:, KLEV)
242	BV(:, 1) = UH(:, 2)
243	BV(:, KLEV + 1) = UH(:, KLEV)
244	! SMOOTHING THE BV HELPS
245	DO LL = 2, KLEV
246	BV(:, LL) = (UH(:, LL + 1) + 2.*UH(:, LL) + UH(:, LL - 1))/4.
247	end DO
248
249	BV = MAX(SQRT(MAX(BV, 0.)), BVSEC)
250	BVLOW = MAX(SQRT(MAX(BVLOW, 0.)), BVSEC)
251
252	! WINDS
253	DO LL = 2, KLEV
254	UH(:, LL) = 0.5*(UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
255	VH(:, LL) = 0.5*(VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
256	end DO
257	UH(:, 1) = 0.
258	VH(:, 1) = 0.
259	UH(:, KLEV + 1) = UU(:, KLEV)
260	VH(:, KLEV + 1) = VV(:, KLEV)
261
262	! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE
263
264	! The mod functions of weird arguments are used to produce the
265	! waves characteristics in an almost stochastic way
266
267	DO JW = 1, NW
268	! Angle
269	DO II = 1, KLON
270	! Angle (0 or PI so far)
271	RAN_NUM_1 = MOD(TT(II, JW)*10., 1.)
272	RAN_NUM_2 = MOD(TT(II, JW)*100., 1.)
273	ZP(II, JW) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
274	*RPI/2.
275	! Horizontal wavenumber amplitude
276	ZK(II, JW) = KMIN + (KMAX - KMIN)*RAN_NUM_2
277	! Horizontal phase speed
278	CPHA = 0.
279	DO JJ = 1, NA
280	RAN_NUM_3 = MOD(TT(II, JW + 3JJ)*2, 1.)
281	CPHA = CPHA + &
282	CMAX2.(RAN_NUM_3 - 0.5)SQRT(3.)/SQRT(NA1.)
283	END DO
284	IF (CPHA .LT. 0.) THEN
285	CPHA = -1.*CPHA
286	ZP(II, JW) = ZP(II, JW) + RPI
287	END IF
288	! Absolute frequency is imposed
289	ZO(II, JW) = CPHA*ZK(II, JW)
290	! Intrinsic frequency is imposed
291	ZO(II, JW) = ZO(II, JW) &
292	+ ZK(II, JW)COS(ZP(II, JW))UH(II, LAUNCH) &
293	+ ZK(II, JW)SIN(ZP(II, JW))VH(II, LAUNCH)
294	! Momentum flux at launch lev
295	RUW0(II, JW) = RUWMAX
296	END DO
297	END DO
298
299	! 4. COMPUTE THE FLUXES
300
301	! 4.1 Vertical velocity at launching altitude to ensure
302	! the correct value to the imposed fluxes.
303
304	DO JW = 1, NW
305
306	! Evaluate intrinsic frequency at launching altitude:
307	ZOP(:, JW) = ZO(:, JW) &
308	- ZK(:, JW)COS(ZP(:, JW))UH(:, LAUNCH) &
309	- ZK(:, JW)SIN(ZP(:, JW))VH(:, LAUNCH)
310
311	! VERSION WITH CONVECTIVE SOURCE
312
313	! Vertical velocity at launch level, value to ensure the
314	! imposed factor related to the convective forcing:
315	! precipitations.
316
317	! tanh limitation to values above prmax:
318	WWP(:, JW) = RUW0(:, JW) &
319	(RD/RCPD/H0RLVTTPRMAXTANH(PREC(:)/PRMAX))**2
320
321	! Factor related to the characteristics of the waves:
322	WWP(:, JW) = WWP(:, JW)ZK(:, JW)*3/KMIN/BVLOW(:) &
323	/MAX(ABS(ZOP(:, JW)), ZOISEC)**3
324
325	! Moderation by the depth of the source (dz here):
326	WWP(:, JW) = WWP(:, JW) &
327	EXP(-BVLOW(:)2/MAX(ABS(ZOP(:, JW)), ZOISEC)2ZK(:, JW)**2 &
328	DZ*2)
329
330	! Put the stress in the right direction:
331	RUWP(:, JW) = ZOP(:, JW)/MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
332	BV(:, LAUNCH)COS(ZP(:, JW))WWP(:, JW)*2
333	RVWP(:, JW) = ZOP(:, JW)/MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
334	BV(:, LAUNCH)SIN(ZP(:, JW))WWP(:, JW)*2
335	end DO
336
337	! 4.2 Uniform values below the launching altitude
338
339	DO LL = 1, LAUNCH
340	RUW(:, LL) = 0
341	RVW(:, LL) = 0
342	DO JW = 1, NW
343	RUW(:, LL) = RUW(:, LL) + RUWP(:, JW)
344	RVW(:, LL) = RVW(:, LL) + RVWP(:, JW)
345	end DO
346	end DO
347
348	! 4.3 Loop over altitudes, with passage from one level to the next
349	! done by i) conserving the EP flux, ii) dissipating a little,
350	! iii) testing critical levels, and vi) testing the breaking.
351
352	DO LL = LAUNCH, KLEV - 1
353	! Warning: all the physics is here (passage from one level
354	! to the next)
355	DO JW = 1, NW
356	ZOM(:, JW) = ZOP(:, JW)
357	WWM(:, JW) = WWP(:, JW)
358	! Intrinsic Frequency
359	ZOP(:, JW) = ZO(:, JW) - ZK(:, JW)COS(ZP(:, JW))UH(:, LL + 1) &
360	- ZK(:, JW)SIN(ZP(:, JW))VH(:, LL + 1)
361
362	! No breaking (Eq.6)
363	! Dissipation (Eq. 8)
364	WWP(:, JW) = WWM(:, JW)EXP(-4.RDISS*PR/(PH(:, LL + 1) &
365	+ PH(:, LL))((BV(:, LL + 1) + BV(:, LL))/2.)*3 &
366	/MAX(ABS(ZOP(:, JW) + ZOM(:, JW))/2., ZOISEC)**4 &
367	ZK(:, JW)3(ZH(:, LL + 1) - ZH(:, LL)))
368
369	! Critical levels (forced to zero if intrinsic frequency changes sign)
370	! Saturation (Eq. 12)
371	WWP(:, JW) = min(WWP(:, JW), MAX(0., &
372	SIGN(1., ZOP(:, JW)ZOM(:, JW)))ABS(ZOP(:, JW))**3 &
373	/BV(:, LL + 1)EXP(-ZH(:, LL + 1)/H0)KMIN**2 &
374	SAT2/ZK(:, JW)*4)
375	end DO
376
377	! Evaluate EP-flux from Eq. 7 and give the right orientation to
378	! the stress
379
380	DO JW = 1, NW
381	RUWP(:, JW) = SIGN(1., ZOP(:, JW))COS(ZP(:, JW))WWP(:, JW)
382	RVWP(:, JW) = SIGN(1., ZOP(:, JW))SIN(ZP(:, JW))WWP(:, JW)
383	end DO
384
385	RUW(:, LL + 1) = 0.
386	RVW(:, LL + 1) = 0.
387
388	DO JW = 1, NW
389	RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:, JW)
390	RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:, JW)
391	EAST_GWSTRESS(:, LL) = EAST_GWSTRESS(:, LL) + MAX(0., RUWP(:, JW))/REAL(NW)
392	WEST_GWSTRESS(:, LL) = WEST_GWSTRESS(:, LL) + MIN(0., RUWP(:, JW))/REAL(NW)
393	end DO
394	end DO
395	! OFFLINE ONLY
396	! PRINT *,'SAT PROFILE:'
397	! DO LL=1,KLEV
398	! PRINT ,ZH(KLON/2,LL)/1000.,SAT(2.+TANH(ZH(KLON/2,LL)/H0-8.))
399	! ENDDO
400
401	! 5 CALCUL DES TENDANCES:
402
403	! 5.1 Rectification des flux au sommet et dans les basses couches
404
405	RUW(:, KLEV + 1) = 0.
406	RVW(:, KLEV + 1) = 0.
407	RUW(:, 1) = RUW(:, LAUNCH)
408	RVW(:, 1) = RVW(:, LAUNCH)
409	DO LL = 1, LAUNCH
410	RUW(:, LL) = RUW(:, LAUNCH + 1)
411	RVW(:, LL) = RVW(:, LAUNCH + 1)
412	EAST_GWSTRESS(:, LL) = EAST_GWSTRESS(:, LAUNCH)
413	WEST_GWSTRESS(:, LL) = WEST_GWSTRESS(:, LAUNCH)
414	end DO
415
416	! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
417	DO LL = 1, KLEV
418	D_U(:, LL) = (1.-DTIME/DELTAT)D_U(:, LL) + DTIME/DELTAT/REAL(NW) &
419	RG*(RUW(:, LL + 1) - RUW(:, LL)) &
420	/(PH(:, LL + 1) - PH(:, LL))*DTIME
421	! NO AR-1 FOR MERIDIONAL TENDENCIES
422	D_V(:, LL) = 1./REAL(NW)* &
423	RG*(RVW(:, LL + 1) - RVW(:, LL)) &
424	/(PH(:, LL + 1) - PH(:, LL))*DTIME
425	END DO
426
427	! Cosmetic: evaluation of the cumulated stress
428	ZUSTR = 0.
429	ZVSTR = 0.
430	DO LL = 1, KLEV
431	ZUSTR = ZUSTR + D_U(:, LL)/RG*(PH(:, LL + 1) - PH(:, LL))/DTIME
432	ZVSTR = ZVSTR + D_V(:, LL)/RG*(PH(:, LL + 1) - PH(:, LL))/DTIME
433	END DO
434
435	END SUBROUTINE FLOTT_GWD_RANDO
436
437	end module lmdz_gwd_precip

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

Download in other formats:

Original Format