source: LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.F90 @ 3613

Last change on this file since 3613 was 3531, checked in by Laurent Fairhead, 5 years ago

Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86
Still some work to be done in phylmd subdirectories
  • 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.
  • Property svn:keywords set to Id
File size: 15.2 KB
Line 
1!
2! $Id: flott_gwd_rando_m.F90 3531 2019-06-06 15:08:45Z acozic $
3!
4module FLOTT_GWD_rando_m
5
6  implicit none
7
8contains
9
10  SUBROUTINE FLOTT_GWD_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, &
11       d_v,east_gwstress,west_gwstress)
12
13    ! Parametrization of the momentum flux deposition due to a discrete
14    ! number of gravity waves.
15    ! Author: F. Lott
16    ! July, 12th, 2012
17    ! Gaussian distribution of the source, source is precipitation
18    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
19
20    !ONLINE:
21      use dimphy, only: klon, klev
22      use assert_m, only: assert
23      USE ioipsl_getin_p_mod, ONLY : getin_p
24      USE vertical_layers_mod, ONLY : presnivs
25      CHARACTER (LEN=20) :: modname='flott_gwd_rando'
26      CHARACTER (LEN=80) :: abort_message
27
28      include "YOMCST.h"
29      include "clesphys.h"
30    ! OFFLINE:
31    ! include "dimensions.h"
32    ! include "dimphy.h"
33    ! END OF DIFFERENCE ONLINE-OFFLINE
34    include "YOEGWD.h"
35
36    ! 0. DECLARATIONS:
37
38    ! 0.1 INPUTS
39    REAL, intent(in)::DTIME ! Time step of the Physics
40    REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels
41    REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s)
42    REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
43    REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
44    REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
45
46    ! 0.2 OUTPUTS
47    REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
48
49    REAL, intent(inout):: d_u(:, :), d_v(:, :)
50    REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
51    REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress
52
53    ! (KLON, KLEV) tendencies on winds
54
55    ! O.3 INTERNAL ARRAYS
56    REAL BVLOW(klon)
57    REAL DZ   !  Characteristic depth of the Source
58
59    INTEGER II, JJ, LL
60
61    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
62
63    REAL DELTAT
64
65    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
66
67    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
68    INTEGER JK, JP, JO, JW
69    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
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(NW, KLON) ! Horizontal wavenumber amplitude
75    REAL ZP(NW, KLON) ! Horizontal wavenumber angle
76    REAL ZO(NW, KLON) ! Absolute frequency !
77
78    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
79    REAL ZOM(NW, KLON), ZOP(NW, KLON)
80
81    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
82    REAL WWM(NW, KLON), WWP(NW, KLON)
83
84    REAL RUW0(NW, KLON) ! Fluxes at launching level
85
86    REAL RUWP(NW, KLON), RVWP(NW, KLON)
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    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
120    LOGICAL, SAVE :: firstcall = .TRUE.
121  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
122
123
124  IF (firstcall) THEN
125    ! Cle introduite pour resoudre un probleme de non reproductibilite
126    ! Le but est de pouvoir tester de revenir a la version precedenete
127    ! A eliminer rapidement
128    CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp)
129    IF (NW+3*NA>=KLEV) THEN
130       abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes'
131       CALL abort_physic (modname,abort_message,1)
132    ENDIF
133    firstcall=.false.
134  ENDIF
135
136
137    !-----------------------------------------------------------------
138
139    ! 1. INITIALISATIONS
140
141    ! 1.1 Basic parameter
142
143    ! Are provided from elsewhere (latent heat of vaporization, dry
144    ! gaz constant for air, gravity constant, heat capacity of dry air
145    ! at constant pressure, earth rotation rate, pi).
146
147    ! 1.2 Tuning parameters of V14
148
149   
150    RDISS = 0.5 ! Diffusion parameter
151    ! ONLINE
152      RUWMAX=GWD_RANDO_RUWMAX
153      SAT=gwd_rando_sat
154    !END ONLINE
155    ! OFFLINE
156    ! RUWMAX= 1.75    ! Launched flux
157    ! SAT=0.25     ! Saturation parameter
158    ! END OFFLINE
159
160    PRMAX = 20. / 24. /3600.
161    ! maximum of rain for which our theory applies (in kg/m^2/s)
162
163 ! Characteristic depth of the source
164    DZ = 1000.
165    XLAUNCH=0.5 ! Parameter that control launching altitude
166    XTROP=0.2 ! Parameter that control tropopause altitude
167    DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)
168    !  OFFLINE
169    !  DELTAT=DTIME
170    !  END OFFLINE
171
172    KMIN = 2.E-5
173    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
174
175    KMAX = 1.E-3 ! Max horizontal wavenumber
176    CMAX = 30. ! Max phase speed velocity
177
178    TR = 240. ! Reference Temperature
179    PR = 101300. ! Reference pressure
180    H0 = RD * TR / RG ! Characteristic vertical scale height
181
182    BVSEC = 5.E-3 ! Security to avoid negative BVF
183    PSEC = 1.E-6 ! Security to avoid division by 0 pressure
184    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
185
186IF (1==0) THEN
187    !ONLINE
188        call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
189         size(vv, 1), size(zustr), size(zvstr), size(d_u, 1), &
190         size(d_v, 1), &
191         size(east_gwstress, 1), size(west_gwstress, 1) /), &
192         "FLOTT_GWD_RANDO klon")
193     call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
194          size(vv, 2), size(d_u, 2), size(d_v, 2), &
195          size(east_gwstress,2), size(west_gwstress,2) /), &
196          "FLOTT_GWD_RANDO klev")
197    !END ONLINE
198ENDIF
199
200    IF(DELTAT < DTIME)THEN
201       abort_message='flott_gwd_rando: deltat < dtime!'
202       CALL abort_physic(modname,abort_message,1)
203    ENDIF
204
205    IF (KLEV < NW) THEN
206       abort_message='flott_gwd_rando: you will have problem with random numbers'
207       CALL abort_physic(modname,abort_message,1)
208    ENDIF
209
210    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
211
212    ! Pressure and Inv of pressure
213    DO LL = 2, KLEV
214       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
215    end DO
216    PH(:, KLEV + 1) = 0.
217    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
218
219    ! Launching altitude
220
221    !Pour revenir a la version non reproductible en changeant le nombre de process
222    IF (gwd_reproductibilite_mpiomp) THEN
223       ! Reprend la formule qui calcule PH en fonction de PP=play
224       DO LL = 2, KLEV
225          HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.)
226       end DO
227       HREF(KLEV + 1) = 0.
228       HREF(1) = 2. * presnivs(1) - HREF(2)
229    ELSE
230       HREF(1:KLEV)=PH(KLON/2,1:KLEV)
231    ENDIF
232
233    LAUNCH=0
234    LTROP =0
235    DO LL = 1, KLEV
236       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
237    ENDDO
238    DO LL = 1, KLEV
239       IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL
240    ENDDO
241    !LAUNCH=22 ; LTROP=33
242!   print*,'LAUNCH=',LAUNCH,'LTROP=',LTROP
243
244    ! Log pressure vert. coordinate
245    DO LL = 1, KLEV + 1
246       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
247    end DO
248
249    ! BV frequency
250    DO LL = 2, KLEV
251       ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
252       UH(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
253            * RD**2 / RCPD / H0**2 + (TT(:, LL) &
254            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
255    end DO
256    BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
257         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
258         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
259
260    UH(:, 1) = UH(:, 2)
261    UH(:, KLEV + 1) = UH(:, KLEV)
262    BV(:, 1) = UH(:, 2)
263    BV(:, KLEV + 1) = UH(:, KLEV)
264    ! SMOOTHING THE BV HELPS
265    DO LL = 2, KLEV
266       BV(:, LL)=(UH(:, LL+1)+2.*UH(:, LL)+UH(:, LL-1))/4.
267    end DO
268
269    BV=MAX(SQRT(MAX(BV, 0.)), BVSEC)
270    BVLOW=MAX(SQRT(MAX(BVLOW, 0.)), BVSEC)
271
272
273    ! WINDS
274    DO LL = 2, KLEV
275       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
276       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
277    end DO
278    UH(:, 1) = 0.
279    VH(:, 1) = 0.
280    UH(:, KLEV + 1) = UU(:, KLEV)
281    VH(:, KLEV + 1) = VV(:, KLEV)
282
283    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE
284
285    ! The mod functions of weird arguments are used to produce the
286    ! waves characteristics in an almost stochastic way
287
288    DO JW = 1, NW
289             ! Angle
290             DO II = 1, KLON
291                ! Angle (0 or PI so far)
292                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
293                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
294                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
295                     * RPI / 2.
296                ! Horizontal wavenumber amplitude
297                ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
298                ! Horizontal phase speed
299                CPHA = 0.
300                DO JJ = 1, NA
301                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
302                    CPHA = CPHA + &
303                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
304                END DO
305                IF (CPHA.LT.0.)  THEN
306                   CPHA = -1.*CPHA
307                   ZP(JW,II) = ZP(JW,II) + RPI
308                ENDIF
309                ! Absolute frequency is imposed
310                ZO(JW, II) = CPHA * ZK(JW, II)
311                ! Intrinsic frequency is imposed
312                ZO(JW, II) = ZO(JW, II) &
313                     + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
314                     + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
315                ! Momentum flux at launch lev
316                RUW0(JW, II) = RUWMAX
317             ENDDO
318    ENDDO
319
320    ! 4. COMPUTE THE FLUXES
321
322    ! 4.1 Vertical velocity at launching altitude to ensure
323    ! the correct value to the imposed fluxes.
324
325    DO JW = 1, NW
326
327       ! Evaluate intrinsic frequency at launching altitude:
328       ZOP(JW, :) = ZO(JW, :) &
329            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
330            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
331
332       ! VERSION WITH CONVECTIVE SOURCE
333
334       ! Vertical velocity at launch level, value to ensure the
335       ! imposed factor related to the convective forcing:
336       ! precipitations.
337
338       ! tanh limitation to values above prmax:
339       WWP(JW, :) = RUW0(JW, :) &
340            * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
341
342       ! Factor related to the characteristics of the waves:
343       WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:)  &
344            / MAX(ABS(ZOP(JW, :)), ZOISEC)**3
345
346       ! Moderation by the depth of the source (dz here):
347       WWP(JW, :) = WWP(JW, :) &
348            * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 * ZK(JW, :)**2 &
349            * DZ**2)
350
351       ! Put the stress in the right direction:
352       RUWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
353            * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
354       RVWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
355            * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
356    end DO
357
358
359    ! 4.2 Uniform values below the launching altitude
360
361    DO LL = 1, LAUNCH
362       RUW(:, LL) = 0
363       RVW(:, LL) = 0
364       DO JW = 1, NW
365          RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
366          RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
367       end DO
368    end DO
369
370    ! 4.3 Loop over altitudes, with passage from one level to the next
371    ! done by i) conserving the EP flux, ii) dissipating a little,
372    ! iii) testing critical levels, and vi) testing the breaking.
373
374    DO LL = LAUNCH, KLEV - 1
375       ! Warning: all the physics is here (passage from one level
376       ! to the next)
377       DO JW = 1, NW
378          ZOM(JW, :) = ZOP(JW, :)
379          WWM(JW, :) = WWP(JW, :)
380          ! Intrinsic Frequency
381          ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
382               - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
383
384          ! No breaking (Eq.6)
385          ! Dissipation (Eq. 8)
386          WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
387               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
388               / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
389               * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
390
391          ! Critical levels (forced to zero if intrinsic frequency changes sign)
392          ! Saturation (Eq. 12)
393          WWP(JW, :) = min(WWP(JW, :), MAX(0., &
394               SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
395               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2  &
396               * SAT**2 / ZK(JW, :)**4)
397       end DO
398
399       ! Evaluate EP-flux from Eq. 7 and give the right orientation to
400       ! the stress
401
402       DO JW = 1, NW
403          RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
404          RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
405       end DO
406
407       RUW(:, LL + 1) = 0.
408       RVW(:, LL + 1) = 0.
409
410       DO JW = 1, NW
411          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
412          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
413          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW)
414          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW)
415       end DO
416    end DO
417! OFFLINE ONLY
418!   PRINT *,'SAT PROFILE:'
419!   DO LL=1,KLEV
420!   PRINT *,ZH(KLON/2,LL)/1000.,SAT*(2.+TANH(ZH(KLON/2,LL)/H0-8.))
421!   ENDDO
422
423    ! 5 CALCUL DES TENDANCES:
424
425    ! 5.1 Rectification des flux au sommet et dans les basses couches
426
427    RUW(:, KLEV + 1) = 0.
428    RVW(:, KLEV + 1) = 0.
429    RUW(:, 1) = RUW(:, LAUNCH)
430    RVW(:, 1) = RVW(:, LAUNCH)
431    DO LL = 1, LAUNCH
432       RUW(:, LL) = RUW(:, LAUNCH+1)
433       RVW(:, LL) = RVW(:, LAUNCH+1)
434       EAST_GWSTRESS(:, LL)  = EAST_GWSTRESS(:, LAUNCH)
435       WEST_GWSTRESS(:, LL)  = WEST_GWSTRESS(:, LAUNCH)
436    end DO
437
438    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
439    DO LL = 1, KLEV
440       D_U(:, LL) = (1.-DTIME/DELTAT) * D_U(:, LL) + DTIME/DELTAT/REAL(NW) * &
441            RG * (RUW(:, LL + 1) - RUW(:, LL)) &
442            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
443       ! NO AR-1 FOR MERIDIONAL TENDENCIES
444       D_V(:, LL) =                                            1./REAL(NW) * &
445            RG * (RVW(:, LL + 1) - RVW(:, LL)) &
446            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
447    ENDDO
448
449    ! Cosmetic: evaluation of the cumulated stress
450    ZUSTR = 0.
451    ZVSTR = 0.
452    DO LL = 1, KLEV
453       ZUSTR = ZUSTR + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
454       ZVSTR = ZVSTR + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
455    ENDDO
456
457
458  END SUBROUTINE FLOTT_GWD_RANDO
459
460end module FLOTT_GWD_rando_m
Note: See TracBrowser for help on using the repository browser.