source: LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90

Last change on this file was 5274, checked in by abarral, 49 minutes ago

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