Changeset 5512


Ignore:
Timestamp:
Jan 28, 2025, 7:07:51 PM (2 days ago)
Author:
yann meurdesoif
Message:

Implement GPU automatic port for :

  • Thermics
  • acama_gwd_rando
  • flott_gwd_rando

YM

Location:
LMDZ6/trunk/libf/phylmd
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90

    r5309 r5512  
    55
    66  USE clesphys_mod_h
    7     implicit none
     7  IMPLICIT NONE
     8  LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
     9  LOGICAL, SAVE :: firstcall = .TRUE.
     10  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
     11 
     12  INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
     13  INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
     14
    815
    916contains
    10 
    11   SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
    12        zustr, zvstr, d_u, 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, A. de la Camara
    17     ! July, 24th, 2014
    18     ! Gaussian distribution of the source, source is vorticity squared
    19     ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
    20     ! Lott et al (JAS, 2010, vol 67, page 157-170)
    21     ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
    22 
    23 !  ONLINE:
    24 USE yoegwd_mod_h
    25         USE yomcst_mod_h
    26 use dimphy, only: klon, klev
    27     use assert_m, only: assert
    28     USE ioipsl_getin_p_mod, ONLY : getin_p
    29     USE vertical_layers_mod, ONLY : presnivs
    30 
    31 
    32 !  OFFLINE:
    33 !   include "dimensions_mod.f90"
    34 !   include "dimphy.h"
    35 !END DIFFERENCE
    36 
    37     ! 0. DECLARATIONS:
    38 
    39     ! 0.1 INPUTS
    40     REAL, intent(in)::DTIME ! Time step of the Physics
    41     REAL, intent(in):: PP(:, :) ! (KLON, KLEV) Pressure at full levels
    42     REAL, intent(in):: ROT(:,:) ! Relative vorticity             
    43     REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels
    44     REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels
    45     REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels
    46     REAL, intent(in):: PLAT(:) ! (KLON) LATITUDE                   
    47 
    48     ! 0.2 OUTPUTS
    49     REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses
    50 
    51     REAL, intent(inout):: d_u(:, :), d_v(:, :)
    52     REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress
    53     REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress
    54     ! (KLON, KLEV) tendencies on winds
    55 
    56     ! O.3 INTERNAL ARRAYS
    57     REAL BVLOW(klon)  !  LOW LEVEL BV FREQUENCY
    58     REAL ROTBA(KLON),CORIO(KLON)  !  BAROTROPIC REL. VORTICITY AND PLANETARY
    59     REAL UZ(KLON, KLEV + 1)
    60 
    61     INTEGER II, JJ, LL
    62 
    63     ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
    64 
    65     REAL DELTAT
    66 
    67     ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
    68 
    69     INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
    70     INTEGER JK, JP, JO, JW
    71     INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
    72     REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
    73     REAL CMIN, CMAX ! Min and Max absolute ph. vel.
    74     REAL CPHA ! absolute PHASE VELOCITY frequency
    75     REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
    76     REAL ZP(NW, KLON) ! Horizontal wavenumber angle
    77     REAL ZO(NW, KLON) ! Absolute frequency !
    78 
    79     ! Waves Intr. freq. at the 1/2 lev surrounding the full level
    80     REAL ZOM(NW, KLON), ZOP(NW, KLON)
    81 
    82     ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    83     REAL WWM(NW, KLON), WWP(NW, KLON)
    84 
    85     REAL RUW0(NW, KLON) ! Fluxes at launching level
    86 
    87     REAL RUWP(NW, KLON), RVWP(NW, KLON)
    88     ! Fluxes X and Y for each waves at 1/2 Levels
    89 
    90     INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
    91 
    92     REAL XLAUNCH ! Controle the launching altitude
    93     REAL XTROP ! SORT of Tropopause altitude
    94     REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
    95     REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
    96 
    97     REAL PRMAX ! Maximum value of PREC, and for which our linear formula
    98     ! for GWs parameterisation apply
    99 
    100     ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
    101 
    102     REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
    103     REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS
    104     REAL RUWFRT,SATFRT
    105 
    106     ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
    107 
    108     REAL H0 ! Characteristic Height of the atmosphere
    109     REAL DZ ! Characteristic depth of the source!
    110     REAL PR, TR ! Reference Pressure and Temperature
    111 
    112     REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
    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 PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
    118     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
    119     REAL BVSEC ! Security to avoid negative BVF
    120 
    121     REAL, DIMENSION(klev+1) ::HREF
    122     LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
    123     LOGICAL, SAVE :: firstcall = .TRUE.
    124   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
    125 
    126     CHARACTER (LEN=20) :: modname='acama_gwd_rando_m'
     17 
     18  SUBROUTINE ACAMA_GWD_rando_first
     19  use dimphy, only: klev
     20  USE ioipsl_getin_p_mod, ONLY : getin_p
     21  IMPLICIT NONE
     22    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
    12723    CHARACTER (LEN=80) :: abort_message
    128 
    129 
    130 
    131   IF (firstcall) THEN
     24 
     25    IF (firstcall) THEN
    13226    ! Cle introduite pour resoudre un probleme de non reproductibilite
    13327    ! Le but est de pouvoir tester de revenir a la version precedenete
    … …  
    14034    firstcall=.false.
    14135!    CALL iophys_ini(dtime)
    142   ENDIF
     36    ENDIF
     37  END SUBROUTINE ACAMA_GWD_rando_first
     38
     39  SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, &
     40       zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress)
     41
     42    ! Parametrization of the momentum flux deposition due to a discrete
     43    ! number of gravity waves.
     44    ! Author: F. Lott, A. de la Camara
     45    ! July, 24th, 2014
     46    ! Gaussian distribution of the source, source is vorticity squared
     47    ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )
     48    ! Lott et al (JAS, 2010, vol 67, page 157-170)
     49    ! Lott et al (JAS, 2012, vol 69, page 2134-2151)
     50
     51!  ONLINE:
     52USE yoegwd_mod_h
     53        USE yomcst_mod_h
     54use dimphy, only: klon, klev
     55    use assert_m, only: assert
     56    USE vertical_layers_mod, ONLY : presnivs
     57    USE lmdz_fake_call, ONLY : fake_call
     58
     59!  OFFLINE:
     60!   include "dimensions_mod.f90"
     61!   include "dimphy.h"
     62!END DIFFERENCE
     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):: ROT(KLON,KLEV) ! Relative vorticity             
     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    REAL, intent(in):: PLAT(KLON) ! (KLON) LATITUDE                   
     74
     75    ! 0.2 OUTPUTS
     76    REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses
     77
     78    REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV)
     79    REAL, intent(inout):: east_gwstress(KLON, KLEV) !  Profile of eastward stress
     80    REAL, intent(inout):: west_gwstress(KLON, KLEV) !  Profile of westward stress
     81    ! (KLON, KLEV) tendencies on winds
     82
     83    ! O.3 INTERNAL ARRAYS
     84    REAL BVLOW(klon)  !  LOW LEVEL BV FREQUENCY
     85    REAL ROTBA(KLON),CORIO(KLON)  !  BAROTROPIC REL. VORTICITY AND PLANETARY
     86    REAL UZ(KLON, KLEV + 1)
     87
     88    INTEGER II, JJ, LL
     89
     90    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
     91
     92    REAL DELTAT
     93
     94    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
     95
     96    INTEGER JK, JP, JO, JW
     97    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
     98    REAL CMIN, CMAX ! Min and Max absolute ph. vel.
     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    REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS
     129    REAL RUWFRT,SATFRT
     130
     131    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
     132
     133    REAL H0 ! Characteristic Height of the atmosphere
     134    REAL DZ ! Characteristic depth of the source!
     135    REAL PR, TR ! Reference Pressure and Temperature
     136
     137    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
     138
     139    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
     140    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
     141    REAL PSEC ! Security to avoid division by 0 pressure
     142    REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
     143    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
     144    REAL BVSEC ! Security to avoid negative BVF
     145
     146    REAL, DIMENSION(klev+1) ::HREF
     147
     148    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
     149    CHARACTER (LEN=80) :: abort_message
     150
    143151
    144152    !-----------------------------------------------------------------
    … …  
    211219!  END ONLINE
    212220
     221    CALL FAKE_CALL(BVLOW) ! to be suppress in future
     222    CALL FAKE_CALL(CORIO) ! to be suppress in future
     223    CALL FAKE_CALL(ROTBA) ! to be suppress in future
     224   
    213225    IF(DELTAT < DTIME)THEN
    214226!       PRINT *, 'flott_gwd_rando: deltat < dtime!'
    … …  
    276288            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
    277289    end DO
    278     BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
     290    BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
    279291         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
    280292         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
    … …  
    345357             DO II = 1, KLON
    346358                ! Angle (0 or PI so far)
    347                 ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
     359                ! ZP(II,JW) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
    348360                !      * RPI / 2.
    349361                ! Angle between 0 and pi
    350                   ZP(JW, II) = MOD(TT(II, JW) * 10., 1.) * RPI
     362                  ZP(II,JW) = MOD(TT(II, JW) * 10., 1.) * RPI
    351363! TEST WITH POSITIVE WAVES ONLY (Part I/II)
    352 !               ZP(JW, II) = 0.
     364!               ZP(II,JW) = 0.
    353365                ! Horizontal wavenumber amplitude
    354                 ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
     366                ZK(II,JW) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
    355367                ! Horizontal phase speed
    356368                CPHA = 0.
    … …  
    361373                IF (CPHA.LT.0.)  THEN
    362374                   CPHA = -1.*CPHA
    363                    ZP(JW,II) = ZP(JW,II) + RPI
     375                   ZP(II,JW) = ZP(II,JW) + RPI
    364376! TEST WITH POSITIVE WAVES ONLY (Part II/II)
    365 !               ZP(JW, II) = 0.
     377!               ZP(II,JW) = 0.
    366378                ENDIF
    367379                CPHA = CPHA + CMIN !we dont allow |c|<1m/s
    368380                ! Absolute frequency is imposed
    369                 ZO(JW, II) = CPHA * ZK(JW, II)
     381                ZO(II, JW) = CPHA * ZK(II,JW)
    370382                ! Intrinsic frequency is imposed
    371                 ZO(JW, II) = ZO(JW, II) &
    372                      + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
    373                      + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH)
     383                ZO(II, JW) = ZO(II, JW) &
     384                     + ZK(II,JW) * COS(ZP(II,JW)) * UH(II, LAUNCH) &
     385                     + ZK(II,JW) * SIN(ZP(II,JW)) * VH(II, LAUNCH)
    374386                ! Momentum flux at launch lev
    375387                ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE
    376388                !  RIGHT IN THE SH (GWD4 after 1990)
    377                   RUW0(JW, II) = 0.
     389                  RUW0(II, JW) = 0.
    378390                 DO JJ = 1, NA
    379                     RUW0(JW, II) = RUW0(JW,II) + &
     391                    RUW0(II, JW) = RUW0(II, JW) + &
    380392         2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.)
    381393                END DO
    382                 RUW0(JW, II) = RUWFRT &
    383                           * EXP(RUW0(JW,II))/1250. &  ! 2 mpa at south pole
     394                RUW0(II,JW) = RUWFRT &
     395                          * EXP(RUW0(II,JW))/1250. &  ! 2 mpa at south pole
    384396       *((1.05+SIN(PLAT(II)*RPI/180.))/(1.01+SIN(PLAT(II)*RPI/180.))-2.05/2.01)
    385                 ! RUW0(JW, II) = RUWFRT
     397                ! RUW0(II,JW) = RUWFRT
    386398             ENDDO
    387399    end DO
    … …  
    397409
    398410       ! Evaluate intrinsic frequency at launching altitude:
    399        ZOP(JW, :) = ZO(JW, :) &
    400             - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
    401             - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH)
     411       ZOP(:, JW) = ZO(:, JW) &
     412            - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LAUNCH) &
     413            - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LAUNCH)
    402414
    403415       ! VERSION WITH FRONTAL SOURCES
    … …  
    406418
    407419       ! tanh limitation for values above CORIO (inertial instability).
    408        ! WWP(JW, :) = RUW0(JW, :) &
    409        WWP(JW, :) = RUWFRT      &
     420       ! WWP(:,JW) = RUW0(:,JW) &
     421       WWP(:,JW) = RUWFRT      &
    410422       !     * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 &
    411423       !    * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) &
    … …  
    413425       !    * (CORIO(:)*CORIO(:)) &
    414426       ! MODERATION BY THE DEPTH OF THE SOURCE (DZ HERE)
    415        !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(JW, :)),ZOISEC)**2 &
    416        !      *ZK(JW, :)**2*DZ**2) &
     427       !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(:,JW)),ZOISEC)**2 &
     428       !      *ZK(:,JW)**2*DZ**2) &
    417429       ! COMPLETE FORMULA:
    418430            !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) &
    … …  
    420432       !  RESTORE DIMENSION OF A FLUX
    421433       !     *RD*TR/PR
    422        !     *1. + RUW0(JW, :)
     434       !     *1. + RUW0(:,JW)
    423435             *1.
    424436
    … …  
    429441       ! Put the stress in the right direction:
    430442
    431         RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    432         RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     443        RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     444        RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    433445
    434446    end DO
    … …  
    440452       RVW(:, LL) = 0
    441453       DO JW = 1, NW
    442           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    443           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
     454          RUW(:, LL) = RUW(:, LL) + RUWP(:,JW)
     455          RVW(:, LL) = RVW(:, LL) + RVWP(:,JW)
    444456       end DO
    445457    end DO
    … …  
    453465       ! to the next)
    454466       DO JW = 1, NW
    455           ZOM(JW, :) = ZOP(JW, :)
    456           WWM(JW, :) = WWP(JW, :)
     467          ZOM(:, JW) = ZOP(:,JW)
     468          WWM(:,JW) = WWP(:,JW)
    457469          ! Intrinsic Frequency
    458           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
    459                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
     470          ZOP(:,JW) = ZO(:, JW) - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LL + 1) &
     471               - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LL + 1)
    460472
    461473          ! No breaking (Eq.6)
    462474          ! Dissipation (Eq. 8)
    463           WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
     475          WWP(:,JW) = WWM(:,JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
    464476               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
    465                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
    466                * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
     477               / MAX(ABS(ZOP(:,JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
     478               * ZK(:,JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
    467479
    468480          ! Critical levels (forced to zero if intrinsic frequency changes sign)
    469481          ! Saturation (Eq. 12)
    470           WWP(JW, :) = min(WWP(JW, :), MAX(0., &
    471                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
     482          WWP(:,JW) = min(WWP(:,JW), MAX(0., &
     483               SIGN(1., ZOP(:,JW) * ZOM(:, JW))) * ABS(ZOP(:,JW))**3 &
    472484          !    / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 &
    473485               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 &
    474486!              *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 &
    475487               *SATFRT**2       &
    476                / ZK(JW, :)**4)
     488               / ZK(:,JW)**4)
    477489       end DO
    478490
    … …  
    481493
    482494       DO JW = 1, NW
    483           RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :)
    484           RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :)
     495          RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW)
     496          RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW)
    485497       end DO
    486498
    … …  
    489501
    490502       DO JW = 1, NW
    491           RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :)
    492           RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :)
    493           EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW)
    494           WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW)
     503          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:,JW)
     504          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:,JW)
     505          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:,JW))/REAL(NW)
     506          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:,JW))/REAL(NW)
    495507       end DO
    496508    end DO

  • LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90

    r5309 r5512  
    66  USE clesphys_mod_h
    77      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 
    678    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
    68     INTEGER JK, JP, JO, JW
    699    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 
    11910    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
    12011    LOGICAL, SAVE :: firstcall = .TRUE.
    121   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
    122 
    123 
    124   IF (firstcall) THEN
     12   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
     13
     14contains
     15
     16  SUBROUTINE FLOTT_GWD_rando_first
     17  use dimphy, only: klev
     18  USE ioipsl_getin_p_mod, ONLY : getin_p
     19  IMPLICIT NONE
     20    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
     21    CHARACTER (LEN=80) :: abort_message
     22 
     23    IF (firstcall) THEN
    12524    ! Cle introduite pour resoudre un probleme de non reproductibilite
    12625    ! Le but est de pouvoir tester de revenir a la version precedenete
    … …  
    13332    firstcall=.false.
    13433  ENDIF
     34  END SUBROUTINE FLOTT_GWD_rando_first
     35
     36
     37  SUBROUTINE FLOTT_GWD_rando(DTIME, PP, tt, uu, vv, prec, zustr, zvstr, d_u, &
     38       d_v,east_gwstress,west_gwstress)
     39
     40    ! Parametrization of the momentum flux deposition due to a discrete
     41    ! number of gravity waves.
     42    ! Author: F. Lott
     43    ! July, 12th, 2012
     44    ! Gaussian distribution of the source, source is precipitation
     45    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
     46
     47    !ONLINE:
     48      USE yomcst_mod_h
     49use dimphy, only: klon, klev
     50      use assert_m, only: assert
     51      USE ioipsl_getin_p_mod, ONLY : getin_p
     52      USE vertical_layers_mod, ONLY : presnivs
     53      USE yoegwd_mod_h
     54      USE lmdz_fake_call, ONLY : fake_call
     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   
     219    CALL FAKE_CALL(BVLOW)  ! to be suppress in future
     220   
    210221    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
    211222
    … …  
    292303                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
    293304                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
    294                 ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
     305                ZP(II, JW) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
    295306                     * RPI / 2.
    296307                ! Horizontal wavenumber amplitude
    297                 ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
     308                ZK(II, JW) = KMIN + (KMAX - KMIN) *RAN_NUM_2
    298309                ! Horizontal phase speed
    299310                CPHA = 0.
    … …  
    305316                IF (CPHA.LT.0.)  THEN
    306317                   CPHA = -1.*CPHA
    307                    ZP(JW,II) = ZP(JW,II) + RPI
     318                   ZP(II, JW) = ZP(II, JW) + RPI
    308319                ENDIF
    309320                ! Absolute frequency is imposed
    310                 ZO(JW, II) = CPHA * ZK(JW, II)
     321                ZO(II, JW) = CPHA * ZK(II, JW)
    311322                ! 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)
     323                ZO(II, JW) = ZO(II, JW) &
     324                     + ZK(II, JW) * COS(ZP(II, JW)) * UH(II, LAUNCH) &
     325                     + ZK(II, JW) * SIN(ZP(II, JW)) * VH(II, LAUNCH)
    315326                ! Momentum flux at launch lev
    316                 RUW0(JW, II) = RUWMAX
     327                RUW0(II, JW) = RUWMAX
    317328             ENDDO
    318329    ENDDO
    … …  
    326337
    327338       ! 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)
     339       ZOP(:,JW) = ZO(:, JW) &
     340            - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LAUNCH) &
     341            - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LAUNCH)
    331342
    332343       ! VERSION WITH CONVECTIVE SOURCE
    … …  
    337348
    338349       ! tanh limitation to values above prmax:
    339        WWP(JW, :) = RUW0(JW, :) &
     350       WWP(:, JW) = RUW0(:, JW) &
    340351            * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
    341352
    342353       ! 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
     354       WWP(:, JW) = WWP(:, JW) * ZK(:, JW)**3 / KMIN / BVLOW(:)  &
     355            / MAX(ABS(ZOP(:, JW)), ZOISEC)**3
    345356
    346357       ! 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 &
     358       WWP(:, JW) = WWP(:, JW) &
     359            * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 * ZK(:, JW)**2 &
    349360            * DZ**2)
    350361
    351362       ! 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
     363       RUWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
     364            * BV(:, LAUNCH) * COS(ZP(:, JW)) * WWP(:, JW)**2
     365       RVWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
     366            * BV(:, LAUNCH) * SIN(ZP(:, JW)) * WWP(:, JW)**2
    356367    end DO
    357368
    … …  
    363374       RVW(:, LL) = 0
    364375       DO JW = 1, NW
    365           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    366           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
     376          RUW(:, LL) = RUW(:, LL) + RUWP(:, JW)
     377          RVW(:, LL) = RVW(:, LL) + RVWP(:, JW)
    367378       end DO
    368379    end DO
    … …  
    376387       ! to the next)
    377388       DO JW = 1, NW
    378           ZOM(JW, :) = ZOP(JW, :)
    379           WWM(JW, :) = WWP(JW, :)
     389          ZOM(:, JW) = ZOP(:,JW)
     390          WWM(:, JW) = WWP(:, JW)
    380391          ! Intrinsic Frequency
    381           ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &
    382                - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1)
     392          ZOP(:, JW) = ZO(:, JW) - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LL + 1) &
     393               - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LL + 1)
    383394
    384395          ! No breaking (Eq.6)
    385396          ! Dissipation (Eq. 8)
    386           WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
     397          WWP(:, JW) = WWM(:, JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
    387398               + 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)))
     399               / MAX(ABS(ZOP(:, JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
     400               * ZK(:, JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
    390401
    391402          ! Critical levels (forced to zero if intrinsic frequency changes sign)
    392403          ! Saturation (Eq. 12)
    393           WWP(JW, :) = min(WWP(JW, :), MAX(0., &
    394                SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &
     404          WWP(:, JW) = min(WWP(:, JW), MAX(0., &
     405               SIGN(1., ZOP(:, JW) * ZOM(:, JW))) * ABS(ZOP(:, JW))**3 &
    395406               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2  &
    396                * SAT**2 / ZK(JW, :)**4)
     407               * SAT**2 / ZK(:, JW)**4)
    397408       end DO
    398409
    … …  
    401412
    402413       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, :)
     414          RUWP(:, JW) = SIGN(1., ZOP(:, JW))*COS(ZP(:, JW)) * WWP(:, JW)
     415          RVWP(:, JW) = SIGN(1., ZOP(:, JW))*SIN(ZP(:, JW)) * WWP(:, JW)
    405416       end DO
    406417
    … …  
    409420
    410421       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)
     422          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:, JW)
     423          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:, JW)
     424          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:, JW))/REAL(NW)
     425          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:, JW))/REAL(NW)
    415426       end DO
    416427    end DO

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_down.f90

    r5390 r5512  
    4343   integer ig,ilay
    4444   real, dimension(ngrid,nlay):: s1,s2,num !coefficients pour la resolution implicite
    45    integer :: iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
    46    
     45   integer :: iflag_impl ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
     46 
    4747   fdn(:,:)=0.
    4848   fup(:,:)=0.
    … …  
    5959   s2(:,:)=0.
    6060   num(:,:)=1.
     61   
     62   iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
    6163
    6264   if ( iflag_thermals_down < 10 ) then

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_dq.f90

    r5434 r5512  
    3838
    3939      integer niter,iter
    40       CHARACTER (LEN=20) :: modname='thermcell_dq'
     40      CHARACTER (LEN=20), PARAMETER :: modname='thermcell_dq'
    4141      CHARACTER (LEN=80) :: abort_message
    4242
    … …  
    190190      real ztimestep
    191191      integer niter,iter
    192       CHARACTER (LEN=20) :: modname='thermcell_dq'
     192      CHARACTER (LEN=20), PARAMETER :: modname='thermcell_dq'
    193193      CHARACTER (LEN=80) :: abort_message
    194194

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_dry.f90

    r5390 r5512  
    3333       REAL linter(ngrid),zlevinter(ngrid)
    3434       INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid)
    35       CHARACTER (LEN=20) :: modname='thermcell_dry'
    36       CHARACTER (LEN=80) :: abort_message
     35       CHARACTER (LEN=20), PARAMETER :: modname='thermcell_dry'
     36       CHARACTER (LEN=80) :: abort_message
    3737       INTEGER l,ig
    3838

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_env.f90

    r5390 r5512  
    5151! Calcul de l'humidite a saturation et de la condensation
    5252
    53    call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
     53   call thermcell_qsat(ngrid, nlay,mask,pplev,pt,po,pqsat)
    5454   do ll=1,nlay
    5555      do ig=1,ngrid

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_flux2.f90

    r5390 r5512  
    1616!---------------------------------------------------------------------------
    1717
    18       USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux
     18      USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux, thermals_fomass_max, thermals_alphamax
    1919      IMPLICIT NONE
    2020     
    … …  
    4848      REAL f_old,ddd0,eee0,ddd,eee,zzz
    4949
    50       REAL,SAVE :: fomass_max=0.5
    51       REAL,SAVE :: alphamax=0.7
    52 !$OMP THREADPRIVATE(fomass_max,alphamax)
    53 
    5450      logical check_debug,labort_physic
    5551
    56       character (len=20) :: modname='thermcell_flux2'
     52      character (len=20), PARAMETER :: modname='thermcell_flux2'
    5753      character (len=80) :: abort_message
    5854
    … …  
    391387        do ig=1,ngrid
    392388           if (zw2(ig,l+1).gt.1.e-10) then
    393            zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
     389           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*thermals_alphamax
    394390           if ( fm(ig,l+1) .gt. zfm) then
    395391              f_old=fm(ig,l+1)
    … …  
    430426            eee0=entr(ig,l)
    431427            ddd0=detr(ig,l)
    432             eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
     428            eee=entr(ig,l)-masse(ig,l)*thermals_fomass_max/ptimestep
    433429            ddd=detr(ig,l)-eee
    434430            if (eee.gt.0.) then
    … …  
    470466                         print*,'detr',detr(ig,l)
    471467                         print*,'masse',masse(ig,l)
    472                          print*,'fomass_max',fomass_max
    473                          print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
     468                         print*,'thermal_fomass_max',thermals_fomass_max
     469                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*thermals_fomass_max/ptimestep
    474470                         print*,'ptimestep',ptimestep
    475471                         print*,'lmax(ig)',lmax(ig)

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_ini.f90

    r5434 r5512  
    11MODULE lmdz_thermcell_ini
     2  USE strings_mod, ONLY : maxlen
    23
    34IMPLICIT NONE
    … …  
    3435integer, protected :: thermals_flag_alim=0       !
    3536integer, protected :: iflag_thermals_tenv=0      !
     37real,    protected :: thermals_fomass_max=0.5    ! Limitation du "vidage" des mailles sur un pas de temps 'thermcell_flux2'
     38real,    protected :: thermals_alphamax=0.7      ! fraction max des thermiques dans 'thermcell_flux2'
    3639
    3740   ! WARNING !!! fact_epsilon is not protected. It can be modified in thermcell_plume*
    … …  
    4750!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power)
    4851!$OMP THREADPRIVATE( mix0, thermals_flag_alim)
    49 !$OMP THREADPRIVATE(iflag_thermals_tenv)
     52!$OMP THREADPRIVATE(thermals_fomass_max)
     53!$OMP THREADPRIVATE(thermals_alphamax)
    5054
    5155integer, protected       :: thermals_subsid_advect_more_than_one=1
    52 character*6, protected :: thermals_subsid_advect_scheme = 'upwind' ! or 'center'
     56character(LEN=maxlen), protected :: thermals_subsid_advect_scheme = 'upwind' ! or 'center'
    5357
    5458!$OMP THREADPRIVATE(thermals_subsid_advect_scheme,thermals_subsid_advect_more_than_one)

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_main.F90

    r5390 r5512  
    140140
    141141
    142       integer,save :: igout=1
    143 !$OMP THREADPRIVATE(igout)
    144       integer,save :: lunout1=6
    145 !$OMP THREADPRIVATE(lunout1)
    146       integer,save :: lev_out=10
    147 !$OMP THREADPRIVATE(lev_out)
     142      integer, parameter :: igout=1
     143      integer, parameter :: lunout1=6
     144      integer, parameter :: lev_out=10
    148145
    149146      real lambda, zf,zf2,var,vardiff,CHI
    … …  
    166163      logical, dimension(ngrid,nlay) :: mask
    167164
    168       character (len=20) :: modname='thermcell_main'
     165      character (len=20), parameter :: modname='thermcell_main'
    169166      character (len=80) :: abort_message
    170167
    … …  
    191188       sorties=.true.
    192189      IF(ngrid.NE.ngrid) THEN
    193          PRINT*
    194190         PRINT*,'STOP dans convadj'
    195191         PRINT*,'ngrid    =',ngrid
    … …  
    240236        !    SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
    241237        ! &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
    242         ! contenu thermcell_env : call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
     238        ! contenu thermcell_env : call thermcell_qsat(ngrid, nlay,mask,pplev,pt,po,pqsat)
    243239        ! contenu thermcell_env : do ll=1,nlay
    244240        ! contenu thermcell_env :    do ig=1,ngrid
    … …  
    272268            enddo
    273269        enddo
    274         call thermcell_qsat(ngrid*nlay,mask,pplev,ptemp_env,p_o,zqsat)
     270        call thermcell_qsat(ngrid, nlay, mask,pplev,ptemp_env,p_o,zqsat)
    275271         
    276272      endif

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume.f90

    r5434 r5512  
    218218
    219219   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    220    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     220   call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    221221    do ig=1,ngrid
    222222!       print*,'active',active(ig),ig,l
    … …  
    351351
    352352   ztemp(:)=zpspsk(:,l)*ztla(:,l)
    353    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
     353   call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    354354   do ig=1,ngrid
    355355      if (activetmp(ig)) then

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume_6A.f90

    r5434 r5512  
    216216
    217217   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    218    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     218   call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    219219    do ig=1,ngrid
    220220!       print*,'active',active(ig),ig,l
    … …  
    556556
    557557   ztemp(:)=zpspsk(:,l)*ztla(:,l)
    558    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
     558   call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    559559   do ig=1,ngrid
    560560      if (activetmp(ig)) then
    … …  
    917917
    918918   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    919    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     919   call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    920920
    921921    do ig=1,ngrid
    … …  
    10051005
    10061006   ztemp(:)=zpspsk(:,l)*ztla(:,l)
    1007    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
     1007   call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    10081008
    10091009   do ig=1,ngrid

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_qsat.f90

    r5390 r5512  
    11MODULE lmdz_thermcell_qsat
     2
     3  REAL, PARAMETER :: DDT0=.01
     4
    25CONTAINS
    36
    4 subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
     7subroutine thermcell_qsat(klon, nlev, active,pplev,ztemp,zqta,zqsat)
    58USE yoethf_mod_h
    69  USE yomcst_mod_h
     10
     11
    712implicit none
    813
    … …  
    1621
    1722! Arguments
    18 INTEGER klon
    19 REAL zpspsk(klon),pplev(klon)
    20 REAL ztemp(klon),zqta(klon),zqsat(klon)
    21 LOGICAL active(klon)
     23INTEGER, INTENT(IN) :: klon
     24INTEGER, INTENT(IN) :: nlev  ! number of vertical to apply qsat
     25REAL zpspsk(klon, nlev),pplev(klon, nlev)
     26REAL ztemp(klon, nlev),zqta(klon,nlev),zqsat(klon,nlev)
     27LOGICAL active(klon, nlev)
    2228
    2329! Variables locales
    2430INTEGER ig,iter
    25 REAL Tbef(klon),DT(klon)
     31REAL Tbef(klon,nlev),DT(klon,nlev)
    2632REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
    2733logical Zsat
    2834REAL RLvCp
    2935
    30 REAL, SAVE :: DDT0=.01
    31 !$OMP THREADPRIVATE(DDT0)
    32 
    33 LOGICAL afaire(klon),tout_converge
    34 
     36LOGICAL afaire(klon, nlev),tout_converge
     37INTEGER :: l
    3538!====================================================================
    3639! INITIALISATIONS
    … …  
    3942RLvCp = RLVTT/RCPD
    4043tout_converge=.false.
    41 afaire(:)=.false.
    42 DT(:)=0.
     44afaire(:,:)=.false.
     45DT(:,:)=0.
    4346
    4447
    … …  
    4851! converge= false des que la convergence est atteinte.
    4952!====================================================================
    50 
    51 do ig=1,klon
    52    if (active(ig)) then
    53                Tbef(ig)=ztemp(ig)
    54                zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
    55                qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
     53do l=1,nlev
     54  do ig=1,klon
     55     if (active(ig,l)) then
     56               Tbef(ig,l)=ztemp(ig,l)
     57               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
     58               qsatbef= R2ES * FOEEW(Tbef(ig,l),zdelta)/pplev(ig,l)
    5659               qsatbef=MIN(0.5,qsatbef)
    5760               zcor=1./(1.-retv*qsatbef)
    5861               qsatbef=qsatbef*zcor
    59                qlbef=max(0.,zqta(ig)-qsatbef)
     62               qlbef=max(0.,zqta(ig,l)-qsatbef)
    6063               DT(ig) = 0.5*RLvCp*qlbef
    61                zqsat(ig)=qsatbef
     64               zqsat(ig,l)=qsatbef
    6265     endif
     66  enddo
    6367enddo
    64 
    6568! Traitement du cas ou il y a condensation mais faible
    6669! On ne condense pas mais on dit que le qsat est le qta
    67 do ig=1,klon
    68    if (active(ig)) then
    69       if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
    70          zqsat(ig)=zqta(ig)
    71       endif
    72    endif
     70do l=1,nlev
     71  do ig=1,klon
     72     if (active(ig,l)) then
     73       if (0.<abs(DT(ig,l)).and.abs(DT(ig,l))<=DDT0) then
     74           zqsat(ig,l)=zqta(ig,l)
     75        endif
     76     endif
     77  enddo
    7378enddo
    7479
    7580do iter=1,10
    76     afaire(:)=abs(DT(:)).gt.DDT0
    77     do ig=1,klon
    78                if (afaire(ig)) then
    79                  Tbef(ig)=Tbef(ig)+DT(ig)
    80                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
    81                  qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
     81    do l=1,nlev
     82      afaire(:,l)=abs(DT(:,l)).gt.DDT0
     83      do ig=1,klon
     84               if (afaire(ig,l)) then
     85                 Tbef(ig,l)=Tbef(ig,l)+DT(ig,l)
     86                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
     87                 qsatbef= R2ES * FOEEW(Tbef(ig,l),zdelta)/pplev(ig,l)
    8288                 qsatbef=MIN(0.5,qsatbef)
    8389                 zcor=1./(1.-retv*qsatbef)
    8490                 qsatbef=qsatbef*zcor
    85                  qlbef=zqta(ig)-qsatbef
    86                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
     91                 qlbef=zqta(ig,l)-qsatbef
     92                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig,l)))
    8793                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    8894                 zcor=1./(1.-retv*qsatbef)
    89                  dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
    90                  num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
     95                 dqsat_dT=FOEDE(Tbef(ig,l),zdelta,zcvm5,qsatbef,zcor)
     96                 num=-Tbef(ig,l)+ztemp(ig,l)+RLvCp*qlbef
    9197                 denom=1.+RLvCp*dqsat_dT
    92                  zqsat(ig) = qsatbef
    93                  DT(ig)=num/denom
     98                 zqsat(ig,l) = qsatbef
     99                 DT(ig,l)=num/denom
    94100               endif
     101      enddo
    95102    enddo
    96103enddo

  • LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

    r5491 r5512  
    928928      ALLOCATE(dv_gwd_rando(klon,klev),dv_gwd_front(klon,klev))
    929929      ALLOCATE(east_gwstress(klon,klev),west_gwstress(klon,klev))
    930       east_gwstress(:,:)=0 !ym missing init
    931       west_gwstress(:,:)=0 !ym missing init
     930      east_gwstress(:,:)=0. !ym missing init
     931      west_gwstress(:,:)=0. !ym missing init
    932932      ALLOCATE(d_t_hin(klon,klev))
    933933      ALLOCATE(d_q_ch4(klon,klev))

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r5506 r5512  
    2020! PLEASE try to follow this rule
    2121
    22     USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
     22    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando, ACAMA_GWD_rando_first
    2323    USE aero_mod
    2424    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
    … …  
    3232    USE dimphy
    3333    USE etat0_limit_unstruct_mod
    34     USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
     34    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando, FLOTT_GWD_rando_first
    3535    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
    3636    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
    … …  
    8686    USE lmdz_atke_turbulence_ini, ONLY : atke_ini
    8787    USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
     88    USE calltherm_mod, ONLY : calltherm
    8889    USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
    8990    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
    … …  
    36883689          ENDIF
    36893690          !>jyg
    3690           CALL calltherm(pdtphys &
     3691          CALL calltherm(itap, pdtphys &
    36913692               ,pplay,paprs,pphi,weak_inversion &
    36923693                        ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
    … …  
    49694970    IF (.not. ok_hines .and. ok_gwd_rando) then
    49704971       ! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
     4972       CALL acama_GWD_rando_first()
    49714973       CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
    49724974            v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
    … …  
    49874989
    49884990    IF (ok_gwd_rando) THEN
     4991       CALL FLOTT_GWD_rando_first()
    49894992       CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
    49904993            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
Note: See TracChangeset for help on using the changeset viewer.