Ignore:
Timestamp:
Jul 28, 2025, 7:23:15 PM (3 weeks ago)
Author:
aborella
Message:

Merge with trunk r5789

Location:
LMDZ6/branches/contrails
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails

  • LMDZ6/branches/contrails/libf/phylmd/flott_gwd_rando_m.f90

    r5309 r5791  
    22! $Id$
    33!
     4!$gpum horizontal klon
    45module FLOTT_GWD_rando_m
    56
    67  USE clesphys_mod_h
    78      implicit none
    8 
    9 contains
    10 
    11   SUBROUTINE FLOTT_GWD_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, &
    12        d_v,east_gwstress,west_gwstress)
    13 
    14     ! Parametrization of the momentum flux deposition due to a discrete
    15     ! number of gravity waves.
    16     ! Author: F. Lott
    17     ! July, 12th, 2012
    18     ! Gaussian distribution of the source, source is precipitation
    19     ! Reference: Lott (JGR, vol 118, page 8897, 2013)
    20 
    21     !ONLINE:
    22       USE yomcst_mod_h
    23 use dimphy, only: klon, klev
    24       use assert_m, only: assert
    25       USE ioipsl_getin_p_mod, ONLY : getin_p
    26       USE vertical_layers_mod, ONLY : presnivs
    27       USE yoegwd_mod_h
    28       CHARACTER (LEN=20) :: modname='flott_gwd_rando'
    29       CHARACTER (LEN=80) :: abort_message
    30 
    31     ! OFFLINE:
    32     ! include "dimensions_mod.f90"
    33     ! include "dimphy.h"
    34     ! END OF DIFFERENCE ONLINE-OFFLINE
    35 
    36     ! 0. DECLARATIONS:
    37 
    38     ! 0.1 INPUTS
    39     REAL, intent(in)::DTIME ! Time step of the Physics
    40     REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels
    41     REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s)
    42     REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
    43     REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
    44     REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
    45 
    46     ! 0.2 OUTPUTS
    47     REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
    48 
    49     REAL, intent(inout):: d_u(:, :), d_v(:, :)
    50     REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
    51     REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress
    52 
    53     ! (KLON, KLEV) tendencies on winds
    54 
    55     ! O.3 INTERNAL ARRAYS
    56     REAL BVLOW(klon)
    57     REAL DZ   !  Characteristic depth of the Source
    58 
    59     INTEGER II, JJ, LL
    60 
    61     ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
    62 
    63     REAL DELTAT
    64 
    65     ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
    66 
    679    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
    68     INTEGER JK, JP, JO, JW
    6910    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
    70     REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
    71     REAL CMAX ! standard deviation of the phase speed distribution
    72     REAL RUWMAX,SAT  ! ONLINE SPECIFIED IN run.def
    73     REAL CPHA ! absolute PHASE VELOCITY frequency
    74     REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
    75     REAL ZP(NW, KLON) ! Horizontal wavenumber angle
    76     REAL ZO(NW, KLON) ! Absolute frequency !
    77 
    78     ! Waves Intr. freq. at the 1/2 lev surrounding the full level
    79     REAL ZOM(NW, KLON), ZOP(NW, KLON)
    80 
    81     ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    82     REAL WWM(NW, KLON), WWP(NW, KLON)
    83 
    84     REAL RUW0(NW, KLON) ! Fluxes at launching level
    85 
    86     REAL RUWP(NW, KLON), RVWP(NW, KLON)
    87     ! Fluxes X and Y for each waves at 1/2 Levels
    88 
    89     INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
    90 
    91     REAL XLAUNCH ! Controle the launching altitude
    92     REAL XTROP ! SORT of Tropopause altitude
    93     REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
    94     REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
    95 
    96     REAL PRMAX ! Maximum value of PREC, and for which our linear formula
    97     ! for GWs parameterisation apply
    98 
    99     ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
    100 
    101     REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
    102 
    103     ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
    104 
    105     REAL H0 ! Characteristic Height of the atmosphere
    106     REAL PR, TR ! Reference Pressure and Temperature
    107 
    108     REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
    109 
    110     REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
    111     REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
    112     REAL PSEC ! Security to avoid division by 0 pressure
    113     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
    114     REAL BVSEC ! Security to avoid negative BVF
    115     REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
    116 
    117     REAL, DIMENSION(klev+1) ::HREF
    118 
    11911    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
    12012    LOGICAL, SAVE :: firstcall = .TRUE.
    121   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
    122 
    123 
    124   IF (firstcall) THEN
     13   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
     14
     15contains
     16
     17  SUBROUTINE FLOTT_GWD_rando_first
     18  use dimphy, only: klev
     19  USE ioipsl_getin_p_mod, ONLY : getin_p
     20  IMPLICIT NONE
     21    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
     22    CHARACTER (LEN=80) :: abort_message
     23 
     24    IF (firstcall) THEN
    12525    ! Cle introduite pour resoudre un probleme de non reproductibilite
    12626    ! Le but est de pouvoir tester de revenir a la version precedenete
    … …  
    13333    firstcall=.false.
    13434  ENDIF
     35  END SUBROUTINE FLOTT_GWD_rando_first
     36
     37
     38  SUBROUTINE FLOTT_GWD_rando(DTIME, PP, tt, uu, vv, prec, zustr, zvstr, d_u, &
     39       d_v,east_gwstress,west_gwstress)
     40
     41    ! Parametrization of the momentum flux deposition due to a discrete
     42    ! number of gravity waves.
     43    ! Author: F. Lott
     44    ! July, 12th, 2012
     45    ! Gaussian distribution of the source, source is precipitation
     46    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
     47
     48    !ONLINE:
     49      USE yomcst_mod_h
     50use dimphy, only: klon, klev
     51      use assert_m, only: assert
     52      USE ioipsl_getin_p_mod, ONLY : getin_p
     53      USE vertical_layers_mod, ONLY : presnivs
     54      USE yoegwd_mod_h
     55
     56      CHARACTER (LEN=20),PARAMETER :: modname='flott_gwd_rando'
     57      CHARACTER (LEN=80) :: abort_message
     58
     59    ! OFFLINE:
     60    ! include "dimensions_mod.f90"
     61    ! include "dimphy.h"
     62    ! END OF DIFFERENCE ONLINE-OFFLINE
     63
     64    ! 0. DECLARATIONS:
     65
     66    ! 0.1 INPUTS
     67    REAL, intent(in)::DTIME ! Time step of the Physics
     68    REAL, intent(in):: pp(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels
     69    REAL, intent(in):: prec(KLON) ! (klon) Precipitation (kg/m^2/s)
     70    REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels
     71    REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels
     72    REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels
     73
     74    ! 0.2 OUTPUTS
     75    REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses
     76
     77    REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV)
     78    REAL, intent(inout):: east_gwstress(KLON, KLEV) !  Profile of eastward stress
     79    REAL, intent(inout):: west_gwstress(KLON, KLEV) !  Profile of westward stress
     80
     81    ! (KLON, KLEV) tendencies on winds
     82
     83    ! O.3 INTERNAL ARRAYS
     84    REAL BVLOW(klon)
     85    REAL DZ   !  Characteristic depth of the Source
     86
     87    INTEGER II, JJ, LL
     88
     89    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
     90
     91    REAL DELTAT
     92
     93    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
     94
     95    INTEGER JK, JP, JO, JW
     96    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
     97    REAL CMAX ! standard deviation of the phase speed distribution
     98    REAL RUWMAX,SAT  ! ONLINE SPECIFIED IN run.def
     99    REAL CPHA ! absolute PHASE VELOCITY frequency
     100    REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude
     101    REAL ZP(KLON, NW) ! Horizontal wavenumber angle
     102    REAL ZO(KLON, NW) ! Absolute frequency !
     103
     104    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
     105    REAL ZOM(KLON, NW), ZOP(KLON, NW)
     106
     107    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
     108    REAL WWM(KLON, NW), WWP(KLON, NW)
     109
     110    REAL RUW0(KLON, NW) ! Fluxes at launching level
     111
     112    REAL RUWP(KLON, NW), RVWP(KLON, NW)
     113    ! Fluxes X and Y for each waves at 1/2 Levels
     114
     115    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
     116
     117    REAL XLAUNCH ! Controle the launching altitude
     118    REAL XTROP ! SORT of Tropopause altitude
     119    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
     120    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
     121
     122    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
     123    ! for GWs parameterisation apply
     124
     125    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
     126
     127    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
     128
     129    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
     130
     131    REAL H0 ! Characteristic Height of the atmosphere
     132    REAL PR, TR ! Reference Pressure and Temperature
     133
     134    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
     135
     136    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
     137    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
     138    REAL PSEC ! Security to avoid division by 0 pressure
     139    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
     140    REAL BVSEC ! Security to avoid negative BVF
     141    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
     142
     143    REAL, DIMENSION(klev+1) ::HREF
    135144
    136145
    … …  
    197206    !END ONLINE
    198207ENDIF
    199 
     208   
    200209    IF(DELTAT < DTIME)THEN
    201210       abort_message='flott_gwd_rando: deltat < dtime!'
    … …  
    207216       CALL abort_physic(modname,abort_message,1)
    208217    ENDIF
    209 
     218   
    210219    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
    211220
    … …  
    292301                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
    293302                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
    294                 ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
     303                ZP(II, JW) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
    295304                     * RPI / 2.
    296305                ! Horizontal wavenumber amplitude
    297                 ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
     306                ZK(II, JW) = KMIN + (KMAX - KMIN) *RAN_NUM_2
    298307                ! Horizontal phase speed
    299308                CPHA = 0.
    … …  
    305314                IF (CPHA.LT.0.)  THEN
    306315                   CPHA = -1.*CPHA
    307                    ZP(JW,II) = ZP(JW,II) + RPI
     316                   ZP(II, JW) = ZP(II, JW) + RPI
    308317                ENDIF
    309318                ! Absolute frequency is imposed
    310                 ZO(JW, II) = CPHA * ZK(JW, II)
     319                ZO(II, JW) = CPHA * ZK(II, JW)
    311320                ! Intrinsic frequency is imposed
    312                 ZO(JW, II) = ZO(JW, II) &
    313                      + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
    314                      + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
     321                ZO(II, JW) = ZO(II, JW) &
     322                     + ZK(II, JW) * COS(ZP(II, JW)) * UH(II, LAUNCH) &
     323                     + ZK(II, JW) * SIN(ZP(II, JW)) * VH(II, LAUNCH)
    315324                ! Momentum flux at launch lev
    316                 RUW0(JW, II) = RUWMAX
     325                RUW0(II, JW) = RUWMAX
    317326             ENDDO
    318327    ENDDO
    … …  
    326335
    327336       ! Evaluate intrinsic frequency at launching altitude:
    328        ZOP(JW, :) = ZO(JW, :) &
    329             - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
    330             - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
     337       ZOP(:,JW) = ZO(:, JW) &
     338            - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LAUNCH) &
     339            - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LAUNCH)
    331340
    332341       ! VERSION WITH CONVECTIVE SOURCE
    … …  
    337346
    338347       ! tanh limitation to values above prmax:
    339        WWP(JW, :) = RUW0(JW, :) &
     348       WWP(:, JW) = RUW0(:, JW) &
    340349            * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
    341350
    342351       ! Factor related to the characteristics of the waves:
    343        WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:)  &
    344             / MAX(ABS(ZOP(JW, :)), ZOISEC)**3
     352       WWP(:, JW) = WWP(:, JW) * ZK(:, JW)**3 / KMIN / BVLOW(:)  &
     353            / MAX(ABS(ZOP(:, JW)), ZOISEC)**3
    345354
    346355       ! Moderation by the depth of the source (dz here):
    347        WWP(JW, :) = WWP(JW, :) &
    348             * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 * ZK(JW, :)**2 &
     356       WWP(:, JW) = WWP(:, JW) &
     357            * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 * ZK(:, JW)**2 &
    349358            * DZ**2)
    350359
    351360       ! Put the stress in the right direction:
    352        RUWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    353             * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
    354        RVWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    355             * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
     361       RUWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
     362            * BV(:, LAUNCH) * COS(ZP(:, JW)) * WWP(:, JW)**2
     363       RVWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
     364            * BV(:, LAUNCH) * SIN(ZP(:, JW)) * WWP(:, JW)**2
    356365    end DO
    357366
    … …  
    363372       RVW(:, LL) = 0
    364373       DO JW = 1, NW
    365           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    366           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
     374          RUW(:, LL) = RUW(:, LL) + RUWP(:, JW)
     375          RVW(:, LL) = RVW(:, LL) + RVWP(:, JW)
    367376       end DO
    368377    end DO
    … …  
    376385       ! to the next)
    377386       DO JW = 1, NW
    378           ZOM(JW, :) = ZOP(JW, :)
    379           WWM(JW, :) = WWP(JW, :)
     387          ZOM(:, JW) = ZOP(:,JW)
     388          WWM(:, JW) = WWP(:, JW)
    380389          ! Intrinsic Frequency
    381           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
    382                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
     390          ZOP(:, JW) = ZO(:, JW) - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LL + 1) &
     391               - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LL + 1)
    383392
    384393          ! No breaking (Eq.6)
    385394          ! Dissipation (Eq. 8)
    386           WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
     395          WWP(:, JW) = WWM(:, JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
    387396               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
    388                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
    389                * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
     397               / MAX(ABS(ZOP(:, JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
     398               * ZK(:, JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
    390399
    391400          ! Critical levels (forced to zero if intrinsic frequency changes sign)
    392401          ! Saturation (Eq. 12)
    393           WWP(JW, :) = min(WWP(JW, :), MAX(0., &
    394                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
     402          WWP(:, JW) = min(WWP(:, JW), MAX(0., &
     403               SIGN(1., ZOP(:, JW) * ZOM(:, JW))) * ABS(ZOP(:, JW))**3 &
    395404               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2  &
    396                * SAT**2 / ZK(JW, :)**4)
     405               * SAT**2 / ZK(:, JW)**4)
    397406       end DO
    398407
    … …  
    401410
    402411       DO JW = 1, NW
    403           RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    404           RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     412          RUWP(:, JW) = SIGN(1., ZOP(:, JW))*COS(ZP(:, JW)) * WWP(:, JW)
     413          RVWP(:, JW) = SIGN(1., ZOP(:, JW))*SIN(ZP(:, JW)) * WWP(:, JW)
    405414       end DO
    406415
    … …  
    409418
    410419       DO JW = 1, NW
    411           RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
    412           RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
    413           EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW)
    414           WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW)
     420          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:, JW)
     421          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:, JW)
     422          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:, JW))/REAL(NW)
     423          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:, JW))/REAL(NW)
    415424       end DO
    416425    end DO
Note: See TracChangeset for help on using the changeset viewer.