source: LMDZ5/trunk/libf/phylmd/acama_gwd_rando_m.F90 @ 4487

Last change on this file since 4487 was 2665, checked in by dcugnet, 8 years ago
  • A (re)startphy.nc file (standard name: "startphy0.nc") can be read by ce0l to get land mask, so mask can be defined (in decreasing priority order) from: 1) "o2a.nc file" if this file is found 2) "startphy0.nc" if this file is found 3) "Relief.nc" otherwise
  • Sub-cell scales parameters for orographic gravity waves can be read from file "oro_params.nc" if the configuration key "read_orop" is TRUE. The effect is to bypass the "grid_noro" routine in ce0l, so that any pre-defined mask (from o2a.nc or startphy0.nc) is then overwritten.
  • The gcm stops if the "limit.nc" records number differs from the current year number of days. A warning is issued in case the gcm calendar does not match the time axis attribute "calendar" (if available) from the "limit.nc" file. This attribute is now added to the "limit.nc" time axis.
  • Few simplifications in grid_noro
  • Few parameters changes in acama_gwd and flott_gwd.
  • Variable d_u can be saved in the outputs.
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.