source: LMDZ6/branches/Amaury_dev/libf/phylmd/flott_gwd_rando_m.F90 @ 5135

Last change on this file since 5135 was 5119, checked in by abarral, 16 months ago

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