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

Last change on this file since 5144 was 5144, checked in by abarral, 7 weeks ago

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