Ignore:
Timestamp:
Jan 22, 2020, 11:41:33 AM (5 years ago)
Author:
dbardet
Message:

Update the nonorographic gravity waves drag parametrization (from the r3599 of Earth s model LMDZ6): 1) add east_ and west_gwdstress variables, 2) delete aleas function: random waves are producted by the MOD function, 3) reproductibility concerning the number of procs is now validated, 4) changing name of some variables like RUW, RVW, ZOP, ZOM,..., 5) tendency of winds due to GW drag are now module variables and written in restart files.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90

    r2173 r2220  
    33IMPLICIT NONE
    44
     5REAL,allocatable,save :: du_nonoro_gwd(:,:) ! Zonal wind tendency due to GWD
     6REAL,allocatable,save :: dv_nonoro_gwd(:,:) ! Meridional wind tendency due to GWD
     7
    58CONTAINS
    69
    7       SUBROUTINE NONORO_GWD_RAN(NLON,NLEV,DTIME, pp,  &
     10      SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, pp,  &
    811                  zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
    9                   zustr,zvstr,d_t, d_u, d_v)
     12                  zustr,zvstr,d_t, d_u, d_v, &
     13                  east_gwstress, west_gwstress)
    1014
    1115    !--------------------------------------------------------------------------------
    … …  
    2125    ! ADAPTED FOR MARS     G.GILLI     02/2016
    2226    !        Revision with F.Forget    06/2016  Variable EP-flux according to
    23     !                                            PBL variation (max velocity thermals)
     27    !                                           PBL variation (max velocity thermals)
     28    ! UPDATED              D.BARDET    01/2020  - reproductibility of the
     29    !                                           launching altitude calculation
     30    !                                           - wave characteristic
     31    !                                           calculation using MOD
     32    !                                           - adding east_gwstress and
     33    !                                           west_gwstress variables   
    2434    !---------------------------------------------------------------------------------
    2535
    26       use dimphy, only: klon, klev
    27 !      use moyzon_mod, only: plevmoy
    2836      use comcstfi_h, only: g, pi, cpp, r
    29 !      use comgeomfi_h, only: latitude
    3037      USE ioipsl_getin_p_mod, ONLY : getin_p
     38      use assert_m, only : assert
     39      use vertical_layers_mod, only : presnivs
    3140      implicit none
    3241
    … …  
    3443      include "paramet.h"
    3544      include "yoegwd.h"
     45      include "callkeys.h"
     46
     47      CHARACTER (LEN=20) :: modname='flott_gwd_rando'
     48      CHARACTER (LEN=80) :: abort_message
     49
    3650
    3751    ! 0. DECLARATIONS:
    3852
    3953    ! 0.1 INPUTS
    40     INTEGER, intent(in):: NLON, NLEV
     54    INTEGER, intent(in):: ngrid, nlayer
    4155    REAL, intent(in):: DTIME ! Time step of the Physics
    42     REAL, intent(in):: zmax_therm(NLON) ! altitude of max velocity thermals (m)
    43     REAL, intent(in):: pp(NLON,NLEV)   ! Pressure at full levels
    44     REAL, intent(in):: pt(NLON,NLEV)   ! Temp at full levels
    45     REAL, intent(in):: pu(NLON,NLEV),pv(NLON,NLEV) ! Hor winds at full levels
    46     REAL,INTENT(in) :: pdt(nlon,nlev) ! tendency on temperature
    47     REAL,INTENT(in) :: pdu(nlon,nlev) ! tendency on zonal wind
    48     REAL,INTENT(in) :: pdv(nlon,nlev) ! tendency on meridional wind
     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
    4963
    5064    ! 0.2 OUTPUTS
    51     REAL, intent(out):: zustr(NLON), zvstr(NLON) ! Surface Stresses
    52     REAL, intent(inout):: d_t(NLON, NLEV)        ! Tendency on Temp.
    53 
    54 
    55     REAL, intent(inout):: d_u(NLON, NLEV), d_v(NLON, NLEV)
    56     ! Tendencies on winds
     65    REAL, intent(out):: zustr(ngrid), zvstr(ngrid) ! Surface Stresses
     66    REAL, intent(inout):: d_t(ngrid, nlayer)        ! Tendency on Temp.
     67    REAL, intent(inout):: d_u(ngrid, nlayer), d_v(ngrid, nlayer) ! tendency on winds
     68    REAL, intent(inout):: east_gwstress(ngrid, nlayer) ! Profile of eastward stress
     69    REAL, intent(inout):: west_gwstress(ngrid, nlayer) ! Profile of westward stress
    5770
    5871    ! O.3 INTERNAL ARRAYS
    59     REAL :: TT(NLON, NLEV)   ! Temp at full levels
    60 
    61     REAL :: UU(NLON, NLEV) , VV(NLON, NLEV)
    62     ! Hor winds at full levels
    63 
    64     INTEGER II, LL
     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
    6578
    6679    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
    6780
    68     REAL DELTAT
     81    REAL, parameter:: DELTAT = 24. * 3600.
    6982   
    7083    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
    7184
    72 !!  NO= values of C, NW= total of GW
    73     INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
    74     !Online output: change NO
    75 !    INTEGER, PARAMETER:: NK = 1, NP = 2, NO = 11, NW = NK * NP * NO
     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
    7689    INTEGER JK, JP, JO, JW
    77     REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
    78     REAL CMIN, CMAX ! Min and Max absolute ph. vel.
    79     REAL CPHA ! absolute PHASE VELOCITY frequency
    80     REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
    81     REAL ZP(NW)       ! Horizontal wavenumber angle       
    82     REAL ZO(NW, KLON) ! Absolute frequency      !
    83 
    84     ! Waves Intr. freq. at the 1/2 lev surrounding the full level
    85     REAL ZOM(NW, KLON), ZOP(NW, KLON)
    86 
    87     ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    88     REAL WWM(NW, KLON), WWP(NW, KLON)
    89     ! Fluxes X and Y for each waves at 1/2 Levels
    90     REAL RUWP(NW, KLON), RVWP(NW, KLON)
    91     REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
    92     REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
    93 
    94     REAL RUW0(NW, KLON) ! Fluxes at launching level
    95     REAL,SAVE :: RUWMAX         ! Max value
     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, parameter:: cmax = 30.   ! Max horizontal absolute phase velocity
     94    REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity   
     95    REAL CPHA                   ! absolute PHASE VELOCITY frequency
     96    REAL ZK(NW, ngrid)           ! Horizontal wavenumber amplitude
     97    REAL ZP(NW, ngrid)           ! Horizontal wavenumber angle       
     98    REAL ZO(NW, ngrid)           ! Absolute frequency
     99
     100    REAL intr_freq_m(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
     101    REAL intr_freq_p(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
     102    REAL wwm(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level below the full level
     103    REAL wwp(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level above the full level
     104    REAL u_epflux_p(nw, ngrid)           ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
     105    REAL v_epflux_p(nw, ngrid)           ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
     106    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) 
     107    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)
     108    REAL epflux_0(nw, ngrid)             ! Fluxes at launching level (previous name: RUW0)
     109    REAL, save :: epflux_max             ! Max EP flux value at launching altitude (previous name: RUWMAX)
    96110    INTEGER LAUNCH      ! Launching altitude
    97     REAL XLAUNCH        ! Control the launching altitude
    98         REAL ZMAXTH_TOP     ! Top of convective layer (approx.)
    99 
    100    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
    101 
    102     REAL SAT  ! saturation parameter
    103     REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
    104 
    105     ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
    106 
    107     REAL H0bis(KLON, KLEV) ! Characteristic Height of the atmosphere
    108     REAL H0 ! Characteristic Height of the atmosphere
    109     REAL PR, TR ! Reference Pressure and Temperature
    110 
    111     REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude (constant H0)
    112     REAL ZHbis(KLON, KLEV + 1) ! Log-pressure altitude (varying H)
    113 
    114     REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
    115     REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
    116     REAL PSEC ! Security to avoid division by 0 pressure
    117     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
    118     REAL BVSEC ! Security to avoid negative BVF
     111    REAL, parameter:: xlaunch = 0.4      ! Control the launching altitude
     112    REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.)
     113
     114
     115    REAL PREC(ngrid)
     116    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
     117
     118
     119    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
     120    REAL, parameter:: sat   = 1.     ! saturation parameter
     121    REAL, parameter:: rdiss = 1.     ! coefficient of dissipation
     122    REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency
     123
     124    ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
     125    REAL H0bis(ngrid, nlayer)          ! Characteristic Height of the atmosphere
     126    REAL, save:: H0                 ! Characteristic Height of the atmosphere
     127    REAL, parameter:: pr = 250      ! Reference pressure [Pa]
     128    REAL, parameter:: tr = 220.     ! Reference temperature [K]
     129    REAL ZH(ngrid, nlayer + 1)         ! Log-pressure altitude (constant H0)
     130    REAL ZHbis(ngrid, nlayer + 1)      ! Log-pressure altitude (varying H)
     131    REAL UH(ngrid, nlayer + 1)         ! zonal wind at 1/2 levels
     132    REAL VH(ngrid, nlayer + 1)         ! meridional wind at 1/2 levels
     133    REAL PH(ngrid, nlayer + 1)         ! Pressure at 1/2 levels
     134    REAL, parameter:: psec = 1.e-6  ! Security to avoid division by 0 pressure
     135    REAL BV(ngrid, nlayer + 1)         ! Brunt Vaisala freq. (BVF) at 1/2 levels
     136    REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV 
     137    REAL HREF(nlayer + 1)             ! Reference altitude for launching alt.
     138
    119139
    120140! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
    … …  
    125145    integer      ieq
    126146
    127 ! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
    128     real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
     147    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
     148
     149
    129150    LOGICAL,SAVE :: firstcall = .true.
    130151
    131 !      REAL ALEAS
    132 !      EXTERNAL ALEAS
    133152
    134153   !-----------------------------------------------------------------
    135154    !  1. INITIALISATIONS
    136155
    137       IF (firstcall) THEN
     156     IF (firstcall) THEN
    138157        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
    139         allocate(d_u_sav(NLON,NLEV),d_v_sav(NLON,NLEV))
    140         d_u_sav(:,:) = 0.
    141         d_v_sav(:,:) = 0.
    142         RUWMAX=1.E-4
    143         CALL getin_p("nonoro_gwd_RUWMAX",RUWMAX)
    144         write(*,*) "nonoro_gwd_ran: RUWMAX=",RUWMAX
    145         firstcall=.false.
    146       ENDIF
     158        epflux_max = 7.E-7 ! Mars' value !!
     159        call getin_p("nonoro_gwd_epflux_max", epflux_max)
     160        write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max
     161        ! Characteristic vertical scale height
     162        H0 = r * tr / g
     163        ! Control
     164        if (deltat .LT. dtime) THEN
     165             call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
     166        endif
     167        if (nlayer .LT. nw) THEN
     168             call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
     169        endif
     170        firstcall = .false.
     171     ENDIF
     172
     173    gwd_convective_source=.false.
    147174
    148175    ! Compute current values of temperature and winds
    … …  
    150177    uu(:,:)=pu(:,:)+dtime*pdu(:,:)
    151178    vv(:,:)=pv(:,:)+dtime*pdv(:,:)
    152    
    153     !    1.1 Basic parameter
    154 
    155     !  PARAMETERS CORRESPONDING TO V3:
    156 !     RUWMAX  = 1.e-4      ! best-case for MCS comparison
    157 !      RUWMAX = 5.e-6
    158 !     RUWMAX  = 1.e-3 
    159 
    160 !     SAT    = 0.25       ! Saturation parameter: Sc in (12)       
    161       SAT    = 1.         ! GG: This is the bestcase value for MCS comparison,
    162 !                         ! but runs for 1 full MY in the thermosphere are instables.                 
    163 !!!   TESTS                   
    164 !      SAT    = 0.85
    165 !
    166 !   RDISS  = 40         ! Diffusion parameter     
    167     RDISS  =  1
    168 !    RDISS = 0.5
    169 
    170     DELTAT=24.*3600.    ! Time scale of the waves (first introduced in 9b)
    171 !    DELTAT=DTIME
    172 
    173     KMIN = 2.E-5        ! Min horizontal wavenumber
    174     KMAX = 7.E-4        ! Max horizontal wavenumber
    175 
    176     !Online output: one value only
    177     if (output) then
    178       KMIN = 6.3E-5
    179       KMAX = 6.3E-5
    180     endif
    181     CMIN = 1.           ! Min phase velocity
    182 !    CMAX = 10.          ! Max phase speed velocity
    183     CMAX = 30.       
    184 !    CMAX = 60.
    185     XLAUNCH=0.4         ! Parameter that control launching altitude P/Ps
    186 !!! TEST: launching GW near the ground
    187 !    XLAUNCH =0.8
    188 
    189 ! MARS above typical convective cell at 250 Pa
    190     PR = 250.           ! Reference pressure  [Pa]   !
    191     TR = 220.           ! Reference Temperature [K]  !
    192 
    193     ZMAXTH_TOP  = 8000. ! Reference max altitude of convective layer [m]
    194          
    195     H0 = r * TR / g   ! Characteristic vertical scale height
    196 
    197     BVSEC  = 1.E-5      ! Security to avoid negative BVF 
    198     PSEC   = 1.E-6      ! Security to avoid division by 0 pressure
    199     ZOISEC = 1.E-6      ! Security FOR 0 INTRINSIC FREQ
    200 
    201     IF(DELTAT.LT.DTIME)THEN
    202        PRINT *,'GWD RANDO: DELTAT LT DTIME!'
    203        STOP
    204     ENDIF
    205 
    206 
    207     IF (NLEV < NW) THEN
    208        PRINT *, 'YOU WILL HAVE PROBLEM WITH RANDOM NUMBERS'
    209        PRINT *, 'FLOTT GWD STOP'
    210        STOP 1
    211     ENDIF
     179
     180
    212181
    213182    !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
    … …  
    218187
    219188    ! Pressure and Inv of pressure, Temperature / at 1/2 level
    220     DO LL = 2, KLEV
     189    DO LL = 2, nlayer
    221190       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
    222191    end DO
    223192
    224     PH(:, KLEV + 1) = 0.
     193    PH(:, nlayer + 1) = 0.
    225194    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
    226195
    227196    ! Launching altitude
    228197
    229     DO LL = 1, NLEV
    230        IF (PH(NLON / 2, LL) / PH(NLON / 2, 1) > XLAUNCH) LAUNCH = LL
    231 !       IF (plevmoy(LL) / plevmoy(1) > XLAUNCH) LAUNCH = LL   ! on venus
     198    !Pour revenir a la version non reproductible en changeant le nombre de
     199    !process
     200    ! Reprend la formule qui calcule PH en fonction de PP=play
     201    DO LL = 2, nlayer
     202       HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
     203    end DO
     204    HREF(nlayer + 1) = 0.
     205    HREF(1) = 2. * presnivs(1) - HREF(2)
     206
     207    LAUNCH=0
     208    DO LL = 1, nlayer
     209       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
    232210    ENDDO
    233211 
    234212    if (output) print*, " WE ARE IN FLOTT GW SCHEME "
    235213   
    236     !print*, "launch=",LAUNCH
    237     !print*,"launch p =",PH(NLON/2, LAUNCH)
    238 
    239214    ! Log pressure vert. coordinate
    240     DO LL = 1, KLEV + 1
     215    DO LL = 1, nlayer + 1
    241216       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
    242      !  print*,'ZH at each level',LL, ZH(NLON/2, LL)
    243217    end DO
    244218
    … …  
    246220    ! altitude above surface
    247221       ZHbis(:,1) = 0.   
    248        DO LL = 2, KLEV + 1
     222       DO LL = 2, nlayer + 1
    249223          H0bis(:, LL-1) = r * TT(:, LL-1) / g
    250224          ZHbis(:, LL) = ZHbis(:, LL-1) &
    … …  
    254228
    255229    ! Winds and BV frequency
    256     DO LL = 2, KLEV
     230    DO LL = 2, nlayer
    257231       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
    258232       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
    259233       ! GG test       
    260        !print*, 'TT, UH, VH, ZH at launch', TT(NLON/2,LAUNCH), UH(NLON/2,LAUNCH),VH(NLON/2, LAUNCH), ZH(NLON/2,LAUNCH) 
     234       !print*, 'TT, UH, VH, ZH at launch', TT(ngrid/2,LAUNCH), UH(ngrid/2,LAUNCH),VH(ngrid/2, LAUNCH), ZH(ngrid/2,LAUNCH)     
    261235       ! BVSEC: BV Frequency
    262236       BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
    … …  
    266240    end DO
    267241       !GG test
    268        !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(NLON/2, LAUNCH) 
     242       !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(ngrid/2, LAUNCH) 
    269243
    270244    BV(:, 1) = BV(:, 2)
    271245    UH(:, 1) = 0.
    272246    VH(:, 1) = 0.
    273     BV(:, KLEV + 1) = BV(:, KLEV)
    274     UH(:, KLEV + 1) = UU(:, KLEV)
    275     VH(:, KLEV + 1) = VV(:, KLEV)
     247    BV(:, nlayer + 1) = BV(:, nlayer)
     248    UH(:, nlayer + 1) = UU(:, nlayer)
     249    VH(:, nlayer + 1) = VV(:, nlayer)
    276250
    277251
    … …  
    283257    ! in a stochastic way
    284258
    285 !! A REVOIR:
    286 !! - utilisation de MOD ou bien de aleas ?
    287 !! - distribution gaussienne des CPHA ? (avec signe ZP qui est ajuste apres)
    288 
    289     JW = 0
    290     DO JP = 1, NP
    291        DO JK = 1, NK
    292           DO JO = 1, NO
    293              JW = JW + 1
     259    DO JW = 1, NW
    294260             !  Angle
    295              ZP(JW) = 2. * pi * REAL(JP - 1) / REAL(NP)
    296              DO II = 1, KLON
     261             DO II = 1, ngrid
     262                ! Angle (0 or PI so far)
     263                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
     264                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
     265                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
     266                     * PI / 2.
    297267                ! Horizontal wavenumber amplitude
    298 !                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
    299                 ZK(JW, II) = KMIN + (KMAX - KMIN) * ALEAS(0.)
     268                ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
    300269                ! Horizontal phase speed
    301 !                CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
    302                 CPHA = CMIN + (CMAX - CMIN) * ALEAS(0.)
     270                CPHA = 0.
     271                DO JJ = 1, NA
     272                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
     273                    CPHA = CPHA + &
     274                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
     275                END DO
     276                IF (CPHA.LT.0.)  THEN
     277                   CPHA = -1.*CPHA
     278                   ZP(JW,II) = ZP(JW,II) + PI
     279                ENDIF
    303280       !Online output: linear
    304281                if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1)
    … …  
    307284                ! Intrinsic frequency  is imposed
    308285                    ZO(JW, II) = ZO(JW, II)      &
    309                   + ZK(JW, II) * COS(ZP(JW)) * UH(II, LAUNCH) &
    310                   + ZK(JW, II) * SIN(ZP(JW)) * VH(II, LAUNCH)
     286                  + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
     287                  + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
    311288                ! Momentum flux at launch lev
    312 !!
    313 !!! CHANGE on MARS GCM: variable EP-FLUX according to:
    314 !     TEST1:   PBL diurnal variation
    315 !!
    316 !                 RUW0(JW, II) = MAX(RUWMAX / 100.,  &
    317 !                               zmax_therm(II)/ZMAXTH_TOP * RUWMAX)
    318          
    319 !!! BELOW: previous version from Flott original version
    320                 RUW0(JW, II) = RUWMAX &
    321                      * ALEAS(0.)     ! &
    322 !!! The factor cos**8:Limit the fluxes to tropical regions
    323 !                    * COS(pi/180.*latitude(II))**8
     289                ! epflux_0(JW, II) = epflux_max / REAL(NW) &
     290                epflux_0(JW, II) = epflux_max &
     291                     * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.)
    324292       !Online output: fixed to max
    325                 if (output) RUW0(JW, II) = RUWMAX
     293                if (output) epflux_0(JW, II) = epflux_max
    326294             ENDDO
    327           end DO
    328        end DO
    329     end DO
     295   end DO
    330296
    331297    ! 4. COMPUTE THE FLUXES
    … …  
    336302    !
    337303    DO JW = 1, NW
    338 
    339304       ! Evaluate intrinsic frequency at launching altitude:
    340        ZOP(JW, :) = ZO(JW, :) &
    341             - ZK(JW, :) * COS(ZP(JW)) * UH(:, LAUNCH) &
    342             - ZK(JW, :) * SIN(ZP(JW)) * VH(:, LAUNCH)
     305       intr_freq_p(JW, :) = ZO(JW, :) &
     306            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
     307            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
     308    end DO
     309
     310    IF (gwd_convective_source) THEN
     311         DO JW = 1, NW
     312       ! VERSION WITH CONVECTIVE SOURCE (designed for Earth)
     313
     314       ! Vertical velocity at launch level, value to ensure the
     315       ! imposed mmt flux factor related to the convective forcing:
     316       ! precipitations.
     317
     318       ! tanh limitation to values above prmax:
     319!       WWP(JW, :) = epflux_0(JW, :) &
     320!            * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
     321!       Here, we neglected the kinetic energy providing of the thermodynamic
     322!       phase change
     323
     324!
     325
     326       ! Factor related to the characteristics of the waves:
     327            WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:)  &
     328                 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**3
     329
     330      ! Moderation by the depth of the source (dz here):
     331            WWP(JW, :) = WWP(JW, :) &
     332                 * EXP(- BVLOW(:)**2 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
     333                 * ZK(JW, :)**2 * DZ**2)
     334
     335      ! Put the stress in the right direction:
     336            u_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
     337                 * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
     338            v_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
     339                 * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
     340         end DO
     341    ELSE ! VERSION WITHOUT CONVECTIVE SOURCE
    343342       ! Vertical velocity at launch level, value to ensure the imposed
    344343       ! mom flux:
    345        WWP(JW, :) = SQRT(ABS(ZOP(JW, :)) / MAX(BV(:, LAUNCH),BVSEC) &
    346             * RUW0(JW,:))
    347        RUWP(JW, :) = COS(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
    348        RVWP(JW, :) = SIN(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
    349 
    350     end DO
    351 
    352     !!!  EARTH VERSION  !!!!
    353     !  4.2 Uniform values below the launching altitude
    354 
    355     DO LL = 1, LAUNCH
    356        RUW(:, LL) = 0
    357        RVW(:, LL) = 0
    358        DO JW = 1, NW
    359           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    360           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
    361        end DO
    362     end DO
    363 
    364     !!!  VENUS VERSION  !!!!
     344         DO JW = 1, NW
     345       ! WW is directly a flux, here, not vertical velocity anymore
     346            WWP(JW, :) = epflux_0(JW,:)
     347            u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
     348            v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
     349
     350         end DO
     351    ENDIF
    365352    !  4.2 Initial flux at launching altitude
    366353
    367     !RUW(:, LAUNCH) = 0
    368     !RVW(:, LAUNCH) = 0
    369     !DO JW = 1, NW
    370     !   RUW(:, LAUNCH) = RUW(:, LAUNCH) + RUWP(JW, :)
    371     !   RVW(:, LAUNCH) = RVW(:, LAUNCH) + RVWP(JW, :)
    372     !end DO
     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
    373360
    374361    !  4.3 Loop over altitudes, with passage from one level to the
    … …  
    379366    !Online output
    380367    if (output) then
    381         ieq=nlon/2+1
     368        ieq=ngrid/2+1
    382369        write(str2,'(i2)') NW+2
    383370        outform="("//str2//"(E12.4,1X))"
    384371        WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
    385                (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW)), JW = 1, NW)
     372               (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW, IEQ)), JW = 1, NW)
    386373    endif
    387374
    388     DO LL = LAUNCH, KLEV - 1
     375    DO LL = LAUNCH, nlayer - 1
    389376
    390377
    … …  
    392379       ! TO THE NEXT)
    393380       DO JW = 1, NW
    394           ZOM(JW, :) = ZOP(JW, :)
     381          intr_freq_m(JW, :) = intr_freq_p(JW, :)
    395382          WWM(JW, :) = WWP(JW, :)
    396383          ! Intrinsic Frequency
    397           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW)) * UH(:, LL + 1) &
    398                - ZK(JW, :) * SIN(ZP(JW)) * VH(:, LL + 1)
     384          intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) &
     385               - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1)
    399386
    400387          WWP(JW, :) = MIN( &
    401388       ! No breaking (Eq.6)
    402389               WWM(JW, :) &
    403                * SQRT(BV(:, LL ) / BV(:, LL+1) &
    404                * ABS(ZOP(JW, :)) / MAX(ABS(ZOM(JW, :)), ZOISEC)) &
    405390      ! Dissipation (Eq. 8):
    406391               * EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
    407392               * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
    408                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
     393               / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 &
    409394               * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), &
    410395       ! Critical levels (forced to zero if intrinsic frequency changes sign)
    411                MAX(0., SIGN(1., ZOP(JW, :) * ZOM(JW, :))) &
     396               MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) &
    412397       ! Saturation (Eq. 12)
    413                * ABS(ZOP(JW, :))**3 /BV(:, LL+1) &
     398               * ABS(intr_freq_p(JW, :))**3 /BV(:, LL+1) &
    414399               * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4) 
    415400       end DO
    … …  
    419404       ! Give the right orientation to the stress
    420405       DO JW = 1, NW
    421           RUWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    422                *BV(:,LL+1)&
    423                * COS(ZP(JW)) *  WWP(JW, :)**2
    424           RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    425                *BV(:,LL+1)&
    426                * SIN(ZP(JW)) *  WWP(JW, :)**2
     406          u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) *  WWP(JW, :)
     407          v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) *  WWP(JW, :)
    427408       end DO
    428409 
    429        !
    430        RUW(:, LL + 1) = 0.
    431        RVW(:, LL + 1) = 0.
     410       u_epflux_tot(:, LL + 1) = 0.
     411       v_epflux_tot(:, LL + 1) = 0.
    432412
    433413       DO JW = 1, NW
    434           RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
    435           RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
     414          u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :)
     415          v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :)
     416          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW)
     417          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW)
    436418       end DO
    437419       !Online output
    438420       if (output) then
    439421         do JW=1,NW
    440             if(RUWP(JW, IEQ).gt.0.) then
    441               RUWP(JW, IEQ) = max(RUWP(JW, IEQ), 1.e-99)
     422            if(u_epflux_p(JW, IEQ).gt.0.) then
     423              u_epflux_p(JW, IEQ) = max(u_epflux_p(JW, IEQ), 1.e-99)
    442424            else
    443               RUWP(JW, IEQ) = min(RUWP(JW, IEQ), -1.e-99)
     425              u_epflux_p(JW, IEQ) = min(u_epflux_p(JW, IEQ), -1.e-99)
    444426            endif
    445427         enddo
    446428                   WRITE(11,outform) ZH(IEQ, LL+1) / 1000., &
    447429                                  ZHbis(IEQ, LL+1) / 1000., &
    448                                   (RUWP(JW, IEQ), JW = 1, NW)
     430                                  (u_epflux_p(JW, IEQ), JW = 1, NW)
    449431       endif
    450432
    … …  
    456438    ! 5.1 Rectification des flux au sommet et dans les basses couches:
    457439
    458     RUW(:, KLEV + 1) = 0.
    459     RVW(:, KLEV + 1) = 0.
    460     RUW(:, 1) = RUW(:, LAUNCH)
    461     RVW(:, 1) = RVW(:, LAUNCH)
    462     DO LL = 2, LAUNCH
    463        RUW(:, LL) = RUW(:, LL - 1) + (RUW(:, LAUNCH + 1) - RUW(:, 1)) * &
    464             (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
    465        RVW(:, LL) = RVW(:, LL - 1) + (RVW(:, LAUNCH + 1) - RVW(:, 1)) * &
    466             (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
    467     end DO
    468 
    469     ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
    470     DO LL = 1, KLEV
    471        d_u(:, LL) = G * (RUW(:, LL + 1) - RUW(:, LL)) &
     440! Attention, ici c'est le total sur toutes les ondes...
     441
     442    u_epflux_tot(:, nlayer + 1) = 0.
     443    v_epflux_tot(:, nlayer + 1) = 0.
     444
     445    ! Here, big change compared to FLott version:
     446    ! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
     447    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
     448    DO LL = 1, max(1,LAUNCH-3)
     449      u_epflux_tot(:, LL) = 0.
     450      v_epflux_tot(:, LL) = 0.
     451    end DO
     452    DO LL = max(2,LAUNCH-2), LAUNCH-1
     453       u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * &
     454            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     455       v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * &
     456            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     457       EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + &
     458            EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
     459            (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     460       WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + &
     461            WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
     462            (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     463    end DO
     464    ! This way, the total flux from GW is zero, but there is a net transport
     465    ! (upward) that should be compensated by circulation
     466    ! and induce additional friction at the surface
     467
     468    !Online output
     469    if (output) then
     470       DO LL = 1, nlayer - 1
     471           WRITE(11,*) ZHbis(IEQ, LL)/1000.,u_epflux_tot(IEQ,LL)
     472       end DO
     473       CLOSE(11)
     474       stop
     475    endif
     476
     477
     478    ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
     479    !---------------------------------------------
     480    DO LL = 1, nlayer
     481       d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
    472482            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
    473        d_v(:, LL) = G * (RVW(:, LL + 1) - RVW(:, LL)) &
     483       d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
    474484            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
    475485    ENDDO
    476486    d_t(:,:) = 0.
    477     ! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
    478     d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) + (1.-DTIME/DELTAT) * d_u_sav(:,:)
    479     d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) + (1.-DTIME/DELTAT) * d_v_sav(:,:)
    480     d_u_sav(:,:) = d_u(:,:)
    481     d_v_sav(:,:) = d_v(:,:)
     487
     488    ! 5.3 Update tendency of wind with the previous (and saved) values
     489    !-----------------------------------------------------------------
     490    d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:)                       &
     491             + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
     492    d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:)                       &
     493             + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
     494    du_nonoro_gwd(:,:) = d_u(:,:)
     495    dv_nonoro_gwd(:,:) = d_v(:,:)
    482496
    483497    ! Cosmetic: evaluation of the cumulated stress
    … …  
    485499    ZUSTR(:) = 0.
    486500    ZVSTR(:) = 0.
    487     DO LL = 1, KLEV
     501    DO LL = 1, nlayer
    488502       ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
    489503       ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
    … …  
    492506  END SUBROUTINE NONORO_GWD_RAN
    493507
    494 !===================================================================
    495 !===================================================================
    496 !===================================================================
    497 !===================================================================
    498 
    499   REAL FUNCTION ALEAS (R)
     508! ===================================================================
     509! Subroutines used to write variables of memory in start files       
     510! ===================================================================
     511
     512  SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
     513
    500514  IMPLICIT NONE
    501 !***BEGIN PROLOGUE  ALEAS
    502 !***PURPOSE  Generate a uniformly distributed random number.
    503 !***LIBRARY   SLATEC (FNLIB)
    504 !***CATEGORY  L6A21
    505 !***TYPE      SINGLE PRECISION (ALEAS-S)
    506 !***KEYWORDS  FNLIB, ALEAS NUMBER, SPECIAL FUNCTIONS, UNIFORM
    507 !***AUTHOR  Fullerton, W., (LANL)
    508 !***DESCRIPTION
    509 !
    510 !      This pseudo-random number generator is portable among a wide
    511 ! variety of computers.  RAND(R) undoubtedly is not as good as many
    512 ! readily available installation dependent versions, and so this
    513 ! routine is not recommended for widespread usage.  Its redeeming
    514 ! feature is that the exact same random numbers (to within final round-
    515 ! off error) can be generated from machine to machine.  Thus, programs
    516 ! that make use of random numbers can be easily transported to and
    517 ! checked in a new environment.
    518 !
    519 !      The random numbers are generated by the linear congruential
    520 ! method described, e.g., by Knuth in Seminumerical Methods (p.9),
    521 ! Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
    522 ! sequence, the I+1 -st number is generated from
    523 !             X(I+1) = (A*X(I) + C) MOD M,
    524 ! where here M = 2**22 = 4194304, C = 1731 and several suitable values
    525 ! of the multiplier A are discussed below.  Both the multiplier A and
    526 ! random number X are represented in double precision as two 11-bit
    527 ! words.  The constants are chosen so that the period is the maximum
    528 ! possible, 4194304.
    529 !
    530 !      In order that the same numbers be generated from machine to
    531 ! machine, it is necessary that 23-bit integers be reducible modulo
    532 ! 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
    533 ! integers be multiplied exactly.  Furthermore, if the restart option
    534 ! is used (where R is between 0 and 1), then the product R*2**22 =
    535 ! R*4194304 must be correct to the nearest integer.
    536 !
    537 !      The first four random numbers should be .0004127026,
    538 ! .6750836372, .1614754200, and .9086198807.  The tenth random number
    539 ! is .5527787209, and the hundredth is .3600893021 .  The thousandth
    540 ! number should be .2176990509 .
    541 !
    542 !      In order to generate several effectively independent sequences
    543 ! with the same generator, it is necessary to know the random number
    544 ! for several widely spaced calls.  The I-th random number times 2**22,
    545 ! where I=K*P/8 and P is the period of the sequence (P = 2**22), is
    546 ! still of the form L*P/8.  In particular we find the I-th random
    547 ! number multiplied by 2**22 is given by
    548 ! I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
    549 ! RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
    550 ! Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
    551 !
    552 !      Several multipliers have been subjected to the spectral test
    553 ! (see Knuth, p. 82).  Four suitable multipliers roughly in order of
    554 ! goodness according to the spectral test are
    555 !    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
    556 !    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
    557 !    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
    558 !    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
    559 !
    560 !      In the table below LOG10(NU(I)) gives roughly the number of
    561 ! random decimal digits in the random numbers considered I at a time.
    562 ! C is the primary measure of goodness.  In both cases bigger is better.
    563 !
    564 !                   LOG10 NU(I)              C(I)
    565 !       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
    566 !
    567 !    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
    568 !    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
    569 !    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
    570 !    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
    571 !   Best
    572 !    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
    573 !
    574 !             Input Argument --
    575 ! R      If R=0., the next random number of the sequence is generated.
    576 !        If R .LT. 0., the last generated number will be returned for
    577 !          possible use in a restart procedure.
    578 !        If R .GT. 0., the sequence of random numbers will start with
    579 !          the seed R mod 1.  This seed is also returned as the value of
    580 !          RAND provided the arithmetic is done exactly.
    581 !
    582 !             Output Value --
    583 ! RAND   a pseudo-random number between 0. and 1.
    584 !
    585 !***REFERENCES  (NONE)
    586 !***ROUTINES CALLED  (NONE)
    587 !***REVISION HISTORY  (YYMMDD)
    588 !   770401  DATE WRITTEN
    589 !   890531  Changed all specific intrinsics to generic.  (WRB)
    590 !   890531  REVISION DATE from Version 3.2
    591 !   891214  Prologue converted to Version 4.0 format.  (BAB)
    592 !***END PROLOGUE  RAND
    593       INTEGER IA1, IA0, IA1MA0, IC, IX1, IX0
    594       SAVE IA1, IA0, IA1MA0, IC, IX1, IX0
    595       DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
    596       DATA IC /1731/
    597       DATA IX1, IX0 /0, 0/
    598       REAL R
    599       INTEGER IY0,IY1
    600 !***FIRST EXECUTABLE STATEMENT  RAND
    601 !
    602 !           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
    603 !                   + IA0*IX0) + IA0*IX0
    604 !
    605       IF (R.EQ.0.) THEN
    606        IY0 = IA0*IX0
    607        IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
    608        IY0 = IY0 + IC
    609        IX0 = MOD (IY0, 2048)
    610        IY1 = IY1 + (IY0-IX0)/2048
    611        IX1 = MOD (IY1, 2048)
    612       ENDIF
    613 
    614       IF (R.GT.0.) THEN
    615        IX1 = MOD(R,1.)*4194304. + 0.5
    616        IX0 = MOD (IX1, 2048)
    617        IX1 = (IX1-IX0)/2048
    618       ENDIF
    619 
    620       ALEAS = IX1*2048 + IX0
    621       ALEAS = ALEAS / 4194304.
    622       RETURN
    623 
    624       END FUNCTION ALEAS
     515
     516      INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
     517      INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
     518
     519         allocate(du_nonoro_gwd(ngrid,nlayer))
     520         allocate(dv_nonoro_gwd(ngrid,nlayer))
     521
     522  END SUBROUTINE ini_nonoro_gwd_ran
     523! ----------------------------------
     524  SUBROUTINE end_nonoro_gwd_ran
     525
     526  IMPLICIT NONE
     527
     528         if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
     529         if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
     530
     531  END SUBROUTINE end_nonoro_gwd_ran
    625532
    626533END MODULE nonoro_gwd_ran_mod
Note: See TracChangeset for help on using the changeset viewer.