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

Last change on this file since 2325 was 2266, checked in by emillour, 5 years ago

Mars GCM:
Some very minor fixes in the physics to be able to compile physics as a library:

  • remove reference to unused .h files located in the dynamics in nonoro_gwd_ran_mod.F
  • fix an OpenMP declaration in time_phylmdz_mod.F90

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