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

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

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