source: LMDZ6/branches/Amaury_dev/libf/phylmd/acama_gwd_rando_m.F90 @ 5159

Last change on this file since 5159 was 5159, checked in by abarral, 3 months ago

Put dimensions.h and paramet.h into modules
  • Property svn:keywords set to Id
File size: 19.1 KB
Line 
1! $Id: acama_gwd_rando_m.F90 5159 2024-08-02 19:58:25Z abarral $
2
3module ACAMA_GWD_rando_m
4
5  IMPLICIT NONE
6
7CONTAINS
8
9  SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
10          zustr, zvstr, d_u, 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, A. de la Camara
15    ! July, 24th, 2014
16    ! Gaussian distribution of the source, source is vorticity squared
17    ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
18    ! Lott et al (JAS, 2010, vol 67, page 157-170)
19    ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
20
21    !  ONLINE:
22    USE dimphy, ONLY: klon, klev
23    USE lmdz_assert, ONLY: assert
24    USE lmdz_ioipsl_getin_p, ONLY: getin_p
25    USE lmdz_vertical_layers, ONLY: presnivs
26    USE lmdz_abort_physic, ONLY: abort_physic
27    USE lmdz_clesphys
28    USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, &
29            GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat
30    USE lmdz_yomcst
31
32    ! 0. DECLARATIONS:
33
34    ! 0.1 INPUTS
35    REAL, INTENT(IN) :: DTIME ! Time step of the Physics
36    REAL, INTENT(IN) :: PP(:, :) ! (KLON, KLEV) Pressure at full levels
37    REAL, INTENT(IN) :: ROT(:, :) ! Relative vorticity
38    REAL, INTENT(IN) :: TT(:, :) ! (KLON, KLEV) Temp at full levels
39    REAL, INTENT(IN) :: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
40    REAL, INTENT(IN) :: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
41    REAL, INTENT(IN) :: PLAT(:) ! (KLON) LATITUDE
42
43    ! 0.2 OUTPUTS
44    REAL, INTENT(OUT) :: zustr(:), zvstr(:) ! (KLON) Surface Stresses
45
46    REAL, INTENT(INOUT) :: d_u(:, :), d_v(:, :)
47    REAL, INTENT(INOUT) :: east_gwstress(:, :) !  Profile of eastward stress
48    REAL, INTENT(INOUT) :: west_gwstress(:, :) !  Profile of westward stress
49    ! (KLON, KLEV) tendencies on winds
50
51    ! O.3 INTERNAL ARRAYS
52    REAL BVLOW(klon)  !  LOW LEVEL BV FREQUENCY
53    REAL ROTBA(KLON), CORIO(KLON)  !  BAROTROPIC REL. VORTICITY AND PLANETARY
54    REAL UZ(KLON, KLEV + 1)
55
56    INTEGER II, JJ, LL
57
58    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
59
60    REAL DELTAT
61
62    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
63
64    INTEGER, PARAMETER :: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
65    INTEGER JK, JP, JO, JW
66    INTEGER, PARAMETER :: NA = 5  !number of realizations to get the phase speed
67    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
68    REAL CMIN, CMAX ! Min and Max absolute ph. vel.
69    REAL CPHA ! absolute PHASE VELOCITY frequency
70    REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
71    REAL ZP(NW, KLON) ! Horizontal wavenumber angle
72    REAL ZO(NW, KLON) ! Absolute frequency !
73
74    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
75    REAL ZOM(NW, KLON), ZOP(NW, KLON)
76
77    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
78    REAL WWM(NW, KLON), WWP(NW, KLON)
79
80    REAL RUW0(NW, KLON) ! Fluxes at launching level
81
82    REAL RUWP(NW, KLON), RVWP(NW, KLON)
83    ! Fluxes X and Y for each waves at 1/2 Levels
84
85    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
86
87    REAL XLAUNCH ! Controle the launching altitude
88    REAL XTROP ! SORT of Tropopause altitude
89    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
90    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
91
92    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
93    ! for GWs parameterisation apply
94
95    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
96
97    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
98    REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS
99    REAL RUWFRT, SATFRT
100
101    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
102
103    REAL H0 ! Characteristic Height of the atmosphere
104    REAL DZ ! Characteristic depth of the source!
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 PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
113    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
114    REAL BVSEC ! Security to avoid negative BVF
115
116    REAL, DIMENSION(klev + 1) :: HREF
117    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp = .TRUE.
118    LOGICAL, SAVE :: firstCALL = .TRUE.
119    !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
120
121    CHARACTER (LEN = 20) :: modname = 'acama_gwd_rando_m'
122    CHARACTER (LEN = 80) :: abort_message
123
124    IF (firstcall) THEN
125      ! Cle introduite pour resoudre un probleme de non reproductibilite
126      ! Le but est de pouvoir tester de revenir a la version precedenete
127      ! A eliminer rapidement
128      CALL getin_p('gwd_reproductibilite_mpiomp', gwd_reproductibilite_mpiomp)
129      IF (NW + 4 * (NA - 1) + NA>=KLEV) THEN
130        abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes'
131        CALL abort_physic (modname, abort_message, 1)
132      ENDIF
133      firstcall = .FALSE.
134      !    CALL iophys_ini(dtime)
135    ENDIF
136
137    !-----------------------------------------------------------------
138
139    ! 1. INITIALISATIONS
140
141    ! 1.1 Basic parameter
142
143    ! Are provided from elsewhere (latent heat of vaporization, dry
144    ! gaz constant for air, gravity constant, heat capacity of dry air
145    ! at constant pressure, earth rotation rate, pi).
146
147    ! 1.2 Tuning parameters of V14
148
149    ! Values for linear in rot (recommended):
150    !   RUWFRT=0.005 ! As RUWMAX but for frontal waves
151    !   SATFRT=1.00  ! As SAT    but for frontal waves
152    ! Values when rot^2 is used
153    !    RUWFRT=0.02  ! As RUWMAX but for frontal waves
154    !    SATFRT=1.00  ! As SAT    but for frontal waves
155    !    CMAX = 30.   ! Characteristic phase speed
156    ! Values when rot^2*EXP(-pi*sqrt(J)) is used
157    !   RUWFRT=2.5   ! As RUWMAX but for frontal waves ~ N0*F0/4*DZ
158    !   SATFRT=0.60   ! As SAT    but for frontal waves
159    RUWFRT = gwd_front_ruwmax
160    SATFRT = gwd_front_sat
161    CMAX = 50.    ! Characteristic phase speed
162    ! Phase speed test
163    !   RUWFRT=0.01
164    !   CMAX = 50.   ! Characteristic phase speed (TEST)
165    ! Values when rot^2 and exp(-m^2*dz^2) are used
166    !   RUWFRT=0.03  ! As RUWMAX but for frontal waves
167    !   SATFRT=1.00  ! As SAT    but for frontal waves
168    ! CRUCIAL PARAMETERS FOR THE WIND FILTERING
169    XLAUNCH = 0.95 ! Parameter that control launching altitude
170    RDISS = 0.5  ! Diffusion parameter
171
172    ! maximum of rain for which our theory applies (in kg/m^2/s)
173
174    DZ = 1000. ! Characteristic depth of the source
175    XTROP = 0.2 ! Parameter that control tropopause altitude
176    DELTAT = 24. * 3600. ! Time scale of the waves (first introduced in 9b)
177    !   DELTAT=DTIME     ! No AR-1 Accumulation, OR OFFLINE
178
179    KMIN = 2.E-5
180    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
181
182    KMAX = 1.E-3  ! Max horizontal wavenumber
183    CMIN = 1.     ! Min phase velocity
184
185    TR = 240. ! Reference Temperature
186    PR = 101300. ! Reference pressure
187    H0 = RD * TR / RG ! Characteristic vertical scale height
188
189    BVSEC = 5.E-3 ! Security to avoid negative BVF
190    PSEC = 1.E-6 ! Security to avoid division by 0 pressure
191    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
192    CORSEC = ROMEGA * 2. * SIN(2. * RPI / 180.)! Security for CORIO
193
194    !  ONLINE
195    CALL assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
196            size(vv, 1), size(rot, 1), size(zustr), size(zvstr), size(d_u, 1), &
197            size(d_v, 1), &
198            size(east_gwstress, 1), size(west_gwstress, 1) /), &
199            "ACAMA_GWD_RANDO klon")
200    CALL assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
201            size(vv, 2), size(d_u, 2), size(d_v, 2), &
202            size(east_gwstress, 2), size(west_gwstress, 2) /), &
203            "ACAMA_GWD_RANDO klev")
204    !  END ONLINE
205
206    IF(DELTAT < DTIME)THEN
207      !       PRINT *, 'flott_gwd_rando: deltat < dtime!'
208      !       STOP 1
209      abort_message = ' deltat < dtime! '
210      CALL abort_physic(modname, abort_message, 1)
211    ENDIF
212
213    IF (KLEV < NW) THEN
214      !       PRINT *, 'flott_gwd_rando: you will have problem with random numbers'
215      !       STOP 1
216      abort_message = ' you will have problem with random numbers'
217      CALL abort_physic(modname, abort_message, 1)
218    ENDIF
219
220    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
221
222    ! Pressure and Inv of pressure
223    DO LL = 2, KLEV
224      PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
225      PHM1(:, LL) = 1. / PH(:, LL)
226    end DO
227
228    PH(:, KLEV + 1) = 0.
229    PHM1(:, KLEV + 1) = 1. / PSEC
230    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
231
232    ! Launching altitude
233
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
257    !   PRINT *,'LAUNCH IN ACAMARA:',LAUNCH
258
259    ! Log pressure vert. coordinate
260    DO LL = 1, KLEV + 1
261      ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
262    end DO
263
264    ! BV frequency
265    DO LL = 2, KLEV
266      ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
267      UH(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
268              * RD**2 / RCPD / H0**2 + (TT(:, LL) &
269              - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
270    end DO
271    BVLOW = 0.5 * (TT(:, LTROP) + TT(:, LAUNCH)) &
272            * RD**2 / RCPD / H0**2 + (TT(:, LTROP) &
273            - TT(:, LAUNCH)) / (ZH(:, LTROP) - ZH(:, LAUNCH)) * RD / H0
274
275    UH(:, 1) = UH(:, 2)
276    UH(:, KLEV + 1) = UH(:, KLEV)
277    BV(:, 1) = UH(:, 2)
278    BV(:, KLEV + 1) = UH(:, KLEV)
279    ! SMOOTHING THE BV HELPS
280    DO LL = 2, KLEV
281      BV(:, LL) = (UH(:, LL + 1) + 2. * UH(:, LL) + UH(:, LL - 1)) / 4.
282    end DO
283
284    BV = MAX(SQRT(MAX(BV, 0.)), BVSEC)
285    BVLOW = MAX(SQRT(MAX(BVLOW, 0.)), BVSEC)
286
287    ! WINDS
288    DO LL = 2, KLEV
289      UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
290      VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
291      UZ(:, LL) = ABS((SQRT(UU(:, LL)**2 + VV(:, LL)**2) &
292              - SQRT(UU(:, LL - 1)**2 + VV(:, LL - 1)**2)) &
293              / (ZH(:, LL) - ZH(:, LL - 1)))
294    end DO
295    UH(:, 1) = 0.
296    VH(:, 1) = 0.
297    UH(:, KLEV + 1) = UU(:, KLEV)
298    VH(:, KLEV + 1) = VV(:, KLEV)
299
300    UZ(:, 1) = UZ(:, 2)
301    UZ(:, KLEV + 1) = UZ(:, KLEV)
302    UZ(:, :) = MAX(UZ(:, :), PSEC)
303
304    ! BAROTROPIC VORTICITY AND INTEGRATED CORIOLIS PARAMETER
305
306    CORIO(:) = MAX(ROMEGA * 2. * ABS(SIN(PLAT(:) * RPI / 180.)), CORSEC)
307    ROTBA(:) = 0.
308    DO LL = 1, KLEV - 1
309      !ROTBA(:) = ROTBA(:) + (ROT(:,LL)+ROT(:,LL+1))/2./RG*(PP(:,LL)-PP(:,LL+1))
310      ! Introducing the complete formula (exp of Richardson number):
311      ROTBA(:) = ROTBA(:) + &
312              !((ROT(:,LL)+ROT(:,LL+1))/2.)**2 &
313              (CORIO(:) * TANH(ABS(ROT(:, LL) + ROT(:, LL + 1)) / 2. / CORIO(:)))**2 &
314                      / RG * (PP(:, LL) - PP(:, LL + 1)) &
315                      * EXP(-RPI * BV(:, LL + 1) / UZ(:, LL + 1)) &
316                      !               * DZ*BV(:,LL+1)/4./ABS(CORIO(:))
317                      * DZ * BV(:, LL + 1) / 4. / 1.E-4           !  Changes after 1991
318      !ARRET
319    ENDDO
320    !   PRINT *,'MAX ROTBA:',MAXVAL(ROTBA)
321    !   ROTBA(:)=(1.*ROTBA(:)  & ! Testing zone
322    !           +0.15*CORIO(:)**2                &
323    !           /(COS(PLAT(:)*RPI/180.)+0.02)  &
324    !           )*DZ*0.01/0.0001/4. ! & ! Testing zone
325    !   MODIF GWD4 AFTER 1985
326    !            *(1.25+SIN(PLAT(:)*RPI/180.))/(1.05+SIN(PLAT(:)*RPI/180.))/1.25
327    !          *1./(COS(PLAT(:)*RPI/180.)+0.02)
328    !    CORIO(:) = MAX(ROMEGA*2.*ABS(SIN(PLAT(:)*RPI/180.)),ZOISEC)/RG*PP(:,1)
329
330    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE
331
332    ! The mod functions of weird arguments are used to produce the
333    ! waves characteristics in an almost stochastic way
334
335    JW = 0
336    DO JW = 1, NW
337      ! Angle
338      DO II = 1, KLON
339        ! Angle (0 or PI so far)
340        ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
341        !      * RPI / 2.
342        ! Angle between 0 and pi
343        ZP(JW, II) = MOD(TT(II, JW) * 10., 1.) * RPI
344        ! TEST WITH POSITIVE WAVES ONLY (Part I/II)
345        !               ZP(JW, II) = 0.
346        ! Horizontal wavenumber amplitude
347        ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
348        ! Horizontal phase speed
349        CPHA = 0.
350        DO JJ = 1, NA
351          CPHA = CPHA + &
352                  CMAX * 2. * (MOD(TT(II, JW + 4 * (JJ - 1) + JJ)**2, 1.) - 0.5) * SQRT(3.) / SQRT(NA * 1.)
353        END DO
354        IF (CPHA<0.)  THEN
355          CPHA = -1. * CPHA
356          ZP(JW, II) = ZP(JW, II) + RPI
357          ! TEST WITH POSITIVE WAVES ONLY (Part II/II)
358          !               ZP(JW, II) = 0.
359        ENDIF
360        CPHA = CPHA + CMIN !we dont allow |c|<1m/s
361        ! Absolute frequency is imposed
362        ZO(JW, II) = CPHA * ZK(JW, II)
363        ! Intrinsic frequency is imposed
364        ZO(JW, II) = ZO(JW, II) &
365                + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
366                + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
367        ! Momentum flux at launch lev
368        ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE
369        !  RIGHT IN THE SH (GWD4 after 1990)
370        RUW0(JW, II) = 0.
371        DO JJ = 1, NA
372          RUW0(JW, II) = RUW0(JW, II) + &
373                  2. * (MOD(TT(II, JW + 4 * (JJ - 1) + JJ)**2, 1.) - 0.5) * SQRT(3.) / SQRT(NA * 1.)
374        END DO
375        RUW0(JW, II) = RUWFRT &
376                * EXP(RUW0(JW, II)) / 1250. &  ! 2 mpa at south pole
377                * ((1.05 + SIN(PLAT(II) * RPI / 180.)) / (1.01 + SIN(PLAT(II) * RPI / 180.)) - 2.05 / 2.01)
378        ! RUW0(JW, II) = RUWFRT
379      ENDDO
380    end DO
381
382    ! 4. COMPUTE THE FLUXES
383
384    ! 4.0
385
386    ! 4.1 Vertical velocity at launching altitude to ensure
387    ! the correct value to the imposed fluxes.
388
389    DO JW = 1, NW
390
391      ! Evaluate intrinsic frequency at launching altitude:
392      ZOP(JW, :) = ZO(JW, :) &
393              - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
394              - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
395
396      ! VERSION WITH FRONTAL SOURCES
397
398      ! Momentum flux at launch level imposed by vorticity sources
399
400      ! tanh limitation for values above CORIO (inertial instability).
401      ! WWP(JW, :) = RUW0(JW, :) &
402      WWP(JW, :) = RUWFRT      &
403              !     * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 &
404              !    * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) &
405              !  CONSTANT FLUX
406              !    * (CORIO(:)*CORIO(:)) &
407              ! MODERATION BY THE DEPTH OF THE SOURCE (DZ HERE)
408              !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(JW, :)),ZOISEC)**2 &
409              !      *ZK(JW, :)**2*DZ**2) &
410              ! COMPLETE FORMULA:
411              !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) &
412              * ROTBA(:) &
413              !  RESTORE DIMENSION OF A FLUX
414              !     *RD*TR/PR
415              !     *1. + RUW0(JW, :)
416              * 1.
417
418      ! Factor related to the characteristics of the waves: NONE
419
420      ! Moderation by the depth of the source (dz here): NONE
421
422      ! Put the stress in the right direction:
423
424      RUWP(JW, :) = SIGN(1., ZOP(JW, :)) * COS(ZP(JW, :)) * WWP(JW, :)
425      RVWP(JW, :) = SIGN(1., ZOP(JW, :)) * SIN(ZP(JW, :)) * WWP(JW, :)
426
427    end DO
428
429    ! 4.2 Uniform values below the launching altitude
430
431    DO LL = 1, LAUNCH
432      RUW(:, LL) = 0
433      RVW(:, LL) = 0
434      DO JW = 1, NW
435        RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
436        RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
437      end DO
438    end DO
439
440    ! 4.3 Loop over altitudes, with passage from one level to the next
441    ! done by i) conserving the EP flux, ii) dissipating a little,
442    ! iii) testing critical levels, and vi) testing the breaking.
443
444    DO LL = LAUNCH, KLEV - 1
445      ! Warning: all the physics is here (passage from one level
446      ! to the next)
447      DO JW = 1, NW
448        ZOM(JW, :) = ZOP(JW, :)
449        WWM(JW, :) = WWP(JW, :)
450        ! Intrinsic Frequency
451        ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
452                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
453
454        ! No breaking (Eq.6)
455        ! Dissipation (Eq. 8)
456        WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
457                + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
458                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
459                * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
460
461        ! Critical levels (forced to zero if intrinsic frequency changes sign)
462        ! Saturation (Eq. 12)
463        WWP(JW, :) = min(WWP(JW, :), MAX(0., &
464                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
465                !    / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 &
466                / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 &
467                !              *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 &
468                * SATFRT**2       &
469                / ZK(JW, :)**4)
470      end DO
471
472      ! Evaluate EP-flux from Eq. 7 and give the right orientation to
473      ! the stress
474
475      DO JW = 1, NW
476        RUWP(JW, :) = SIGN(1., ZOP(JW, :)) * COS(ZP(JW, :)) * WWP(JW, :)
477        RVWP(JW, :) = SIGN(1., ZOP(JW, :)) * SIN(ZP(JW, :)) * WWP(JW, :)
478      end DO
479
480      RUW(:, LL + 1) = 0.
481      RVW(:, LL + 1) = 0.
482
483      DO JW = 1, NW
484        RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
485        RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
486        EAST_GWSTRESS(:, LL) = EAST_GWSTRESS(:, LL) + MAX(0., RUWP(JW, :)) / FLOAT(NW)
487        WEST_GWSTRESS(:, LL) = WEST_GWSTRESS(:, LL) + MIN(0., RUWP(JW, :)) / FLOAT(NW)
488      end DO
489    end DO
490
491    ! 5 CALCUL DES TENDANCES:
492
493    ! 5.1 Rectification des flux au sommet et dans les basses couches
494
495    RUW(:, KLEV + 1) = 0.
496    RVW(:, KLEV + 1) = 0.
497    RUW(:, 1) = RUW(:, LAUNCH)
498    RVW(:, 1) = RVW(:, LAUNCH)
499    DO LL = 1, LAUNCH
500      RUW(:, LL) = RUW(:, LAUNCH + 1)
501      RVW(:, LL) = RVW(:, LAUNCH + 1)
502      EAST_GWSTRESS(:, LL) = EAST_GWSTRESS(:, LAUNCH)
503      WEST_GWSTRESS(:, LL) = WEST_GWSTRESS(:, LAUNCH)
504    end DO
505
506    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
507    DO LL = 1, KLEV
508      D_U(:, LL) = (1. - DTIME / DELTAT) * D_U(:, LL) + DTIME / DELTAT / REAL(NW) * &
509              RG * (RUW(:, LL + 1) - RUW(:, LL)) &
510              / (PH(:, LL + 1) - PH(:, LL)) * DTIME
511      !  NO AR1 FOR MERIDIONAL TENDENCIES
512      !      D_V(:, LL) = (1.-DTIME/DELTAT) * D_V(:, LL) + DTIME/DELTAT/REAL(NW) * &
513      D_V(:, LL) = 1. / REAL(NW) * &
514              RG * (RVW(:, LL + 1) - RVW(:, LL)) &
515              / (PH(:, LL + 1) - PH(:, LL)) * DTIME
516    ENDDO
517
518    ! Cosmetic: evaluation of the cumulated stress
519    ZUSTR = 0.
520    ZVSTR = 0.
521    DO LL = 1, KLEV
522      ZUSTR = ZUSTR + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL)) / DTIME
523      !      ZVSTR = ZVSTR + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
524    ENDDO
525    ! COSMETICS TO VISUALIZE ROTBA
526    ZVSTR = ROTBA
527
528  END SUBROUTINE ACAMA_GWD_RANDO
529
530END MODULE ACAMA_GWD_rando_m
Note: See TracBrowser for help on using the repository browser.