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

Last change on this file since 2578 was 2578, checked in by romain.vande, 3 years ago

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