source: trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90 @ 2538

Last change on this file since 2538 was 2400, checked in by dbardet, 5 years ago

Add kstar parameter to control kmin value. kmi=1/lambda_max, to ensure the subgrid scale characteristic, we have to constrain the GWs wavelength by the size of the mesh
File size: 21.8 KB
Line 
1MODULE nonoro_gwd_ran_mod
2
3IMPLICIT NONE
4
5REAL,allocatable,save :: du_nonoro_gwd(:,:) ! Zonal wind tendency due to GWD
6REAL,allocatable,save :: dv_nonoro_gwd(:,:) ! Meridional wind tendency due to GWD
7REAL,ALLOCATABLE,SAVE :: east_gwstress(:,:) ! Profile of eastward stress
8REAL,ALLOCATABLE,SAVE :: west_gwstress(:,:) ! Profile of westward stress
9
10CONTAINS
11
12      SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, pp,  &
13                  zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
14                  zustr,zvstr,d_t, d_u, d_v)
15
16    !--------------------------------------------------------------------------------
17    ! Parametrization of the momentum flux deposition due to a discrete
18    ! number of gravity waves.
19    ! F. Lott
20    ! Version 14, Gaussian distribution of the source
21    ! LMDz model online version     
22    ! ADAPTED FOR VENUS /  F. LOTT + S. LEBONNOIS
23    ! Version adapted on 03/04/2013:
24    !      - input flux compensated in the deepest layers
25    !                           
26    ! ADAPTED FOR MARS     G.GILLI     02/2016
27    !        Revision with F.Forget    06/2016  Variable EP-flux according to
28    !                                           PBL variation (max velocity thermals)
29    ! UPDATED              D.BARDET    01/2020  - reproductibility of the
30    !                                           launching altitude calculation
31    !                                           - wave characteristic
32    !                                           calculation using MOD
33    !                                           - adding east_gwstress and
34    !                                           west_gwstress variables   
35    !---------------------------------------------------------------------------------
36
37      use comcstfi_h, only: g, pi, cpp, r
38      USE ioipsl_getin_p_mod, ONLY : getin_p
39      use assert_m, only : assert
40      use vertical_layers_mod, only : presnivs
41      use geometry_mod, only: cell_area
42      implicit none
43
44      include "yoegwd.h"
45      include "callkeys.h"
46
47      CHARACTER (LEN=20) :: modname='flott_gwd_rando'
48      CHARACTER (LEN=80) :: abort_message
49
50
51    ! 0. DECLARATIONS:
52
53    ! 0.1 INPUTS
54    INTEGER, intent(in):: ngrid, nlayer
55    REAL, intent(in):: DTIME ! Time step of the Physics
56    REAL, intent(in):: zmax_therm(ngrid) ! altitude of max velocity thermals (m)
57    REAL, intent(in):: pp(ngrid,nlayer)   ! Pressure at full levels
58    REAL, intent(in):: pt(ngrid,nlayer)   ! Temp at full levels
59    REAL, intent(in):: pu(ngrid,nlayer),pv(ngrid,nlayer) ! Hor winds at full levels
60    REAL,INTENT(in) :: pdt(ngrid,nlayer) ! tendency on temperature
61    REAL,INTENT(in) :: pdu(ngrid,nlayer) ! tendency on zonal wind
62    REAL,INTENT(in) :: pdv(ngrid,nlayer) ! tendency on meridional wind
63
64    ! 0.2 OUTPUTS
65    REAL, intent(out):: zustr(ngrid), zvstr(ngrid) ! Surface Stresses
66    REAL, intent(out):: d_t(ngrid, nlayer)        ! Tendency on Temp.
67    REAL, intent(out):: d_u(ngrid, nlayer), d_v(ngrid, nlayer) ! tendency on winds
68
69    ! O.3 INTERNAL ARRAYS
70    REAL :: TT(ngrid, nlayer)   ! Temp at full levels
71    REAL :: UU(ngrid, nlayer) , VV(ngrid, nlayer) ! Hor winds at full levels
72    REAL :: BVLOW(ngrid)
73    REAL :: DZ
74
75    INTEGER II, JJ, LL
76
77    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
78
79    REAL, parameter:: DELTAT = 24. * 3600.
80   
81    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
82
83    INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers
84    INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed
85    INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed
86    INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves
87    INTEGER JK, JP, JO, JW
88    INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed
89    REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber
90    REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber
91    REAL kstar                     ! Min horizontal wavenumber constrained by horizontal resolution
92    REAL, parameter:: cmax = 30.   ! Max horizontal absolute phase velocity
93    REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity   
94    REAL CPHA                    ! absolute PHASE VELOCITY frequency
95    REAL ZK(NW, ngrid)           ! Horizontal wavenumber amplitude
96    REAL ZP(NW, ngrid)           ! Horizontal wavenumber angle       
97    REAL ZO(NW, ngrid)           ! Absolute frequency
98
99    REAL intr_freq_m(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
100    REAL intr_freq_p(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
101    REAL wwm(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level below the full level
102    REAL wwp(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level above the full level
103    REAL u_epflux_p(nw, ngrid)           ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
104    REAL v_epflux_p(nw, ngrid)           ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
105    REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW) 
106    REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW)
107    REAL epflux_0(nw, ngrid)             ! Fluxes at launching level (previous name: RUW0)
108    REAL, save :: epflux_max             ! Max EP flux value at launching altitude (previous name: RUWMAX)
109    INTEGER LAUNCH      ! Launching altitude
110    REAL, parameter:: xlaunch = 0.4      ! Control the launching altitude
111    REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.)
112
113
114    REAL PREC(ngrid)
115    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
116
117
118    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
119    REAL, parameter:: sat   = 1.     ! saturation parameter
120    REAL, parameter:: rdiss = 1.     ! coefficient of dissipation
121    REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency
122
123    ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
124    REAL H0bis(ngrid, nlayer)          ! Characteristic Height of the atmosphere
125    REAL, save:: H0                 ! Characteristic Height of the atmosphere
126    REAL, parameter:: pr = 250      ! Reference pressure [Pa]
127    REAL, parameter:: tr = 220.     ! Reference temperature [K]
128    REAL ZH(ngrid, nlayer + 1)         ! Log-pressure altitude (constant H0)
129    REAL ZHbis(ngrid, nlayer + 1)      ! Log-pressure altitude (varying H)
130    REAL UH(ngrid, nlayer + 1)         ! zonal wind at 1/2 levels
131    REAL VH(ngrid, nlayer + 1)         ! meridional wind at 1/2 levels
132    REAL PH(ngrid, nlayer + 1)         ! Pressure at 1/2 levels
133    REAL, parameter:: psec = 1.e-6  ! Security to avoid division by 0 pressure
134    REAL BV(ngrid, nlayer + 1)         ! Brunt Vaisala freq. (BVF) at 1/2 levels
135    REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV 
136    REAL HREF(nlayer + 1)             ! Reference altitude for launching alt.
137
138
139! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
140    logical,save :: output=.false.
141 ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
142    character*14 outform
143    character*2  str2
144    integer      ieq
145
146    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
147
148
149    LOGICAL,SAVE :: firstcall = .true.
150
151
152   !-----------------------------------------------------------------
153    !  1. INITIALISATIONS
154
155     IF (firstcall) THEN
156        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
157        epflux_max = 7.E-7 ! Mars' value !!
158        call getin_p("nonoro_gwd_epflux_max", epflux_max)
159        write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max
160        ! Characteristic vertical scale height
161        H0 = r * tr / g
162        ! Control
163        if (deltat .LT. dtime) THEN
164             call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
165        endif
166        if (nlayer .LT. nw) THEN
167             call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
168        endif
169        firstcall = .false.
170     ENDIF
171
172    gwd_convective_source=.false.
173
174    ! Compute current values of temperature and winds
175    tt(:,:)=pt(:,:)+dtime*pdt(:,:)
176    uu(:,:)=pu(:,:)+dtime*pdu(:,:)
177    vv(:,:)=pv(:,:)+dtime*pdv(:,:)
178
179
180
181    !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
182    !-------------------------------------------------------------
183
184    !Online output
185    if (output) OPEN(11,file="impact-gwno.dat")
186
187    ! Pressure and Inv of pressure, Temperature / at 1/2 level
188    DO LL = 2, nlayer
189       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
190    end DO
191
192    PH(:, nlayer + 1) = 0.
193    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
194
195    ! Launching altitude
196
197    !Pour revenir a la version non reproductible en changeant le nombre de
198    !process
199    ! Reprend la formule qui calcule PH en fonction de PP=play
200    DO LL = 2, nlayer
201       HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
202    end DO
203    HREF(nlayer + 1) = 0.
204    HREF(1) = 2. * presnivs(1) - HREF(2)
205
206    LAUNCH=0
207    DO LL = 1, nlayer
208       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
209    ENDDO
210 
211    if (output) print*, " WE ARE IN FLOTT GW SCHEME "
212   
213    ! Log pressure vert. coordinate
214    DO LL = 1, nlayer + 1
215       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
216    end DO
217
218    if (output) then
219    ! altitude above surface
220       ZHbis(:,1) = 0.   
221       DO LL = 2, nlayer + 1
222          H0bis(:, LL-1) = r * TT(:, LL-1) / g
223          ZHbis(:, LL) = ZHbis(:, LL-1) &
224           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
225       end DO
226    endif
227
228    ! Winds and BV frequency
229    DO LL = 2, nlayer
230       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
231       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
232       ! GG test       
233       !print*, 'TT, UH, VH, ZH at launch', TT(ngrid/2,LAUNCH), UH(ngrid/2,LAUNCH),VH(ngrid/2, LAUNCH), ZH(ngrid/2,LAUNCH)     
234       ! BVSEC: BV Frequency
235       BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
236            * r**2 / cpp / H0**2 + (TT(:, LL) &
237            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * r / H0
238       BV(:,LL) =SQRT(MAX(BVSEC,BV(:,LL)))
239    end DO
240       !GG test
241       !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(ngrid/2, LAUNCH) 
242
243    BV(:, 1) = BV(:, 2)
244    UH(:, 1) = 0.
245    VH(:, 1) = 0.
246    BV(:, nlayer + 1) = BV(:, nlayer)
247    UH(:, nlayer + 1) = UU(:, nlayer)
248    VH(:, nlayer + 1) = VV(:, nlayer)
249
250
251    ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
252    !-------------------------------------------
253
254    ! The mod function of here a weird arguments
255    ! are used to produce the waves characteristics
256    ! in a stochastic way
257
258    DO JW = 1, NW
259             !  Angle
260             DO II = 1, ngrid
261                ! Angle (0 or PI so far)
262                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
263                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
264                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
265                     * PI / 2.
266                ! Horizontal wavenumber amplitude
267                ! From Venus model: TN+GG April/June 2020 (rev2381)
268                ! "Individual waves are not supposed to occupy the
269                ! entire domain, but only a fraction of it" (Lott+2012)
270               
271                !ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
272                KSTAR = PI/SQRT(cell_area(II))
273                ZK(JW, II) = MAX(KMIN,KSTAR) + (KMAX - MAX(KMIN,KSTAR)) *RAN_NUM_2
274                ! Horizontal phase speed
275                CPHA = 0.
276                DO JJ = 1, NA
277                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
278                    CPHA = CPHA + &
279                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
280                END DO
281                IF (CPHA.LT.0.)  THEN
282                   CPHA = -1.*CPHA
283                   ZP(JW,II) = ZP(JW,II) + PI
284                ENDIF
285       !Online output: linear
286                if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1)
287                ! Intrinsic frequency
288                ZO(JW, II) = CPHA * ZK(JW, II)
289                ! Intrinsic frequency  is imposed
290                    ZO(JW, II) = ZO(JW, II)      &
291                  + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
292                  + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
293                ! Momentum flux at launch lev
294                ! epflux_0(JW, II) = epflux_max / REAL(NW) &
295                epflux_0(JW, II) = epflux_max &
296                     * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.)
297       !Online output: fixed to max
298                if (output) epflux_0(JW, II) = epflux_max
299             ENDDO
300   end DO
301
302    ! 4. COMPUTE THE FLUXES
303    !--------------------------
304
305    !  4.1  Vertical velocity at launching altitude to ensure
306    !       the correct value to the imposed fluxes.
307    !
308    DO JW = 1, NW
309       ! Evaluate intrinsic frequency at launching altitude:
310       intr_freq_p(JW, :) = ZO(JW, :) &
311            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
312            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
313    end DO
314
315    IF (gwd_convective_source) THEN
316         DO JW = 1, NW
317       ! VERSION WITH CONVECTIVE SOURCE (designed for Earth)
318
319       ! Vertical velocity at launch level, value to ensure the
320       ! imposed mmt flux factor related to the convective forcing:
321       ! precipitations.
322
323       ! tanh limitation to values above prmax:
324!       WWP(JW, :) = epflux_0(JW, :) &
325!            * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
326!       Here, we neglected the kinetic energy providing of the thermodynamic
327!       phase change
328
329!
330
331       ! Factor related to the characteristics of the waves:
332            WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:)  &
333                 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**3
334
335      ! Moderation by the depth of the source (dz here):
336            WWP(JW, :) = WWP(JW, :) &
337                 * EXP(- BVLOW(:)**2 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
338                 * ZK(JW, :)**2 * DZ**2)
339
340      ! Put the stress in the right direction:
341            u_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
342                 * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
343            v_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
344                 * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
345         end DO
346    ELSE ! VERSION WITHOUT CONVECTIVE SOURCE
347       ! Vertical velocity at launch level, value to ensure the imposed
348       ! mom flux:
349         DO JW = 1, NW
350       ! WW is directly a flux, here, not vertical velocity anymore
351            WWP(JW, :) = epflux_0(JW,:)
352            u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
353            v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
354
355         end DO
356    ENDIF
357    !  4.2 Initial flux at launching altitude
358
359    u_epflux_tot(:, LAUNCH) = 0
360    v_epflux_tot(:, LAUNCH) = 0
361    DO JW = 1, NW
362       u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :)
363       v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :)
364    end DO
365
366    !  4.3 Loop over altitudes, with passage from one level to the
367    !      next done by i) conserving the EP flux, ii) dissipating
368    !      a little, iii) testing critical levels, and vi) testing
369    !      the breaking.
370
371    !Online output
372    if (output) then
373        ieq=ngrid/2+1
374        write(str2,'(i2)') NW+2
375        outform="("//str2//"(E12.4,1X))"
376        WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
377               (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW, IEQ)), JW = 1, NW)
378    endif
379
380    DO LL = LAUNCH, nlayer - 1
381
382
383       !  W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL
384       ! TO THE NEXT)
385       DO JW = 1, NW
386          intr_freq_m(JW, :) = intr_freq_p(JW, :)
387          WWM(JW, :) = WWP(JW, :)
388          ! Intrinsic Frequency
389          intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) &
390               - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1)
391
392          WWP(JW, :) = MIN( &
393       ! No breaking (Eq.6)
394               WWM(JW, :) &
395      ! Dissipation (Eq. 8):
396               * EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
397               * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
398               / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 &
399               * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), &
400       ! Critical levels (forced to zero if intrinsic frequency changes sign)
401               MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) &
402       ! Saturation (Eq. 12)
403               * ABS(intr_freq_p(JW, :))**3 /BV(:, LL+1) &
404               * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4) 
405       end DO
406
407       ! END OF W(KB)ARNING
408       ! Evaluate EP-flux from Eq. 7 and
409       ! Give the right orientation to the stress
410       DO JW = 1, NW
411          u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) *  WWP(JW, :)
412          v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) *  WWP(JW, :)
413       end DO
414 
415       u_epflux_tot(:, LL + 1) = 0.
416       v_epflux_tot(:, LL + 1) = 0.
417
418       DO JW = 1, NW
419          u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :)
420          v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :)
421          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW)
422          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW)
423       end DO
424       !Online output
425       if (output) then
426         do JW=1,NW
427            if(u_epflux_p(JW, IEQ).gt.0.) then
428              u_epflux_p(JW, IEQ) = max(u_epflux_p(JW, IEQ), 1.e-99)
429            else
430              u_epflux_p(JW, IEQ) = min(u_epflux_p(JW, IEQ), -1.e-99)
431            endif
432         enddo
433                   WRITE(11,outform) ZH(IEQ, LL+1) / 1000., &
434                                  ZHbis(IEQ, LL+1) / 1000., &
435                                  (u_epflux_p(JW, IEQ), JW = 1, NW)
436       endif
437
438    end DO
439
440    ! 5 CALCUL DES TENDANCES:
441    !------------------------
442
443    ! 5.1 Rectification des flux au sommet et dans les basses couches:
444
445! Attention, ici c'est le total sur toutes les ondes...
446
447    u_epflux_tot(:, nlayer + 1) = 0.
448    v_epflux_tot(:, nlayer + 1) = 0.
449
450    ! Here, big change compared to FLott version:
451    ! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
452    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
453    DO LL = 1, max(1,LAUNCH-3)
454      u_epflux_tot(:, LL) = 0.
455      v_epflux_tot(:, LL) = 0.
456    end DO
457    DO LL = max(2,LAUNCH-2), LAUNCH-1
458       u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * &
459            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
460       v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * &
461            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
462       EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + &
463            EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
464            (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
465       WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + &
466            WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
467            (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
468    end DO
469    ! This way, the total flux from GW is zero, but there is a net transport
470    ! (upward) that should be compensated by circulation
471    ! and induce additional friction at the surface
472
473    !Online output
474    if (output) then
475       DO LL = 1, nlayer - 1
476           WRITE(11,*) ZHbis(IEQ, LL)/1000.,u_epflux_tot(IEQ,LL)
477       end DO
478       CLOSE(11)
479       CALL abort_physic("nonoro_gwd_ran","stop here, the work is done",0)
480    endif
481
482
483    ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
484    !---------------------------------------------
485    DO LL = 1, nlayer
486       d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
487            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
488       d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
489            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
490    ENDDO
491    d_t(:,:) = 0.
492
493    ! 5.3 Update tendency of wind with the previous (and saved) values
494    !-----------------------------------------------------------------
495    d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:)                       &
496             + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
497    d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:)                       &
498             + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
499    du_nonoro_gwd(:,:) = d_u(:,:)
500    dv_nonoro_gwd(:,:) = d_v(:,:)
501
502    ! Cosmetic: evaluation of the cumulated stress
503
504    ZUSTR(:) = 0.
505    ZVSTR(:) = 0.
506    DO LL = 1, nlayer
507       ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
508       ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
509    ENDDO
510
511  END SUBROUTINE NONORO_GWD_RAN
512
513! ========================================================
514! Subroutines used to allocate/deallocate module variables       
515! ========================================================
516
517  SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
518
519  IMPLICIT NONE
520
521      INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
522      INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
523
524         allocate(du_nonoro_gwd(ngrid,nlayer))
525         allocate(dv_nonoro_gwd(ngrid,nlayer))
526         allocate(east_gwstress(ngrid,nlayer))
527         east_gwstress(:,:)=0
528         allocate(west_gwstress(ngrid,nlayer))
529         west_gwstress(:,:)=0
530
531  END SUBROUTINE ini_nonoro_gwd_ran
532! ----------------------------------
533  SUBROUTINE end_nonoro_gwd_ran
534
535  IMPLICIT NONE
536
537         if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
538         if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
539         if (allocated(east_gwstress)) deallocate(east_gwstress)
540         if (allocated(west_gwstress)) deallocate(west_gwstress)
541
542  END SUBROUTINE end_nonoro_gwd_ran
543
544END MODULE nonoro_gwd_ran_mod
Note: See TracBrowser for help on using the repository browser.