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

Last change on this file since 5119 was 5119, checked in by abarral, 4 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.