source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/flott_gwd_rando_m.F90 @ 5440

Last change on this file since 5440 was 1938, checked in by lguez, 11 years ago

New run-time parameters: ok_gwd_rando, ruwmax and sat. Logical
variable ok_gwd_rando activates the new stochastic parameterization of
the drag of non-orographic gravity waves, by F. Lott. Default value of
ok_gwd_rando is false. A run with ok_gwd_rando creates new variables
in restartphy.nc: du_gwd_rando and dv_gwd_rando. Possible output in
history files of new variables du_gwd_rando, dv_gwd_rando and
vstr_gwd_rando.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory.
File size: 12.3 KB
Line 
1module FLOTT_GWD_rando_m
2
3  implicit none
4
5contains
6
7  SUBROUTINE FLOTT_GWD_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, &
8       d_v)
9
10    ! Parametrization of the momentum flux deposition due to a discrete
11    ! number of gravity waves.
12    ! Author: F. Lott
13    ! July, 12th, 2012
14    ! Gaussian distribution of the source, source is precipitation
15    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
16
17    use dimphy, only: klon, klev
18    use assert_m, only: assert
19
20    include "YOMCST.h"
21    include "clesphys.h"
22    include "YOEGWD.h"
23
24    ! 0. DECLARATIONS:
25
26    ! 0.1 INPUTS
27    REAL, intent(in)::DTIME ! Time step of the Physics
28    REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels
29    REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s)
30    REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
31    REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
32    REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
33
34    ! 0.2 OUTPUTS
35    REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
36
37    REAL, intent(inout):: d_u(:, :), d_v(:, :)
38    ! (KLON, KLEV) tendencies on winds
39
40    ! O.3 INTERNAL ARRAYS
41    REAL BVLOW(klon)
42
43    INTEGER II, LL
44
45    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
46
47    REAL DELTAT
48
49    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
50
51    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
52    INTEGER JK, JP, JO, JW
53    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
54    REAL CMIN, CMAX ! Min and Max absolute ph. vel.
55    REAL CPHA ! absolute PHASE VELOCITY frequency
56    REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
57    REAL ZP(NW, KLON) ! Horizontal wavenumber angle
58    REAL ZO(NW, KLON) ! Absolute frequency !
59
60    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
61    REAL ZOM(NW, KLON), ZOP(NW, KLON)
62
63    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
64    REAL WWM(NW, KLON), WWP(NW, KLON)
65
66    REAL RUW0(NW, KLON) ! Fluxes at launching level
67
68    REAL RUWP(NW, KLON), RVWP(NW, KLON)
69    ! Fluxes X and Y for each waves at 1/2 Levels
70
71    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
72
73    REAL XLAUNCH ! Controle the launching altitude
74    REAL XTROP ! SORT of Tropopause altitude
75    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
76    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
77
78    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
79    ! for GWs parameterisation apply
80
81    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
82
83    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
84
85    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
86
87    REAL H0 ! Characteristic Height of the atmosphere
88    REAL DZ ! Characteristic depth of the source!
89    REAL PR, TR ! Reference Pressure and Temperature
90
91    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
92
93    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
94    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
95    REAL PSEC ! Security to avoid division by 0 pressure
96    REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
97    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
98    REAL BVSEC ! Security to avoid negative BVF
99
100    !-----------------------------------------------------------------
101
102    ! 1. INITIALISATIONS
103
104    ! 1.1 Basic parameter
105
106    ! Are provided from elsewhere (latent heat of vaporization, dry
107    ! gaz constant for air, gravity constant, heat capacity of dry air
108    ! at constant pressure, earth rotation rate, pi).
109
110    ! 1.2 Tuning parameters of V14
111
112    RDISS = 1. ! Diffusion parameter
113
114    PRMAX = 20. / 24. /3600.
115    ! maximum of rain for which our theory applies (in kg/m^2/s)
116
117    DZ = 1000. ! Characteristic depth of the source
118    XLAUNCH=0.5 ! Parameter that control launching altitude
119    XTROP=0.2 ! Parameter that control tropopause altitude
120    DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)
121
122    KMIN = 2.E-5
123    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
124
125    KMAX = 1.E-3 ! Max horizontal wavenumber
126    CMIN = 1. ! Min phase velocity
127    CMAX = 50. ! Max phase speed velocity
128
129    TR = 240. ! Reference Temperature
130    PR = 101300. ! Reference pressure
131    H0 = RD * TR / RG ! Characteristic vertical scale height
132
133    BVSEC = 5.E-3 ! Security to avoid negative BVF
134    PSEC = 1.E-6 ! Security to avoid division by 0 pressure
135    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
136
137    call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
138         size(vv, 1), size(prec), size(zustr), size(zvstr), size(d_u, 1), &
139         size(d_v, 1)/), "FLOTT_GWD_RANDO klon")
140    call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
141         size(vv, 2), size(d_u, 2), size(d_v, 2)/), "FLOTT_GWD_RANDO klev")
142
143    IF(DELTAT < DTIME)THEN
144       PRINT *, 'flott_gwd_rando: deltat < dtime!'
145       STOP 1
146    ENDIF
147
148    IF (KLEV < NW) THEN
149       PRINT *, 'flott_gwd_rando: you will have problem with random numbers'
150       STOP 1
151    ENDIF
152
153    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
154
155    ! Pressure and Inv of pressure
156    DO LL = 2, KLEV
157       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
158       PHM1(:, LL) = 1. / PH(:, LL)
159    end DO
160
161    PH(:, KLEV + 1) = 0.
162    PHM1(:, KLEV + 1) = 1. / PSEC
163    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
164
165    ! Launching altitude
166
167    LAUNCH=0
168    LTROP =0
169    DO LL = 1, KLEV
170       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XLAUNCH) LAUNCH = LL
171    ENDDO
172    DO LL = 1, KLEV
173       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XTROP) LTROP = LL
174    ENDDO
175
176    ! Log pressure vert. coordinate
177    DO LL = 1, KLEV + 1
178       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
179    end DO
180
181    ! BV frequency
182    DO LL = 2, KLEV
183       ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
184       UH(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
185            * RD**2 / RCPD / H0**2 + (TT(:, LL) &
186            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
187    end DO
188    BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
189         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
190         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
191
192    UH(:, 1) = UH(:, 2)
193    UH(:, KLEV + 1) = UH(:, KLEV)
194    BV(:, 1) = UH(:, 2)
195    BV(:, KLEV + 1) = UH(:, KLEV)
196    ! SMOOTHING THE BV HELPS
197    DO LL = 2, KLEV
198       BV(:, LL)=(UH(:, LL+1)+2.*UH(:, LL)+UH(:, LL-1))/4.
199    end DO
200
201    BV=MAX(SQRT(MAX(BV, 0.)), BVSEC)
202    BVLOW=MAX(SQRT(MAX(BVLOW, 0.)), BVSEC)
203
204    ! WINDS
205    DO LL = 2, KLEV
206       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
207       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
208    end DO
209    UH(:, 1) = 0.
210    VH(:, 1) = 0.
211    UH(:, KLEV + 1) = UU(:, KLEV)
212    VH(:, KLEV + 1) = VV(:, KLEV)
213
214    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE
215
216    ! The mod functions of weird arguments are used to produce the
217    ! waves characteristics in an almost stochastic way
218
219    JW = 0
220    DO JP = 1, NP
221       DO JK = 1, NK
222          DO JO = 1, NO
223             JW = JW + 1
224             ! Angle
225             DO II = 1, KLON
226                ! Angle (0 or PI so far)
227                ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
228                     * RPI / 2.
229                ! Horizontal wavenumber amplitude
230                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
231                ! Horizontal phase speed
232                CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
233                ! Absolute frequency is imposed
234                ZO(JW, II) = CPHA * ZK(JW, II)
235                ! Intrinsic frequency is imposed
236                ZO(JW, II) = ZO(JW, II) &
237                     + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
238                     + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
239                ! Momentum flux at launch lev
240                RUW0(JW, II) = RUWMAX
241             ENDDO
242          end DO
243       end DO
244    end DO
245
246    ! 4. COMPUTE THE FLUXES
247
248    ! 4.1 Vertical velocity at launching altitude to ensure
249    ! the correct value to the imposed fluxes.
250
251    DO JW = 1, NW
252
253       ! Evaluate intrinsic frequency at launching altitude:
254       ZOP(JW, :) = ZO(JW, :) &
255            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
256            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
257
258       ! VERSION WITH CONVECTIVE SOURCE
259
260       ! Vertical velocity at launch level, value to ensure the
261       ! imposed factor related to the convective forcing:
262       ! precipitations.
263
264       ! tanh limitation to values above prmax:
265       WWP(JW, :) = RUW0(JW, :) &
266            * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC / PRMAX))**2
267
268       ! Factor related to the characteristics of the waves:
269       WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW &
270            / MAX(ABS(ZOP(JW, :)), ZOISEC)**3
271
272       ! Moderation by the depth of the source (dz here):
273       WWP(JW, :) = WWP(JW, :) &
274            * EXP(- BVLOW**2 / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 * ZK(JW, :)**2 &
275            * DZ**2)
276
277       ! Put the stress in the right direction:
278       RUWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
279            * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
280       RVWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
281            * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
282    end DO
283
284    ! 4.2 Uniform values below the launching altitude
285
286    DO LL = 1, LAUNCH
287       RUW(:, LL) = 0
288       RVW(:, LL) = 0
289       DO JW = 1, NW
290          RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
291          RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
292       end DO
293    end DO
294
295    ! 4.3 Loop over altitudes, with passage from one level to the next
296    ! done by i) conserving the EP flux, ii) dissipating a little,
297    ! iii) testing critical levels, and vi) testing the breaking.
298
299    DO LL = LAUNCH, KLEV - 1
300       ! Warning: all the physics is here (passage from one level
301       ! to the next)
302       DO JW = 1, NW
303          ZOM(JW, :) = ZOP(JW, :)
304          WWM(JW, :) = WWP(JW, :)
305          ! Intrinsic Frequency
306          ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
307               - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
308
309          ! No breaking (Eq.6)
310          ! Dissipation (Eq. 8)
311          WWP(JW, :) = WWM(JW, :) * EXP(- 2. * RDISS * PR / (PH(:, LL + 1) &
312               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
313               / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
314               * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
315
316          ! Critical levels (forced to zero if intrinsic frequency changes sign)
317          ! Saturation (Eq. 12)
318          WWP(JW, :) = min(WWP(JW, :), MAX(0., &
319               SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
320               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SAT**2 * KMIN**2 &
321               / ZK(JW, :)**4)
322       end DO
323
324       ! Evaluate EP-flux from Eq. 7 and give the right orientation to
325       ! the stress
326
327       DO JW = 1, NW
328          RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
329          RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
330       end DO
331
332       RUW(:, LL + 1) = 0.
333       RVW(:, LL + 1) = 0.
334
335       DO JW = 1, NW
336          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
337          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
338       end DO
339    end DO
340
341    ! 5 CALCUL DES TENDANCES:
342
343    ! 5.1 Rectification des flux au sommet et dans les basses couches
344
345    RUW(:, KLEV + 1) = 0.
346    RVW(:, KLEV + 1) = 0.
347    RUW(:, 1) = RUW(:, LAUNCH)
348    RVW(:, 1) = RVW(:, LAUNCH)
349    DO LL = 1, LAUNCH
350       RUW(:, LL) = RUW(:, LAUNCH+1)
351       RVW(:, LL) = RVW(:, LAUNCH+1)
352    end DO
353
354    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
355    DO LL = 1, KLEV
356       D_U(:, LL) = (1.-DTIME/DELTAT) * D_U(:, LL) + DTIME/DELTAT/REAL(NW) * &
357            RG * (RUW(:, LL + 1) - RUW(:, LL)) &
358            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
359       D_V(:, LL) = (1.-DTIME/DELTAT) * D_V(:, LL) + DTIME/DELTAT/REAL(NW) * &
360            RG * (RVW(:, LL + 1) - RVW(:, LL)) &
361            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
362    ENDDO
363
364    ! Cosmetic: evaluation of the cumulated stress
365    ZUSTR = 0.
366    ZVSTR = 0.
367    DO LL = 1, KLEV
368       ZUSTR = ZUSTR + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
369       ZVSTR = ZVSTR + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
370    ENDDO
371
372  END SUBROUTINE FLOTT_GWD_RANDO
373
374end module FLOTT_GWD_rando_m
Note: See TracBrowser for help on using the repository browser.