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

Last change on this file since 2800 was 2682, checked in by emillour, 3 years ago

Mars GCM:
Update Non-Orographic GW scheme: cmax is not longer fixed but evaluated
as the zonal wind at launch altitude.
JL

File size: 25.6 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, cpnew, rnew, 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    ! UPDATED              J.LIU       12/2021  The rho (density) at the specific
38    !                                           locations is introduced. The equation
39    !                                           of EP-flux is corrected by adding the
40    !                                           term of density at launch (source)
41    !                                           altitude.
42    ! 
43    !---------------------------------------------------------------------------------
44
45      use comcstfi_h, only: g, pi, r
46      USE ioipsl_getin_p_mod, ONLY : getin_p
47      use vertical_layers_mod, only : presnivs
48      use geometry_mod, only: cell_area
49#ifdef CPP_XIOS
50     use xios_output_mod, only: send_xios_field
51#endif     
52     
53      implicit none
54      include "callkeys.h"
55
56      CHARACTER (LEN=20) :: modname='nonoro_gwd_ran'
57      CHARACTER (LEN=80) :: abort_message
58
59
60    ! 0. DECLARATIONS:
61
62    ! 0.1 INPUTS
63    INTEGER, intent(in):: ngrid          ! number of atmospheric columns
64    INTEGER, intent(in):: nlayer         ! number of atmospheric layers
65    REAL, intent(in):: DTIME             ! Time step of the Physics(s)
66    REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m)
67    REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere
68    REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere
69    REAL, intent(in):: pp(ngrid,nlayer)  ! Pressure at full levels(Pa)
70    REAL, intent(in):: pt(ngrid,nlayer)  ! Temp at full levels(K)
71    REAL, intent(in):: pu(ngrid,nlayer)  ! Zonal wind at full levels(m/s)
72    REAL, intent(in):: pv(ngrid,nlayer)  ! Meridional winds at full levels(m/s)
73    REAL,INTENT(in) :: pdt(ngrid,nlayer) ! Tendency on temperature (K/s)
74    REAL,INTENT(in) :: pdu(ngrid,nlayer) ! Tendency on zonal wind (m/s/s)
75    REAL,INTENT(in) :: pdv(ngrid,nlayer) ! Tendency on meridional wind (m/s/s)
76
77    ! 0.2 OUTPUTS
78    REAL, intent(out):: zustr(ngrid)       ! Zonal surface stress
79    REAL, intent(out):: zvstr(ngrid)       ! Meridional surface stress
80    REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature (K/s) due to gravity waves (not used set to zero)
81    REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind (m/s/s) due to gravity waves
82    REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves
83
84    ! 0.3 INTERNAL ARRAYS
85    REAL :: TT(ngrid, nlayer)   ! Temperature at full levels
86    REAL :: RHO(ngrid, nlayer)  ! Mass density at full levels
87    REAL :: UU(ngrid, nlayer)   ! Zonal wind at full levels
88    REAL :: VV(ngrid, nlayer)   ! Meridional winds at full levels
89    REAL :: BVLOW(ngrid)        ! initial N at each grid (not used)
90
91    INTEGER II, JJ, LL
92
93    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
94    REAL, parameter:: DELTAT = 24. * 3600.
95   
96    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
97    INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers
98    INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed
99    INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed
100    INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed
101    INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves
102    INTEGER JK, JP, JO, JW      ! Loop index for NK,NP,NO, and NW
103
104    REAL, save :: kmax             ! Max horizontal wavenumber (lambda min,lambda=2*pi/kmax),kmax=N/u, u(=30~50) zonal wind velocity at launch altitude
105!$OMP THREADPRIVATE(kmax)
106    REAL, save :: kmin             ! Min horizontal wavenumber (lambda max = 314 km,lambda=2*pi/kmin)
107!$OMP THREADPRIVATE(kmin)
108    REAL kstar                     ! Min horizontal wavenumber constrained by horizontal resolution
109
110    REAL :: max_k(ngrid)           ! max_k=max(kstar,kmin)
111    REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity (not used)   
112    REAL CPHA(ngrid)                      ! absolute PHASE VELOCITY frequency
113    REAL ZK(NW, ngrid)             ! Horizontal wavenumber amplitude
114    REAL ZP(NW, ngrid)             ! Horizontal wavenumber angle       
115    REAL ZO(NW, ngrid)             ! Absolute frequency
116
117    REAL intr_freq_m(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
118    REAL intr_freq_p(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
119    REAL wwm(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level below the full level
120    REAL wwp(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level above the full level
121    REAL u_epflux_p(nw, ngrid)           ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
122    REAL v_epflux_p(nw, ngrid)           ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
123    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) 
124    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)
125    REAL epflux_0(nw, ngrid)             ! Fluxes at launching level (previous name: RUW0)
126    REAL, save :: epflux_max             ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable)
127!$OMP THREADPRIVATE(epflux_max)
128    INTEGER LAUNCH                       ! Launching altitude
129    REAL, save :: xlaunch                ! Control the launching altitude by pressure
130!$OMP THREADPRIVATE(xlaunch)
131    REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx. not used)
132    REAL cmax(ngrid,nlayer)             ! Max horizontal absolute phase velocity (at the maginitide of zonal wind u at the launch altitude)
133
134    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
135    REAL, save :: sat                 ! saturation parameter(tunable)
136!$OMP THREADPRIVATE(sat)
137    REAL, save :: rdiss               ! dissipation coefficient (tunable)
138!$OMP THREADPRIVATE(rdiss)
139    REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic frequency
140
141    ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
142    REAL H0bis(ngrid, nlayer)          ! Characteristic Height of the atmosphere (specific locations)
143    REAL, save:: H0                    ! Characteristic Height of the atmosphere (constant)
144!$OMP THREADPRIVATE(H0)
145    REAL, parameter:: pr = 250         ! Reference pressure [Pa]
146    REAL, parameter:: tr = 220.        ! Reference temperature [K]
147    REAL ZH(ngrid, nlayer + 1)         ! Log-pressure altitude (constant H0)
148    REAL ZHbis(ngrid, nlayer + 1)      ! Log-pressure altitude (varying H0bis)
149    REAL UH(ngrid, nlayer + 1)         ! zonal wind at 1/2 levels
150    REAL VH(ngrid, nlayer + 1)         ! meridional wind at 1/2 levels
151    REAL PH(ngrid, nlayer + 1)         ! Pressure at 1/2 levels
152    REAL, parameter:: psec = 1.e-20    ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure)
153    REAL BV(ngrid, nlayer + 1)         ! Brunt Vaisala freq. (N^2) at 1/2 levels
154    REAL, parameter:: bvsec = 1.e-5    ! Security to avoid negative BV 
155    REAL HREF(nlayer + 1)              ! Reference pressure for launching altitude
156
157
158    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
159
160
161    LOGICAL,SAVE :: firstcall = .true.
162!$OMP THREADPRIVATE(firstcall)
163
164!-----------------------------------------------------------------------------------------------------------------------
165!  1. INITIALISATIONS
166!-----------------------------------------------------------------------------------------------------------------------
167     IF (firstcall) THEN
168        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
169        epflux_max = 7.E-7 ! Mars' value !!
170        call getin_p("nonoro_gwd_epflux_max", epflux_max)
171        write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max
172        sat = 1. ! default gravity waves saturation value !!
173        call getin_p("nonoro_gwd_sat", sat)
174        write(*,*) "nonoro_gwd_ran: sat=", sat     
175!        cmax = 50. ! default gravity waves phase velocity value !!
176!        call getin_p("nonoro_gwd_cmax", cmax)
177!        write(*,*) "nonoro_gwd_ran: cmax=", cmax
178        rdiss=0.01 ! default coefficient of dissipation !!
179        call getin_p("nonoro_gwd_rdiss", rdiss)
180        write(*,*) "nonoro_gwd_ran: rdiss=", rdiss
181        kmax=7.E-4 ! default Max horizontal wavenumber !!
182        call getin_p("nonoro_gwd_kmax", kmax)
183        write(*,*) "nonoro_gwd_ran: kmax=", kmax
184        kmin=2.e-5 ! default Max horizontal wavenumber !!
185        call getin_p("nonoro_gwd_kmin", kmin)
186        write(*,*) "nonoro_gwd_ran: kmin=", kmin
187        xlaunch=0.4 ! default GW luanch altitude controller !!
188        call getin_p("nonoro_gwd_xlaunch", xlaunch)
189        write(*,*) "nonoro_gwd_ran: xlaunch=", xlaunch
190        ! Characteristic vertical scale height
191        H0 = r * tr / g
192        ! Control
193        if (deltat .LT. dtime) THEN
194             call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
195        endif
196        if (nlayer .LT. nw) THEN
197             call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
198        endif
199        firstcall = .false.
200     ENDIF
201
202! Compute current values of temperature and winds
203    tt(:,:)=pt(:,:)+dtime*pdt(:,:)
204    uu(:,:)=pu(:,:)+dtime*pdu(:,:)
205    vv(:,:)=pv(:,:)+dtime*pdv(:,:)
206! Compute the real mass density by rho=p/R(T)T
207     DO ll=1,nlayer
208       DO ii=1,ngrid
209            rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll))
210       ENDDO
211     ENDDO
212!    print*,'epflux_max just after firstcall:', epflux_max
213
214!-----------------------------------------------------------------------------------------------------------------------
215!  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
216!-----------------------------------------------------------------------------------------------------------------------
217    ! Pressure and inverse of pressure at 1/2 level
218    DO LL = 2, nlayer
219       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
220    end DO
221    PH(:, nlayer + 1) = 0.
222    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
223
224    call writediagfi(ngrid,'nonoro_pp','nonoro_pp', 'm',3,PP(:,1:nlayer))
225    call writediagfi(ngrid,'nonoro_ph','nonoro_ph', 'm',3,PH(:,1:nlayer))
226
227    ! Launching level for reproductible case
228    !Pour revenir a la version non reproductible en changeant le nombre de
229    !process
230    ! Reprend la formule qui calcule PH en fonction de PP=play
231    DO LL = 2, nlayer
232       HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
233    end DO
234    HREF(nlayer + 1) = 0.
235    HREF(1) = 2. * presnivs(1) - HREF(2)
236
237    LAUNCH=0
238    DO LL = 1, nlayer
239       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
240    ENDDO
241       
242    ! Log pressure vert. coordinate
243       ZH(:,1) = 0.
244    DO LL = 2, nlayer + 1
245       !ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
246        H0bis(:, LL-1) = r * TT(:, LL-1) / g
247          ZH(:, LL) = ZH(:, LL-1) &
248           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
249    end DO
250        ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC))
251
252    call writediagfi(ngrid,'nonoro_zh','nonoro_zh', 'm',3,ZH(:,2:nlayer+1))
253
254    ! Winds and Brunt Vaisala frequency
255    DO LL = 2, nlayer
256       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1))                          ! Zonal wind
257       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1))                          ! Meridional wind
258       ! Brunt Vaisala frequency (=g/T*[dT/dz + g/cp] )
259       BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1)))                     &
260        *((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+ &
261        G / cpnew(:,LL))
262
263       BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive
264       BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small
265    end DO
266    BV(:, 1) = BV(:, 2)
267    UH(:, 1) = 0.
268    VH(:, 1) = 0.
269    BV(:, nlayer + 1) = BV(:, nlayer)
270    UH(:, nlayer + 1) = UU(:, nlayer)
271    VH(:, nlayer + 1) = VV(:, nlayer)
272    cmax(:,launch)=UU(:, launch)
273    DO ii=1,ngrid
274       KSTAR = PI/SQRT(cell_area(II))
275       MAX_K(II)=MAX(kmin,kstar)
276    ENDDO
277    call writediagfi(ngrid,'nonoro_bv','nonoro_bv', 'm',3,BV(:,2:nlayer+1))
278
279!-----------------------------------------------------------------------------------------------------------------------
280! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
281!-----------------------------------------------------------------------------------------------------------------------
282! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way
283    DO JW = 1, NW
284             !  Angle
285             DO II = 1, ngrid
286                ! Angle (0 or PI so far)
287                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
288                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
289                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.)* PI / 2.
290               
291                ! Horizontal wavenumber amplitude
292                ! From Venus model: TN+GG April/June 2020 (rev2381)
293                ! "Individual waves are not supposed to occupy the entire domain, but only a fraction of it" (Lott+2012)
294                ! ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
295                KSTAR = PI/SQRT(cell_area(II))         
296                ZK(JW, II) = MAX_K(II) + (KMAX - MAX_K(II)) *RAN_NUM_2
297               
298               ! Horizontal phase speed
299! this computation allows to get a gaussian distribution centered on 0 (right ?)
300! then cmin is not useful, and it favors small values.
301                CPHA(:) = 0.
302                DO JJ = 1, NA
303                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
304                    CPHA(ii) = CPHA(ii) + 2.*CMAX(ii,launch)*      &
305                           (RAN_NUM_3 -0.5)*                       &
306                           SQRT(3.)/SQRT(NA*1.)
307                END DO
308                IF (CPHA(ii).LT.0.)  THEN
309                   CPHA(ii) = -1.*CPHA(ii)
310                   ZP(JW,II) = ZP(JW,II) + PI
311                ENDIF
312! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
313!           ran_num_3 = mod(tt(ii, jw)**2, 1.)
314!           cpha = cmin + (cmax - cmin) * ran_num_3
315
316                ! Intrinsic frequency
317                ZO(JW, II) = CPHA(II) * ZK(JW, II)
318                ! Intrinsic frequency  is imposed
319                ZO(JW, II) = ZO(JW, II)                                            &
320                            + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH)        &
321                            + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
322               
323                ! Momentum flux at launch level
324                epflux_0(JW, II) = epflux_max                                      &
325                                  * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.)
326              ENDDO
327   end DO
328!     print*,'epflux_0 just after waves charac. ramdon:', epflux_0
329
330!-----------------------------------------------------------------------------------------------------------------------
331! 4. COMPUTE THE FLUXES
332!-----------------------------------------------------------------------------------------------------------------------
333    !  4.1  Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes.
334    !------------------------------------------------------
335    DO JW = 1, NW
336       ! Evaluate intrinsic frequency at launching altitude:
337       intr_freq_p(JW, :) = ZO(JW, :)                                    &
338                            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
339                            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
340    end DO
341
342! VERSION WITHOUT CONVECTIVE SOURCE
343       ! Vertical velocity at launch level, value to ensure the imposed
344       ! mom flux:
345         DO JW = 1, NW
346       ! WW is directly a flux, here, not vertical velocity anymore
347            WWP(JW, :) = epflux_0(JW,:)
348            u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
349            v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
350         end DO
351   
352    !  4.2 Initial flux at launching altitude
353    !------------------------------------------------------
354    u_epflux_tot(:, LAUNCH) = 0
355    v_epflux_tot(:, LAUNCH) = 0
356    DO JW = 1, NW
357       u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :)
358       v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :)
359    end DO
360
361    !  4.3 Loop over altitudes, with passage from one level to the next done by:
362    !----------------------------------------------------------------------------
363    !    i) conserving the EP flux,
364    !    ii) dissipating a little,
365    !    iii) testing critical levels,
366    !    iv) testing the breaking.
367    !----------------------------------------------------------------------------
368    DO LL = LAUNCH, nlayer - 1
369       !  W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL TO THE NEXT)
370       DO JW = 1, NW
371          intr_freq_m(JW, :) = intr_freq_p(JW, :)
372          WWM(JW, :) = WWP(JW, :)
373          ! Intrinsic Frequency
374          intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1)     &
375                               - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1)
376          ! Vertical velocity in flux formulation
377          WWP(JW, :) = MIN(                                                              &
378                     ! No breaking (Eq.6)
379                     WWM(JW, :)                                                          &
380                     ! Dissipation (Eq. 8)(real rho used here rather than pressure
381                     ! parametration because the original code has a bug if the density of
382                     ! the planet at the launch altitude not approximate 1):                     !
383                     * EXP(- RDISS*2./rho(:, LL + 1)                                     &
384                     * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3                             &
385                     / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 &
386                     * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))),                      &
387                     ! Critical levels (forced to zero if intrinsic frequency
388                     ! changes sign)
389                     MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :)))          &
390                     ! Saturation (Eq. 12) (rho at launch altitude is imposed by J.Liu.
391                     ! Same reason with Eq. 8)
392                     * ABS(intr_freq_p(JW, :))**3 / (2.*BV(:, LL+1))                     &
393                     * rho(:,launch)*exp(-zh(:, ll + 1) / H0)                            &
394                     * SAT**2 *KMIN**2 / ZK(JW, :)**4) 
395       end DO
396       ! END OF W(KB)ARNING
397
398       ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress
399       DO JW = 1, NW
400          u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) *  WWP(JW, :)
401          v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) *  WWP(JW, :)
402       end DO
403 
404       u_epflux_tot(:, LL + 1) = 0.
405       v_epflux_tot(:, LL + 1) = 0.
406       DO JW = 1, NW
407          u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :)
408          v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :)
409          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW)
410          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW)
411       end DO       
412    end DO ! DO LL = LAUNCH, nlayer - 1
413
414!-----------------------------------------------------------------------------------------------------------------------
415! 5. TENDENCY CALCULATIONS
416!-----------------------------------------------------------------------------------------------------------------------
417
418    ! 5.1 Flow rectification at the top and in the low layers
419    ! --------------------------------------------------------
420    ! Warning, this is the total on all GW
421    u_epflux_tot(:, nlayer + 1) = 0.
422    v_epflux_tot(:, nlayer + 1) = 0.
423
424    ! Here, big change compared to FLott version:
425    ! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
426    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
427    DO LL = 1, max(1,LAUNCH-3)
428      u_epflux_tot(:, LL) = 0.
429      v_epflux_tot(:, LL) = 0.
430    end DO
431    DO LL = max(2,LAUNCH-2), LAUNCH-1
432       u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) *          &
433                             (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
434       v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) *          &
435                             (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
436       EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) +                                   &
437                             EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/            &
438                             (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
439       WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) +                                   &
440                             WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/            &
441                             (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
442    end DO
443    ! This way, the total flux from GW is zero, but there is a net transport
444    ! (upward) that should be compensated by circulation
445    ! and induce additional friction at the surface
446    call writediagfi(ngrid,'nonoro_u_epflux_tot','nonoro_u_epflux_tot', '',3,u_epflux_tot(:,2:nlayer+1))
447    call writediagfi(ngrid,'nonoro_v_epflux_tot','nonoro_v_epflux_tot', '',3,v_epflux_tot(:,2:nlayer+1))
448
449    ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
450    !---------------------------------------------
451    DO LL = 1, nlayer
452       !Notice here the d_u and d_v are tendency (For Mars) but not increment(Venus).
453       d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
454                    / (PH(:, LL + 1) - PH(:, LL))
455       d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
456                    / (PH(:, LL + 1) - PH(:, LL))
457    ENDDO
458    d_t(:,:) = 0.
459    call writediagfi(ngrid,'nonoro_d_u','nonoro_d_u', '',3,d_u)
460    call writediagfi(ngrid,'nonoro_d_v','nonoro_d_v', '',3,d_v)
461
462    ! 5.3 Update tendency of wind with the previous (and saved) values
463    !-----------------------------------------------------------------
464    d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:)                       &
465               + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
466    d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:)                       &
467               + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
468    du_nonoro_gwd(:,:) = d_u(:,:)
469    dv_nonoro_gwd(:,:) = d_v(:,:)
470
471    call writediagfi(ngrid,'du_nonoro_gwd','du_nonoro_gwd', '',3,du_nonoro_gwd)
472    call writediagfi(ngrid,'dv_nonoro_gwd','dv_nonoro_gwd', '',3,dv_nonoro_gwd)
473   
474    ! Cosmetic: evaluation of the cumulated stress
475    ZUSTR(:) = 0.
476    ZVSTR(:) = 0.
477    DO LL = 1, nlayer
478       ZUSTR(:) = ZUSTR(:) + D_U(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME
479       ZVSTR(:) = ZVSTR(:) + D_V(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME
480    ENDDO
481
482!#ifdef CPP_XIOS
483!     call send_xios_field("du_nonoro", d_u)
484!     call send_xios_field("dv_nonoro", d_v)
485!     call send_xios_field("east_gwstress", east_gwstress)
486!     call send_xios_field("west_gwstress", west_gwstress)
487!#endif
488
489  END SUBROUTINE NONORO_GWD_RAN
490
491
492
493! ========================================================
494! Subroutines used to allocate/deallocate module variables       
495! ========================================================
496  SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
497
498  IMPLICIT NONE
499
500      INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
501      INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
502
503         allocate(du_nonoro_gwd(ngrid,nlayer))
504         allocate(dv_nonoro_gwd(ngrid,nlayer))
505         allocate(east_gwstress(ngrid,nlayer))
506         east_gwstress(:,:)=0
507         allocate(west_gwstress(ngrid,nlayer))
508         west_gwstress(:,:)=0
509
510  END SUBROUTINE ini_nonoro_gwd_ran
511! ----------------------------------
512  SUBROUTINE end_nonoro_gwd_ran
513
514  IMPLICIT NONE
515
516         if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
517         if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
518         if (allocated(east_gwstress)) deallocate(east_gwstress)
519         if (allocated(west_gwstress)) deallocate(west_gwstress)
520
521  END SUBROUTINE end_nonoro_gwd_ran
522
523END MODULE nonoro_gwd_ran_mod
Note: See TracBrowser for help on using the repository browser.