source: LMDZ5/branches/IPSLCM6.0.8/libf/phylmd/acama_gwd_rando_m.F90 @ 4947

Last change on this file since 4947 was 2720, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2664:2719 into testing branch

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