source: LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90 @ 5421

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

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