source: LMDZ5/branches/testing/libf/phylmd/flott_gwd_rando_m.F90 @ 3525

Last change on this file since 3525 was 2720, checked in by Laurent Fairhead, 8 years ago

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