source: LMDZ6/branches/LMDZ-QUEST/libf/phylmd/flott_gwd_rando_m.F90

Last change on this file was 3200, checked in by Laurent Fairhead, 7 years ago

Inclusion of r3198 from trunk
Retour vers l'insensibilite au decoupage en sous domaine.
Les routines gwd_rando incluait le calcul de niveaux de reference
sur la base d'un profile pris au milieu du domaine (en klon/2).
Rempace par un test en presnivs.

Une autre intercation entre routines concernant la tke a fait apparaitre
que la tke n'?\195?\169tait pas passee correctement au niveau klev+1 au moment
du regroupement des mailles sous les sous surface.

Ces changements garantissent la convergence numerique si
addtkeoro=0
iflag_pbl<12
et
ok_gwd_rando=n
La convergence n'est pas garantie pour les dernieres versions des physiq.def
mais les differences devraient etre mineures.

FH
  • 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.
File size: 15.0 KB
Line 
1module FLOTT_GWD_rando_m
2
3  implicit none
4
5contains
6
7  SUBROUTINE FLOTT_GWD_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, &
8       d_v,east_gwstress,west_gwstress)
9
10    ! Parametrization of the momentum flux deposition due to a discrete
11    ! number of gravity waves.
12    ! Author: F. Lott
13    ! July, 12th, 2012
14    ! Gaussian distribution of the source, source is precipitation
15    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
16
17    !ONLINE:
18      use dimphy, only: klon, klev
19      use assert_m, only: assert
20      USE ioipsl_getin_p_mod, ONLY : getin_p
21      USE vertical_layers_mod, ONLY : presnivs
22
23      include "YOMCST.h"
24      include "clesphys.h"
25    ! OFFLINE:
26    ! include "dimensions.h"
27    ! include "dimphy.h"
28    ! END OF DIFFERENCE ONLINE-OFFLINE
29    include "YOEGWD.h"
30
31    ! 0. DECLARATIONS:
32
33    ! 0.1 INPUTS
34    REAL, intent(in)::DTIME ! Time step of the Physics
35    REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels
36    REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s)
37    REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
38    REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
39    REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
40
41    ! 0.2 OUTPUTS
42    REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
43
44    REAL, intent(inout):: d_u(:, :), d_v(:, :)
45    REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
46    REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress
47
48    ! (KLON, KLEV) tendencies on winds
49
50    ! O.3 INTERNAL ARRAYS
51    REAL BVLOW(klon)
52    REAL DZ   !  Characteristic depth of the Source
53
54    INTEGER II, JJ, LL
55
56    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
57
58    REAL DELTAT
59
60    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
61
62    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
63    INTEGER JK, JP, JO, JW
64    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
65    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
66    REAL CMAX ! standard deviation of the phase speed distribution
67    REAL RUWMAX,SAT  ! ONLINE SPECIFIED IN run.def
68    REAL CPHA ! absolute PHASE VELOCITY frequency
69    REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
70    REAL ZP(NW, KLON) ! Horizontal wavenumber angle
71    REAL ZO(NW, KLON) ! Absolute frequency !
72
73    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
74    REAL ZOM(NW, KLON), ZOP(NW, KLON)
75
76    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
77    REAL WWM(NW, KLON), WWP(NW, KLON)
78
79    REAL RUW0(NW, KLON) ! Fluxes at launching level
80
81    REAL RUWP(NW, KLON), RVWP(NW, KLON)
82    ! Fluxes X and Y for each waves at 1/2 Levels
83
84    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
85
86    REAL XLAUNCH ! Controle the launching altitude
87    REAL XTROP ! SORT of Tropopause altitude
88    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
89    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
90
91    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
92    ! for GWs parameterisation apply
93
94    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
95
96    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
97
98    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
99
100    REAL H0 ! Characteristic Height of the atmosphere
101    REAL PR, TR ! Reference Pressure and Temperature
102
103    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
104
105    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
106    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
107    REAL PSEC ! Security to avoid division by 0 pressure
108    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
109    REAL BVSEC ! Security to avoid negative BVF
110    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
111
112    REAL, DIMENSION(klev+1) ::HREF
113
114    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
115    LOGICAL, SAVE :: firstcall = .TRUE.
116  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
117
118    CHARACTER (LEN=20) :: modname='flott_gwd_rando'
119    CHARACTER (LEN=80) :: abort_message
120
121
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   
149    RDISS = 0.5 ! Diffusion parameter
150    ! ONLINE
151      RUWMAX=GWD_RANDO_RUWMAX
152      SAT=gwd_rando_sat
153    !END ONLINE
154    ! OFFLINE
155    ! RUWMAX= 1.75    ! Launched flux
156    ! SAT=0.25     ! Saturation parameter
157    ! END OFFLINE
158
159    PRMAX = 20. / 24. /3600.
160    ! maximum of rain for which our theory applies (in kg/m^2/s)
161
162 ! Characteristic depth of the source
163    DZ = 1000.
164    XLAUNCH=0.5 ! Parameter that control launching altitude
165    XTROP=0.2 ! Parameter that control tropopause altitude
166    DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)
167    !  OFFLINE
168    !  DELTAT=DTIME
169    !  END OFFLINE
170
171    KMIN = 2.E-5
172    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
173
174    KMAX = 1.E-3 ! Max horizontal wavenumber
175    CMAX = 30. ! Max phase speed velocity
176
177    TR = 240. ! Reference Temperature
178    PR = 101300. ! Reference pressure
179    H0 = RD * TR / RG ! Characteristic vertical scale height
180
181    BVSEC = 5.E-3 ! Security to avoid negative BVF
182    PSEC = 1.E-6 ! Security to avoid division by 0 pressure
183    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
184
185IF (1==0) THEN
186    !ONLINE
187        call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
188         size(vv, 1), size(zustr), size(zvstr), size(d_u, 1), &
189         size(d_v, 1), &
190         size(east_gwstress, 1), size(west_gwstress, 1) /), &
191         "FLOTT_GWD_RANDO klon")
192     call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
193          size(vv, 2), size(d_u, 2), size(d_v, 2), &
194          size(east_gwstress,2), size(west_gwstress,2) /), &
195          "FLOTT_GWD_RANDO klev")
196    !END ONLINE
197ENDIF
198
199    IF(DELTAT < DTIME)THEN
200       PRINT *, 'flott_gwd_rando: deltat < dtime!'
201       STOP 1
202    ENDIF
203
204    IF (KLEV < NW) THEN
205       PRINT *, 'flott_gwd_rando: you will have problem with random numbers'
206       STOP 1
207    ENDIF
208
209    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
210
211    ! Pressure and Inv of pressure
212    DO LL = 2, KLEV
213       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
214    end DO
215    PH(:, KLEV + 1) = 0.
216    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
217
218    ! Launching altitude
219
220    !Pour revenir a la version non reproductible en changeant le nombre de process
221    IF (gwd_reproductibilite_mpiomp) THEN
222       ! Reprend la formule qui calcule PH en fonction de PP=play
223       DO LL = 2, KLEV
224          HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.)
225       end DO
226       HREF(KLEV + 1) = 0.
227       HREF(1) = 2. * presnivs(1) - HREF(2)
228    ELSE
229       HREF(1:KLEV)=PH(KLON/2,1:KLEV)
230    ENDIF
231
232    LAUNCH=0
233    LTROP =0
234    DO LL = 1, KLEV
235       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
236    ENDDO
237    DO LL = 1, KLEV
238       IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL
239    ENDDO
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.LT.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
455  END SUBROUTINE FLOTT_GWD_RANDO
456
457end module FLOTT_GWD_rando_m
Note: See TracBrowser for help on using the repository browser.