source: LMDZ5/trunk/libf/phylmd/flott_gwd_rando_m.F90 @ 4984

Last change on this file since 4984 was 2665, checked in by dcugnet, 8 years ago
  • A (re)startphy.nc file (standard name: "startphy0.nc") can be read by ce0l to get land mask, so mask can be defined (in decreasing priority order) from: 1) "o2a.nc file" if this file is found 2) "startphy0.nc" if this file is found 3) "Relief.nc" otherwise
  • Sub-cell scales parameters for orographic gravity waves can be read from file "oro_params.nc" if the configuration key "read_orop" is TRUE. The effect is to bypass the "grid_noro" routine in ce0l, so that any pre-defined mask (from o2a.nc or startphy0.nc) is then overwritten.
  • The gcm stops if the "limit.nc" records number differs from the current year number of days. A warning is issued in case the gcm calendar does not match the time axis attribute "calendar" (if available) from the "limit.nc" file. This attribute is now added to the "limit.nc" time axis.
  • Few simplifications in grid_noro
  • Few parameters changes in acama_gwd and flott_gwd.
  • Variable d_u can be saved in the outputs.
  • 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.