flott_gwd_rando_m.F90

Go to the documentation of this file.

module flott_gwd_rando_m

  implicit none

contains

  SUBROUTINE flott_gwd_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, &
       d_v,east_gwstress,west_gwstress)

    ! Parametrization of the momentum flux deposition due to a discrete
    ! number of gravity waves. 
    ! Author: F. Lott
    ! July, 12th, 2012
    ! Gaussian distribution of the source, source is precipitation
    ! Reference: Lott (JGR, vol 118, page 8897, 2013)

    !ONLINE:
      use dimphy, only: klon, klev
      use assert_m, only: assert
      include "YOMCST.h"
      include "clesphys.h"
    ! OFFLINE:
    ! include "dimensions.h"
    ! include "dimphy.h"
    ! END OF DIFFERENCE ONLINE-OFFLINE
    include "YOEGWD.h"

    ! 0. DECLARATIONS:

    ! 0.1 INPUTS
    REAL, intent(in)::DTIME ! Time step of the Physics
    REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels
    REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s) 
    REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels 
    REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
    REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels

    ! 0.2 OUTPUTS
    REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses

    REAL, intent(inout):: d_u(:, :), d_v(:, :) 
    REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
    REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress 

    ! (KLON, KLEV) tendencies on winds

    ! O.3 INTERNAL ARRAYS
    REAL BVLOW(klon)
    REAL DZ   !  Characteristic depth of the Source

    INTEGER II, JJ, LL

    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED

    REAL DELTAT

    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS

    INTEGER, PARAMETER:: NK = 2, np = 2, no = 2, nw = nk * np * no
    INTEGER JK, JP, JO, JW
    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
    REAL CMAX ! standard deviation of the phase speed distribution
    REAL RUWMAX,SAT  ! ONLINE SPECIFIED IN run.def
    REAL CPHA ! absolute PHASE VELOCITY frequency
    REAL ZK(nw, klon) ! Horizontal wavenumber amplitude
    REAL ZP(nw, klon) ! Horizontal wavenumber angle 
    REAL ZO(nw, klon) ! Absolute frequency !

    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
    REAL ZOM(nw, klon), ZOP(nw, klon)

    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    REAL WWM(nw, klon), WWP(nw, klon)

    REAL RUW0(nw, klon) ! Fluxes at launching level

    REAL RUWP(nw, klon), RVWP(nw, klon)
    ! Fluxes X and Y for each waves at 1/2 Levels

    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude

    REAL XLAUNCH ! Controle the launching altitude
    REAL XTROP ! SORT of Tropopause altitude 
    REAL RUW(klon, klev + 1) ! Flux x at semi levels
    REAL RVW(klon, klev + 1) ! Flux y at semi levels

    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
    ! for GWs parameterisation apply

    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS

    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ

    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE

    REAL H0 ! Characteristic Height of the atmosphere
    REAL PR, TR ! Reference Pressure and Temperature

    REAL ZH(klon, klev + 1) ! Log-pressure altitude

    REAL UH(klon, klev + 1), VH(klon, klev + 1) ! Winds at 1/2 levels
    REAL PH(klon, klev + 1) ! Pressure at 1/2 levels
    REAL PSEC ! Security to avoid division by 0 pressure
    REAL PHM1(klon, klev + 1) ! 1/Press at 1/2 levels
    REAL BV(klon, klev + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
    REAL BVSEC ! Security to avoid negative BVF

    !-----------------------------------------------------------------

    ! 1. INITIALISATIONS

    ! 1.1 Basic parameter

    ! Are provided from elsewhere (latent heat of vaporization, dry
    ! gaz constant for air, gravity constant, heat capacity of dry air
    ! at constant pressure, earth rotation rate, pi).

    ! 1.2 Tuning parameters of V14

    
    rdiss = 1. ! Diffusion parameter
    ! ONLINE 
      ruwmax=gwd_rando_ruwmax
      sat=gwd_rando_sat
    !END ONLINE
    ! OFFLINE
    ! RUWMAX= 1.75    ! Launched flux
    ! SAT=0.25     ! Saturation parameter
    ! END OFFLINE

    prmax = 20. / 24. /3600.
    ! maximum of rain for which our theory applies (in kg/m^2/s)

 ! Characteristic depth of the source
    dz = 1000.
    xlaunch=0.5 ! Parameter that control launching altitude
    xtrop=0.2 ! Parameter that control tropopause altitude
    deltat=24.*3600. ! Time scale of the waves (first introduced in 9b)
    !  OFFLINE
    !  DELTAT=DTIME
    !  END OFFLINE

    kmin = 2.e-5
    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)

    kmax = 1.e-3 ! Max horizontal wavenumber
    cmax = 30. ! Max phase speed velocity

    tr = 240. ! Reference Temperature
    pr = 101300. ! Reference pressure
    h0 = rd * tr / rg ! Characteristic vertical scale height

    bvsec = 5.e-3 ! Security to avoid negative BVF 
    psec = 1.e-6 ! Security to avoid division by 0 pressure
    zoisec = 1.e-6 ! Security FOR 0 INTRINSIC FREQ

    !ONLINE
        call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
         size(vv, 1), size(zustr), size(zvstr), size(d_u, 1), &
         size(d_v, 1), &
         size(east_gwstress, 1), size(west_gwstress, 1) /), &
         "FLOTT_GWD_RANDO klon")
     call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
          size(vv, 2), size(d_u, 2), size(d_v, 2), &
          size(east_gwstress,2), size(west_gwstress,2) /), &
          "FLOTT_GWD_RANDO klev")
    !END ONLINE

    IF(deltat < dtime)THEN
       print *, 'flott_gwd_rando: deltat < dtime!'
       stop 1
    ENDIF

    IF (klev < nw) THEN
       print *, 'flott_gwd_rando: you will have problem with random numbers'
       stop 1
    ENDIF

    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS

    ! Pressure and Inv of pressure
    DO ll = 2, klev
       ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.)
       phm1(:, ll) = 1. / ph(:, ll) 
    end DO

    ph(:, klev + 1) = 0. 
    phm1(:, klev + 1) = 1. / psec
    ph(:, 1) = 2. * pp(:, 1) - ph(:, 2)

    ! Launching altitude

    launch=0
    ltrop =0
    DO ll = 1, klev
       IF (ph(klon / 2, ll) / ph(klon / 2, 1) > xlaunch) launch = ll
    ENDDO
    DO ll = 1, klev
       IF (ph(klon / 2, ll) / ph(klon / 2, 1) > xtrop) ltrop = ll
    ENDDO

    ! Log pressure vert. coordinate
    DO ll = 1, klev + 1 
       zh(:, ll) = h0 * log(pr / (ph(:, ll) + psec))
    end DO

    ! BV frequency
    DO ll = 2, klev
       ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
       uh(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1)) &
            * rd**2 / rcpd / h0**2 + (tt(:, ll) &
            - tt(:, ll - 1)) / (zh(:, ll) - zh(:, ll - 1)) * rd / h0
    end DO
    bvlow(:) = 0.5 * (tt(:, ltrop )+ tt(:, launch)) &
         * rd**2 / rcpd / h0**2 + (tt(:, ltrop ) &
         - tt(:, launch))/(zh(:, ltrop )- zh(:, launch)) * rd / h0

    uh(:, 1) = uh(:, 2)
    uh(:, klev + 1) = uh(:, klev)
    bv(:, 1) = uh(:, 2)
    bv(:, klev + 1) = uh(:, klev)
    ! SMOOTHING THE BV HELPS
    DO ll = 2, klev
       bv(:, ll)=(uh(:, ll+1)+2.*uh(:, ll)+uh(:, ll-1))/4.
    end DO

    bv=max(sqrt(max(bv, 0.)), bvsec)
    bvlow=max(sqrt(max(bvlow, 0.)), bvsec)


    ! WINDS
    DO ll = 2, klev
       uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind
       vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind
    end DO
    uh(:, 1) = 0.
    vh(:, 1) = 0.
    uh(:, klev + 1) = uu(:, klev)
    vh(:, klev + 1) = vv(:, klev)

    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE

    ! The mod functions of weird arguments are used to produce the
    ! waves characteristics in an almost stochastic way

    jw = 0
    DO jp = 1, np
       DO jk = 1, nk
          DO jo = 1, no
             jw = jw + 1
             ! Angle
             DO ii = 1, klon
                ! Angle (0 or PI so far)
                zp(jw, ii) = (sign(1., 0.5 - mod(tt(ii, jw) * 10., 1.)) + 1.) &
                     * rpi / 2.
                ! Horizontal wavenumber amplitude
                zk(jw, ii) = kmin + (kmax - kmin) * mod(tt(ii, jw) * 100., 1.)
                ! Horizontal phase speed
                cpha = 0.
                DO jj = 1, na
                    cpha = cpha + &
         cmax*2.*(mod(tt(ii, jw+3*jj)**2, 1.)-0.5)*sqrt(3.)/sqrt(na*1.)
                END DO
                IF (cpha.LT.0.)  THEN
                   cpha = -1.*cpha
                   zp(jw,ii) = zp(jw,ii) + rpi
                ENDIF
                ! Absolute frequency is imposed
                zo(jw, ii) = cpha * zk(jw, ii) 
                ! Intrinsic frequency is imposed
                zo(jw, ii) = zo(jw, ii) &
                     + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) &
                     + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch)
                ! Momentum flux at launch lev 
                ruw0(jw, ii) = ruwmax
             ENDDO
          end DO
       end DO
    end DO

    ! 4. COMPUTE THE FLUXES

    ! 4.1 Vertical velocity at launching altitude to ensure 
    ! the correct value to the imposed fluxes.

    DO jw = 1, nw

       ! Evaluate intrinsic frequency at launching altitude:
       zop(jw, :) = zo(jw, :) &
            - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) &
            - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) 

       ! VERSION WITH CONVECTIVE SOURCE

       ! Vertical velocity at launch level, value to ensure the
       ! imposed factor related to the convective forcing:
       ! precipitations.

       ! tanh limitation to values above prmax:
       wwp(jw, :) = ruw0(jw, :) &
            * (rd / rcpd / h0 * rlvtt * prmax * tanh(prec(:) / prmax))**2

       ! Factor related to the characteristics of the waves:
       wwp(jw, :) = wwp(jw, :) * zk(jw, :)**3 / kmin / bvlow(:)  &
            / max(abs(zop(jw, :)), zoisec)**3 

       ! Moderation by the depth of the source (dz here):
       wwp(jw, :) = wwp(jw, :) &
            * exp(- bvlow(:)**2 / max(abs(zop(jw, :)), zoisec)**2 * zk(jw, :)**2 &
            * dz**2)

       ! Put the stress in the right direction:
       ruwp(jw, :) = zop(jw, :) / max(abs(zop(jw, :)), zoisec)**2 &
            * bv(:, launch) * cos(zp(jw, :)) * wwp(jw, :)**2
       rvwp(jw, :) = zop(jw, :) / max(abs(zop(jw, :)), zoisec)**2 &
            * bv(:, launch) * sin(zp(jw, :)) * wwp(jw, :)**2
    end DO


    ! 4.2 Uniform values below the launching altitude

    DO ll = 1, launch
       ruw(:, ll) = 0
       rvw(:, ll) = 0
       DO jw = 1, nw
          ruw(:, ll) = ruw(:, ll) + ruwp(jw, :)
          rvw(:, ll) = rvw(:, ll) + rvwp(jw, :)
       end DO
    end DO

    ! 4.3 Loop over altitudes, with passage from one level to the next
    ! done by i) conserving the EP flux, ii) dissipating a little,
    ! iii) testing critical levels, and vi) testing the breaking.

    DO ll = launch, klev - 1
       ! Warning: all the physics is here (passage from one level
       ! to the next)
       DO jw = 1, nw
          zom(jw, :) = zop(jw, :)
          wwm(jw, :) = wwp(jw, :)
          ! Intrinsic Frequency
          zop(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) &
               - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1) 

          ! No breaking (Eq.6)
          ! Dissipation (Eq. 8)
          wwp(jw, :) = wwm(jw, :) * exp(- 2. * rdiss * pr / (ph(:, ll + 1) &
               + ph(:, ll)) * ((bv(:, ll + 1) + bv(:, ll)) / 2.)**3 &
               / max(abs(zop(jw, :) + zom(jw, :)) / 2., zoisec)**4 &
               * zk(jw, :)**3 * (zh(:, ll + 1) - zh(:, ll)))

          ! Critical levels (forced to zero if intrinsic frequency changes sign)
          ! Saturation (Eq. 12)
          wwp(jw, :) = min(wwp(jw, :), max(0., &
               sign(1., zop(jw, :) * zom(jw, :))) * abs(zop(jw, :))**3 &
               / bv(:, ll + 1) * exp(- zh(:, ll + 1) / h0) * kmin**2  &
               * sat**2 / zk(jw, :)**4)
       end DO

       ! Evaluate EP-flux from Eq. 7 and give the right orientation to
       ! the stress

       DO jw = 1, nw
          ruwp(jw, :) = sign(1., zop(jw, :))*cos(zp(jw, :)) * wwp(jw, :)
          rvwp(jw, :) = sign(1., zop(jw, :))*sin(zp(jw, :)) * wwp(jw, :)
       end DO

       ruw(:, ll + 1) = 0.
       rvw(:, ll + 1) = 0.

       DO jw = 1, nw
          ruw(:, ll + 1) = ruw(:, ll + 1) + ruwp(jw, :) 
          rvw(:, ll + 1) = rvw(:, ll + 1) + rvwp(jw, :) 
          east_gwstress(:, ll)=east_gwstress(:, ll)+max(0.,ruwp(jw,:))/float(nw)
          west_gwstress(:, ll)=west_gwstress(:, ll)+min(0.,ruwp(jw,:))/float(nw)
       end DO
    end DO
! OFFLINE ONLY
!   PRINT *,'SAT PROFILE:'
!   DO LL=1,KLEV
!   PRINT *,ZH(KLON/2,LL)/1000.,SAT*(2.+TANH(ZH(KLON/2,LL)/H0-8.))
!   ENDDO

    ! 5 CALCUL DES TENDANCES:

    ! 5.1 Rectification des flux au sommet et dans les basses couches

    ruw(:, klev + 1) = 0.
    rvw(:, klev + 1) = 0.
    ruw(:, 1) = ruw(:, launch)
    rvw(:, 1) = rvw(:, launch)
    DO ll = 1, launch
       ruw(:, ll) = ruw(:, launch+1)
       rvw(:, ll) = rvw(:, launch+1)
       east_gwstress(:, ll)  = east_gwstress(:, launch)
       west_gwstress(:, ll)  = west_gwstress(:, launch)
    end DO

    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
    DO ll = 1, klev
       d_u(:, ll) = (1.-dtime/deltat) * d_u(:, ll) + dtime/deltat/REAL(NW) * &
            RG * (ruw(:, ll + 1) - ruw(:, ll)) &
            / (ph(:, ll + 1) - ph(:, ll)) * DTIME
       ! NO AR-1 FOR MERIDIONAL TENDENCIES
       d_v(:, ll) =                                            1./REAL(NW) * &
            RG * (rvw(:, ll + 1) - rvw(:, ll)) &
            / (ph(:, ll + 1) - ph(:, ll)) * DTIME
    ENDDO

    ! Cosmetic: evaluation of the cumulated stress
    zustr = 0.
    zvstr = 0.
    DO ll = 1, klev
       zustr = zustr + d_u(:, ll) / rg * (ph(:, ll + 1) - ph(:, ll))/dtime
       zvstr = zvstr + d_v(:, ll) / rg * (ph(:, ll + 1) - ph(:, ll))/dtime
    ENDDO

  END SUBROUTINE flott_gwd_rando

end module flott_gwd_rando_m