Ignore:
Timestamp:
Jul 29, 2024, 5:47:53 PM (3 months ago)
Author:
abarral
Message:

Put YOEGWD.h, FCTTRE.h into modules

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/coef_diff_turb_mod.F90

    r5139 r5143  
    1 
    21MODULE coef_diff_turb_mod
    32
    4 ! This module contains some procedures for calculation of the coefficients of the
    5 ! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion
    6 ! at surface(cdrag)
     3  ! This module contains some procedures for calculation of the coefficients of the
     4  ! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion
     5  ! at surface(cdrag)
    76
    87  IMPLICIT NONE
    9  
     8
    109CONTAINS
    1110
    12 !****************************************************************************************
     11  !****************************************************************************************
    1312
    1413  SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
    15        ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
    16        ycoefm, ycoefh ,yq2, yeps, ydrgpro)
    17  
     14          ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
     15          ycoefm, ycoefh, yq2, yeps, ydrgpro)
     16
    1817    USE dimphy
    1918    USE indice_sol_mod
    … …  
    2120    USE lmdz_clesphys
    2221    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
    23 
    24 ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
    25 ! atmosphere
    26 ! NB! No values are calculated between surface and the first model layer.
    27 !     ycoefm(:,1) and ycoefh(:,1) are not valid !!!
    28 
    29 
    30 ! Input arguments
    31 !****************************************************************************************
    32     REAL, INTENT(IN)                           :: dtime
    33     INTEGER, INTENT(IN)                        :: nsrf, knon
    34     INTEGER, DIMENSION(klon), INTENT(IN)       :: ni
    35     REAL, DIMENSION(klon,klev+1), INTENT(IN)   :: ypaprs
    36     REAL, DIMENSION(klon,klev), INTENT(IN)     :: ypplay
    37     REAL, DIMENSION(klon,klev), INTENT(IN)     :: yu, yv
    38     REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
    39     REAL, DIMENSION(klon), INTENT(IN)          :: yts, yqsurf
    40     REAL, DIMENSION(klon), INTENT(IN)          :: ycdragm
    41 !FC
    42     REAL, DIMENSION(klon,klev), INTENT(IN)     :: ydrgpro
    43 
    44 
    45 ! InOutput arguments
    46 !****************************************************************************************
    47     REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2
    48  
    49 ! Output arguments
    50 !****************************************************************************************
    51     REAL, DIMENSION(klon,klev+1), INTENT(OUT)  :: yeps
    52     REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefh
    53     REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
    54 
    55 ! Other local variables
    56 !****************************************************************************************
    57     INTEGER                                    :: k, i, j
    58     REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
    59     REAL, DIMENSION(klon,klev+1)               :: yzlev, q2diag, ykmm, ykmn, ykmq
    60     REAL, DIMENSION(klon)                      :: yustar
    61 
    62 ! Include
    63 !****************************************************************************************
    64     INCLUDE "YOETHF.h"
     22    USE lmdz_YOETHF
     23
     24    ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
     25    ! atmosphere
     26    ! NB! No values are calculated between surface and the first model layer.
     27    !     ycoefm(:,1) and ycoefh(:,1) are not valid !!!
     28
     29
     30    ! Input arguments
     31    !****************************************************************************************
     32    REAL, INTENT(IN) :: dtime
     33    INTEGER, INTENT(IN) :: nsrf, knon
     34    INTEGER, DIMENSION(klon), INTENT(IN) :: ni
     35    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: ypaprs
     36    REAL, DIMENSION(klon, klev), INTENT(IN) :: ypplay
     37    REAL, DIMENSION(klon, klev), INTENT(IN) :: yu, yv
     38    REAL, DIMENSION(klon, klev), INTENT(IN) :: yq, yt
     39    REAL, DIMENSION(klon), INTENT(IN) :: yts, yqsurf
     40    REAL, DIMENSION(klon), INTENT(IN) :: ycdragm
     41    !FC
     42    REAL, DIMENSION(klon, klev), INTENT(IN) :: ydrgpro
     43
     44
     45    ! InOutput arguments
     46    !****************************************************************************************
     47    REAL, DIMENSION(klon, klev + 1), INTENT(INOUT) :: yq2
     48
     49    ! Output arguments
     50    !****************************************************************************************
     51    REAL, DIMENSION(klon, klev + 1), INTENT(OUT) :: yeps
     52    REAL, DIMENSION(klon, klev), INTENT(OUT) :: ycoefh
     53    REAL, DIMENSION(klon, klev), INTENT(OUT) :: ycoefm
     54
     55    ! Other local variables
     56    !****************************************************************************************
     57    INTEGER :: k, i, j
     58    REAL, DIMENSION(klon, klev) :: ycoefm0, ycoefh0, yzlay, yteta
     59    REAL, DIMENSION(klon, klev + 1) :: yzlev, q2diag, ykmm, ykmn, ykmq
     60    REAL, DIMENSION(klon) :: yustar
     61
     62    ! Include
     63    !****************************************************************************************
    6564    INCLUDE "YOMCST.h"
    66 
    6765
    6866    ykmm = 0 !ym missing init
    6967    ykmn = 0 !ym missing init
    7068    ykmq = 0 !ym missing init
    71    
    72     yeps(:,:) = 0.
    73 
    74 !****************************************************************************************   
    75 ! Calcul de coefficients de diffusion turbulent de l'atmosphere :
    76 ! ycoefm(:,2:klev), ycoefh(:,2:klev)
    77 
    78 !****************************************************************************************   
     69
     70    yeps(:, :) = 0.
     71
     72    !****************************************************************************************
     73    ! Calcul de coefficients de diffusion turbulent de l'atmosphere :
     74    ! ycoefm(:,2:klev), ycoefh(:,2:klev)
     75
     76    !****************************************************************************************
    7977
    8078    CALL coefkz(nsrf, knon, ypaprs, ypplay, &
    81          ksta, ksta_ter, &
    82          yts, yu, yv, yt, yq, &
    83          yqsurf, &
    84          ycoefm, ycoefh)
    85  
    86 !****************************************************************************************
    87 ! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere :
    88 ! ycoefm(:,2:klev), ycoefh(:,2:klev)
    89 
    90 !****************************************************************************************
     79            ksta, ksta_ter, &
     80            yts, yu, yv, yt, yq, &
     81            yqsurf, &
     82            ycoefm, ycoefh)
     83
     84    !****************************************************************************************
     85    ! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere :
     86    ! ycoefm(:,2:klev), ycoefh(:,2:klev)
     87
     88    !****************************************************************************************
    9189
    9290    IF (iflag_pbl==1) THEN
    93        CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
    94             ycoefm0, ycoefh0)
    95 
    96        DO k = 2, klev
    97           DO i = 1, knon
    98              ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
    99              ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
    100           ENDDO
    101        ENDDO
     91      CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
     92              ycoefm0, ycoefh0)
     93
     94      DO k = 2, klev
     95        DO i = 1, knon
     96          ycoefm(i, k) = MAX(ycoefm(i, k), ycoefm0(i, k))
     97          ycoefh(i, k) = MAX(ycoefh(i, k), ycoefh0(i, k))
     98        ENDDO
     99      ENDDO
    102100    ENDIF
    103101
    104  
    105 !**************************************************************************************** 
    106 ! Calcul d'une diffusion minimale pour les conditions tres stables
    107 
    108 !****************************************************************************************
     102
     103    !****************************************************************************************
     104    ! Calcul d'une diffusion minimale pour les conditions tres stables
     105
     106    !****************************************************************************************
    109107    IF (ok_kzmin) THEN
    110        CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, &
    111             ycoefm0,ycoefh0)
    112        
    113        DO k = 2, klev
    114           DO i = 1, knon
    115              ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
    116              ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
    117           ENDDO
    118        ENDDO
    119        
     108      CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycdragm, &
     109              ycoefm0, ycoefh0)
     110
     111      DO k = 2, klev
     112        DO i = 1, knon
     113          ycoefm(i, k) = MAX(ycoefm(i, k), ycoefm0(i, k))
     114          ycoefh(i, k) = MAX(ycoefh(i, k), ycoefh0(i, k))
     115        ENDDO
     116      ENDDO
     117
    120118    ENDIF
    121119
    122  
    123 !****************************************************************************************
    124 ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
    125 
    126 !****************************************************************************************
     120
     121    !****************************************************************************************
     122    ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
     123
     124    !****************************************************************************************
    127125
    128126    IF (iflag_pbl>=3) THEN
    129127
    130        yzlay(1:knon,1)= &
    131             RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) &
    132             *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
    133        DO k=2,klev
    134           DO i = 1, knon
    135              yzlay(i,k)= &
    136                   yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) &
    137                   /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
    138           END DO
    139        END DO
    140 
    141        DO k=1,klev
    142           DO i = 1, knon
    143              yteta(i,k)= &
    144                   yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA &
    145                   *(1.+0.61*yq(i,k))
    146           END DO
    147        END DO
    148 
    149        yzlev(1:knon,1)=0.
    150        yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
    151        DO k=2,klev
    152           DO i = 1, knon
    153              yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1))
    154           END DO
    155        END DO
    156 
    157 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    158 !!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
    159 !!$! bug sur les coefficients de surface :
    160 !!$!          ycdragh(1:knon) = ycoefm(1:knon,1)
    161 !!$!          ycdragm(1:knon) = ycoefh(1:knon,1)
    162 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    163 
    164        ! Normalement, on peut passer dans les codes avec knon=0
    165        ! Mais ca fait planter le replay.
    166        ! En attendant une réécriture, on a joute des if (Fredho)
    167        IF ( klon>1 .OR. (klon==1 .AND. knon==1) ) THEN
    168           CALL ustarhb(knon,klev,knon,yu,yv,ycdragm, yustar)
    169        endif
    170      
    171        IF (prt_level > 9) THEN
    172           WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon)
    173        ENDIF
    174          
    175 !   iflag_pbl peut etre utilise comme longuer de melange
    176        IF (iflag_pbl>=31) THEN
    177           IF ( klon>1 .OR. (klon==1 .AND. knon==1) ) THEN
    178           CALL vdif_kcay(knon,klev,knon,dtime,RG,RD,ypaprs,yt, &
    179                yzlev,yzlay,yu,yv,yteta, &
    180                ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
    181                iflag_pbl)
    182           endif
    183        ELSE IF (iflag_pbl<20) THEN
    184           CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &
    185                yzlev,yzlay,yu,yv,yteta, &
    186                ycdragm,yq2,yeps,ykmm,ykmn,ykmq,yustar, &
    187                iflag_pbl,ydrgpro)
    188 !FC
    189        ENDIF
    190        
    191        ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
    192        ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
    193                
     128      yzlay(1:knon, 1) = &
     129              RD * yt(1:knon, 1) / (0.5 * (ypaprs(1:knon, 1) + ypplay(1:knon, 1))) &
     130                      * (ypaprs(1:knon, 1) - ypplay(1:knon, 1)) / RG
     131      DO k = 2, klev
     132        DO i = 1, knon
     133          yzlay(i, k) = &
     134                  yzlay(i, k - 1) + RD * 0.5 * (yt(i, k - 1) + yt(i, k)) &
     135                          / ypaprs(i, k) * (ypplay(i, k - 1) - ypplay(i, k)) / RG
     136        END DO
     137      END DO
     138
     139      DO k = 1, klev
     140        DO i = 1, knon
     141          yteta(i, k) = &
     142                  yt(i, k) * (ypaprs(i, 1) / ypplay(i, k))**RKAPPA &
     143                          * (1. + 0.61 * yq(i, k))
     144        END DO
     145      END DO
     146
     147      yzlev(1:knon, 1) = 0.
     148      yzlev(1:knon, klev + 1) = 2. * yzlay(1:knon, klev) - yzlay(1:knon, klev - 1)
     149      DO k = 2, klev
     150        DO i = 1, knon
     151          yzlev(i, k) = 0.5 * (yzlay(i, k) + yzlay(i, k - 1))
     152        END DO
     153      END DO
     154
     155      !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     156      !!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un
     157      !!$! bug sur les coefficients de surface :
     158      !!$!          ycdragh(1:knon) = ycoefm(1:knon,1)
     159      !!$!          ycdragm(1:knon) = ycoefh(1:knon,1)
     160      !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     161
     162      ! Normalement, on peut passer dans les codes avec knon=0
     163      ! Mais ca fait planter le replay.
     164      ! En attendant une réécriture, on a joute des if (Fredho)
     165      IF (klon>1 .OR. (klon==1 .AND. knon==1)) THEN
     166        CALL ustarhb(knon, klev, knon, yu, yv, ycdragm, yustar)
     167      endif
     168
     169      IF (prt_level > 9) THEN
     170        WRITE(lunout, *) 'USTAR = ', (yustar(i), i = 1, knon)
     171      ENDIF
     172
     173      !   iflag_pbl peut etre utilise comme longuer de melange
     174      IF (iflag_pbl>=31) THEN
     175        IF (klon>1 .OR. (klon==1 .AND. knon==1)) THEN
     176          CALL vdif_kcay(knon, klev, knon, dtime, RG, RD, ypaprs, yt, &
     177                  yzlev, yzlay, yu, yv, yteta, &
     178                  ycdragm, yq2, q2diag, ykmm, ykmn, yustar, &
     179                  iflag_pbl)
     180        endif
     181      ELSE IF (iflag_pbl<20) THEN
     182        CALL yamada4(ni, nsrf, knon, dtime, RG, RD, ypaprs, yt, &
     183                yzlev, yzlay, yu, yv, yteta, &
     184                ycdragm, yq2, yeps, ykmm, ykmn, ykmq, yustar, &
     185                iflag_pbl, ydrgpro)
     186        !FC
     187      ENDIF
     188
     189      ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
     190      ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
     191
    194192    ELSE
    195        ! No TKE for Standard Physics
    196        yq2=0.
     193      ! No TKE for Standard Physics
     194      yq2 = 0.
    197195    ENDIF !(iflag_pbl.ge.3)
    198196
    199197  END SUBROUTINE coef_diff_turb
    200198
    201 !****************************************************************************************
     199  !****************************************************************************************
    202200
    203201  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
    204        ksta, ksta_ter, &
    205        ts, &
    206        u,v,t,q, &
    207        qsurf, &
    208        pcfm, pcfh)
    209    
     202          ksta, ksta_ter, &
     203          ts, &
     204          u, v, t, q, &
     205          qsurf, &
     206          pcfm, pcfh)
     207
    210208    USE dimphy
    211209    USE indice_sol_mod
    212210    USE lmdz_print_control, ONLY: prt_level, lunout
    213211    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
    214  
    215 !======================================================================
    216 ! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
    217 !           (une version strictement identique a l'ancien modele)
    218 ! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
    219 !        coefficients d'echange turbulent dans l'atmosphere.
    220 ! Arguments:
    221 ! nsrf-----input-I- indicateur de la nature du sol
    222 ! knon-----input-I- nombre de points a traiter
    223 ! paprs----input-R- pregssion a chaque intercouche (en Pa)
    224 ! pplay----input-R- pression au milieu de chaque couche (en Pa)
    225 ! ts-------input-R- temperature du sol (en Kelvin)
    226 ! u--------input-R- vitesse u
    227 ! v--------input-R- vitesse v
    228 ! t--------input-R- temperature (K)
    229 ! q--------input-R- vapeur d'eau (kg/kg)
    230 
    231 ! pcfm-----output-R- coefficients a calculer (vitesse)
    232 ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
    233 !======================================================================
    234     INCLUDE "YOETHF.h"
     212    USE lmdz_YOETHF
     213    USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
     214
     215    !======================================================================
     216    ! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
     217    !           (une version strictement identique a l'ancien modele)
     218    ! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
     219    !        coefficients d'echange turbulent dans l'atmosphere.
     220    ! Arguments:
     221    ! nsrf-----input-I- indicateur de la nature du sol
     222    ! knon-----input-I- nombre de points a traiter
     223    ! paprs----input-R- pregssion a chaque intercouche (en Pa)
     224    ! pplay----input-R- pression au milieu de chaque couche (en Pa)
     225    ! ts-------input-R- temperature du sol (en Kelvin)
     226    ! u--------input-R- vitesse u
     227    ! v--------input-R- vitesse v
     228    ! t--------input-R- temperature (K)
     229    ! q--------input-R- vapeur d'eau (kg/kg)
     230
     231    ! pcfm-----output-R- coefficients a calculer (vitesse)
     232    ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
     233    !======================================================================
    235234    INCLUDE "YOMCST.h"
    236     INCLUDE "FCTTRE.h"
    237 
    238 ! Arguments:
    239 
    240     INTEGER, INTENT(IN)                      :: knon, nsrf
    241     REAL, INTENT(IN)                         :: ksta, ksta_ter
    242     REAL, DIMENSION(klon), INTENT(IN)        :: ts
    243     REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
    244     REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
    245     REAL, DIMENSION(klon,klev), INTENT(IN)   :: u, v, t, q
    246     REAL, DIMENSION(klon), INTENT(IN)        :: qsurf
    247 
    248     REAL, DIMENSION(klon,klev), INTENT(OUT)  :: pcfm, pcfh
    249 
    250 ! Local variables:
    251 
    252     INTEGER, DIMENSION(klon)    :: itop ! numero de couche du sommet de la couche limite
    253 
    254 ! Quelques constantes et options:
    255 
    256     REAL, PARAMETER :: cepdu2=0.1**2
    257     REAL, PARAMETER :: CKAP=0.4
    258     REAL, PARAMETER :: cb=5.0
    259     REAL, PARAMETER :: cc=5.0
    260     REAL, PARAMETER :: cd=5.0
    261     REAL, PARAMETER :: clam=160.0
    262     REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau
    263     LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide
    264     REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
    265     REAL, PARAMETER :: prandtl=0.4
     235
     236    ! Arguments:
     237
     238    INTEGER, INTENT(IN) :: knon, nsrf
     239    REAL, INTENT(IN) :: ksta, ksta_ter
     240    REAL, DIMENSION(klon), INTENT(IN) :: ts
     241    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
     242    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
     243    REAL, DIMENSION(klon, klev), INTENT(IN) :: u, v, t, q
     244    REAL, DIMENSION(klon), INTENT(IN) :: qsurf
     245
     246    REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
     247
     248    ! Local variables:
     249
     250    INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite
     251
     252    ! Quelques constantes et options:
     253
     254    REAL, PARAMETER :: cepdu2 = 0.1**2
     255    REAL, PARAMETER :: CKAP = 0.4
     256    REAL, PARAMETER :: cb = 5.0
     257    REAL, PARAMETER :: cc = 5.0
     258    REAL, PARAMETER :: cd = 5.0
     259    REAL, PARAMETER :: clam = 160.0
     260    REAL, PARAMETER :: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
     261    LOGICAL, PARAMETER :: richum = .TRUE. ! utilise le nombre de Richardson humide
     262    REAL, PARAMETER :: ric = 0.4 ! nombre de Richardson critique
     263    REAL, PARAMETER :: prandtl = 0.4
    266264    REAL kstable ! diffusion minimale (situation stable)
    267265    ! GKtest
    268266    ! PARAMETER (kstable=1.0e-10)
    269 !IM: 261103     REAL kstable_ter, kstable_sinon
    270 !IM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
    271 !IM: 261103     PARAMETER (kstable_ter = 1.0e-8)
    272 !IM: 261103   PARAMETER (kstable_ter = 1.0e-10)
    273 !IM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
     267    !IM: 261103     REAL kstable_ter, kstable_sinon
     268    !IM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
     269    !IM: 261103     PARAMETER (kstable_ter = 1.0e-8)
     270    !IM: 261103   PARAMETER (kstable_ter = 1.0e-10)
     271    !IM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
    274272    ! fin GKtest
    275     REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
     273    REAL, PARAMETER :: mixlen = 35.0 ! constante controlant longueur de melange
    276274    INTEGER isommet ! le sommet de la couche limite
    277     LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
    278     LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere
    279 
    280 ! Variables locales:
     275    LOGICAL, PARAMETER :: tvirtu = .TRUE. ! calculer Ri d'une maniere plus performante
     276    LOGICAL, PARAMETER :: opt_ec = .FALSE.! formule du Centre Europeen dans l'atmosphere
     277
     278    ! Variables locales:
    281279    INTEGER i, k !IM 120704
    282     REAL zgeop(klon,klev)
     280    REAL zgeop(klon, klev)
    283281    REAL zmgeom(klon)
    284282    REAL zri(klon)
    … …  
    288286    REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
    289287    REAL z2geomf, zalh2, zalm2, zscfh, zscfm
    290     REAL, PARAMETER :: t_coup=273.15
    291     LOGICAL, PARAMETER :: check=.FALSE.
    292 
    293 ! contre-gradient pour la chaleur sensible: Kelvin/metre
     288    REAL, PARAMETER :: t_coup = 273.15
     289    LOGICAL, PARAMETER :: check = .FALSE.
     290
     291    ! contre-gradient pour la chaleur sensible: Kelvin/metre
    294292    REAL gamt(2:klev)
    295293
    296     LOGICAL, SAVE :: appel1er=.TRUE.
     294    LOGICAL, SAVE :: appel1er = .TRUE.
    297295    !$OMP THREADPRIVATE(appel1er)
    298296
    299 ! Fonctions thermodynamiques et fonctions d'instabilite
     297    ! Fonctions thermodynamiques et fonctions d'instabilite
    300298    REAL fsta, fins, x
    301299
    302     fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
    303     fins(x) = SQRT(1.0-18.0*x)
    304 
    305     isommet=klev
    306      
     300    fsta(x) = 1.0 / (1.0 + 10.0 * x * (1 + 8.0 * x))
     301    fins(x) = SQRT(1.0 - 18.0 * x)
     302
     303    isommet = klev
     304
    307305    IF (appel1er) THEN
    308        IF (prt_level > 9) THEN
    309           WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
    310           WRITE(lunout,*)'coefkz, richum:', richum
    311           IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
    312           WRITE(lunout,*)'coefkz, isommet:', isommet
    313           WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
    314           appel1er = .FALSE.
    315        ENDIF
     306      IF (prt_level > 9) THEN
     307        WRITE(lunout, *)'coefkz, opt_ec:', opt_ec
     308        WRITE(lunout, *)'coefkz, richum:', richum
     309        IF (richum) WRITE(lunout, *)'coefkz, ratqs:', ratqs
     310        WRITE(lunout, *)'coefkz, isommet:', isommet
     311        WRITE(lunout, *)'coefkz, tvirtu:', tvirtu
     312        appel1er = .FALSE.
     313      ENDIF
    316314    ENDIF
    317315
    318 ! Initialiser les sorties
     316    ! Initialiser les sorties
    319317
    320318    DO k = 1, klev
    321        DO i = 1, knon
    322           pcfm(i,k) = 0.0
    323           pcfh(i,k) = 0.0
    324        ENDDO
     319      DO i = 1, knon
     320        pcfm(i, k) = 0.0
     321        pcfh(i, k) = 0.0
     322      ENDDO
    325323    ENDDO
    326324    DO i = 1, knon
    327        itop(i) = 0
    328     ENDDO
    329 
    330 ! Prescrire la valeur de contre-gradient
     325      itop(i) = 0
     326    ENDDO
     327
     328    ! Prescrire la valeur de contre-gradient
    331329
    332330    IF (iflag_pbl==1) THEN
    333        DO k = 3, klev
    334           gamt(k) = -1.0E-03
    335        ENDDO
    336        gamt(2) = -2.5E-03
     331      DO k = 3, klev
     332        gamt(k) = -1.0E-03
     333      ENDDO
     334      gamt(2) = -2.5E-03
    337335    ELSE
    338        DO k = 2, klev
    339           gamt(k) = 0.0
    340        ENDDO
     336      DO k = 2, klev
     337        gamt(k) = 0.0
     338      ENDDO
    341339    ENDIF
    342 !IM cf JLD/ GKtest
    343     IF ( nsrf /= is_oce ) THEN
    344 !IM 261103     kstable = kstable_ter
    345        kstable = ksta_ter
     340    !IM cf JLD/ GKtest
     341    IF (nsrf /= is_oce) THEN
     342      !IM 261103     kstable = kstable_ter
     343      kstable = ksta_ter
    346344    ELSE
    347 !IM 261103     kstable = kstable_sinon
    348        kstable = ksta
     345      !IM 261103     kstable = kstable_sinon
     346      kstable = ksta
    349347    ENDIF
    350 !IM cf JLD/ GKtest fin
    351 
    352 ! Calculer les geopotentiels de chaque couche
     348    !IM cf JLD/ GKtest fin
     349
     350    ! Calculer les geopotentiels de chaque couche
    353351
    354352    DO i = 1, knon
    355        zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
    356             * (paprs(i,1)-pplay(i,1))
     353      zgeop(i, 1) = RD * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) &
     354              * (paprs(i, 1) - pplay(i, 1))
    357355    ENDDO
    358356    DO k = 2, klev
    359        DO i = 1, knon
    360           zgeop(i,k) = zgeop(i,k-1) &
    361                + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
    362                * (pplay(i,k-1)-pplay(i,k))
    363        ENDDO
    364     ENDDO
    365 
    366 ! Calculer les coefficients turbulents dans l'atmosphere
     357      DO i = 1, knon
     358        zgeop(i, k) = zgeop(i, k - 1) &
     359                + RD * 0.5 * (t(i, k - 1) + t(i, k)) / paprs(i, k) &
     360                        * (pplay(i, k - 1) - pplay(i, k))
     361      ENDDO
     362    ENDDO
     363
     364    ! Calculer les coefficients turbulents dans l'atmosphere
    367365
    368366    DO i = 1, knon
    369        itop(i) = isommet
    370     ENDDO
    371 
     367      itop(i) = isommet
     368    ENDDO
    372369
    373370    DO k = 2, isommet
    374        DO i = 1, knon
    375           zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
    376                +(v(i,k)-v(i,k-1))**2)
    377           zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
    378           zdphi =zmgeom(i) / 2.0
    379           zt = (t(i,k)+t(i,k-1)) * 0.5
    380           zq = (q(i,k)+q(i,k-1)) * 0.5
    381 
    382 ! Calculer Qs et dQs/dT:
    383 
    384           IF (thermcep) THEN
    385              zdelta = MAX(0.,SIGN(1.,RTT-zt))
    386              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) &
    387                   + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
    388              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
    389              zqs = MIN(0.5,zqs)
    390              zcor = 1./(1.-RETV*zqs)
    391              zqs = zqs*zcor
    392              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
     371      DO i = 1, knon
     372        zdu2 = MAX(cepdu2, (u(i, k) - u(i, k - 1))**2 &
     373                + (v(i, k) - v(i, k - 1))**2)
     374        zmgeom(i) = zgeop(i, k) - zgeop(i, k - 1)
     375        zdphi = zmgeom(i) / 2.0
     376        zt = (t(i, k) + t(i, k - 1)) * 0.5
     377        zq = (q(i, k) + q(i, k - 1)) * 0.5
     378
     379        ! Calculer Qs et dQs/dT:
     380
     381        IF (thermcep) THEN
     382          zdelta = MAX(0., SIGN(1., RTT - zt))
     383          zcvm5 = R5LES * RLVTT / RCPD / (1.0 + RVTMP2 * zq) * (1. - zdelta) &
     384                  + R5IES * RLSTT / RCPD / (1.0 + RVTMP2 * zq) * zdelta
     385          zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
     386          zqs = MIN(0.5, zqs)
     387          zcor = 1. / (1. - RETV * zqs)
     388          zqs = zqs * zcor
     389          zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
     390        ELSE
     391          IF (zt < t_coup) THEN
     392            zqs = qsats(zt) / pplay(i, k)
     393            zdqs = dqsats(zt, zqs)
    393394          ELSE
    394              IF (zt < t_coup) THEN
    395                 zqs = qsats(zt) / pplay(i,k)
    396                 zdqs = dqsats(zt,zqs)
    397              ELSE
    398                 zqs = qsatl(zt) / pplay(i,k)
    399                 zdqs = dqsatl(zt,zqs)
    400              ENDIF
     395            zqs = qsatl(zt) / pplay(i, k)
     396            zdqs = dqsatl(zt, zqs)
    401397          ENDIF
    402 
    403 !           calculer la fraction nuageuse (processus humide):
    404 
    405           IF (zq /= 0.) THEN
    406              zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
    407           else
    408              zfr = 0.
    409           end if
    410           zfr = MAX(0.0,MIN(1.0,zfr))
    411           IF (.NOT.richum) zfr = 0.0
    412 
    413 !           calculer le nombre de Richardson:
    414 
    415           IF (tvirtu) THEN
    416              ztvd =( t(i,k) &
    417                   + zdphi/RCPD/(1.+RVTMP2*zq) &
    418                   *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
    419                   )*(1.+RETV*q(i,k))
    420              ztvu =( t(i,k-1) &
    421                   - zdphi/RCPD/(1.+RVTMP2*zq) &
    422                   *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
    423                   )*(1.+RETV*q(i,k-1))
    424              zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
    425              zri(i) = zri(i) &
    426                   + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
    427                   *(paprs(i,k)/101325.0)**RKAPPA &
    428                   /(zdu2*0.5*(ztvd+ztvu))
    429 
    430           ELSE ! calcul de Ridchardson compatible LMD5
    431 
    432              zri(i) =(RCPD*(t(i,k)-t(i,k-1)) &
    433                   -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) &
    434                   *(pplay(i,k)-pplay(i,k-1)) &
    435                   )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
    436              zri(i) = zri(i) + &
    437                   zmgeom(i)*zmgeom(i)*gamt(k)/RG &
    438                   *(paprs(i,k)/101325.0)**RKAPPA &
    439                   /(zdu2*0.5*(t(i,k-1)+t(i,k)))
     398        ENDIF
     399
     400        !           calculer la fraction nuageuse (processus humide):
     401
     402        IF (zq /= 0.) THEN
     403          zfr = (zq + ratqs * zq - zqs) / (2.0 * ratqs * zq)
     404        else
     405          zfr = 0.
     406        end if
     407        zfr = MAX(0.0, MIN(1.0, zfr))
     408        IF (.NOT.richum) zfr = 0.0
     409
     410        !           calculer le nombre de Richardson:
     411
     412        IF (tvirtu) THEN
     413          ztvd = (t(i, k) &
     414                  + zdphi / RCPD / (1. + RVTMP2 * zq) &
     415                          * ((1. - zfr) + zfr * (1. + RLVTT * zqs / RD / zt) / (1. + zdqs)) &
     416                  ) * (1. + RETV * q(i, k))
     417          ztvu = (t(i, k - 1) &
     418                  - zdphi / RCPD / (1. + RVTMP2 * zq) &
     419                          * ((1. - zfr) + zfr * (1. + RLVTT * zqs / RD / zt) / (1. + zdqs)) &
     420                  ) * (1. + RETV * q(i, k - 1))
     421          zri(i) = zmgeom(i) * (ztvd - ztvu) / (zdu2 * 0.5 * (ztvd + ztvu))
     422          zri(i) = zri(i) &
     423                  + zmgeom(i) * zmgeom(i) / RG * gamt(k) &
     424                          * (paprs(i, k) / 101325.0)**RKAPPA &
     425                          / (zdu2 * 0.5 * (ztvd + ztvu))
     426
     427        ELSE ! calcul de Ridchardson compatible LMD5
     428
     429          zri(i) = (RCPD * (t(i, k) - t(i, k - 1)) &
     430                  - RD * 0.5 * (t(i, k) + t(i, k - 1)) / paprs(i, k) &
     431                          * (pplay(i, k) - pplay(i, k - 1)) &
     432                  ) * zmgeom(i) / (zdu2 * 0.5 * RCPD * (t(i, k - 1) + t(i, k)))
     433          zri(i) = zri(i) + &
     434                  zmgeom(i) * zmgeom(i) * gamt(k) / RG &
     435                          * (paprs(i, k) / 101325.0)**RKAPPA &
     436                          / (zdu2 * 0.5 * (t(i, k - 1) + t(i, k)))
     437        ENDIF
     438
     439        !           finalement, les coefficients d'echange sont obtenus:
     440
     441        zcdn = SQRT(zdu2) / zmgeom(i) * RG
     442
     443        IF (opt_ec) THEN
     444          z2geomf = zgeop(i, k - 1) + zgeop(i, k)
     445          zalm2 = (0.5 * ckap / RG * z2geomf &
     446                  / (1. + 0.5 * ckap / rg / clam * z2geomf))**2
     447          zalh2 = (0.5 * ckap / rg * z2geomf &
     448                  / (1. + 0.5 * ckap / RG / (clam * SQRT(1.5 * cd)) * z2geomf))**2
     449          IF (zri(i)<0.0) THEN  ! situation instable
     450            zscf = ((zgeop(i, k) / zgeop(i, k - 1))**(1. / 3.) - 1.)**3 &
     451                    / (zmgeom(i) / RG)**3 / (zgeop(i, k - 1) / RG)
     452            zscf = SQRT(-zri(i) * zscf)
     453            zscfm = 1.0 / (1.0 + 3.0 * cb * cc * zalm2 * zscf)
     454            zscfh = 1.0 / (1.0 + 3.0 * cb * cc * zalh2 * zscf)
     455            pcfm(i, k) = zcdn * zalm2 * (1. - 2.0 * cb * zri(i) * zscfm)
     456            pcfh(i, k) = zcdn * zalh2 * (1. - 3.0 * cb * zri(i) * zscfh)
     457          ELSE ! situation stable
     458            zscf = SQRT(1. + cd * zri(i))
     459            pcfm(i, k) = zcdn * zalm2 / (1. + 2.0 * cb * zri(i) / zscf)
     460            pcfh(i, k) = zcdn * zalh2 / (1. + 3.0 * cb * zri(i) * zscf)
    440461          ENDIF
    441 
    442 !           finalement, les coefficients d'echange sont obtenus:
    443 
    444           zcdn=SQRT(zdu2) / zmgeom(i) * RG
    445 
    446           IF (opt_ec) THEN
    447              z2geomf=zgeop(i,k-1)+zgeop(i,k)
    448              zalm2=(0.5*ckap/RG*z2geomf &
    449                   /(1.+0.5*ckap/rg/clam*z2geomf))**2
    450              zalh2=(0.5*ckap/rg*z2geomf &
    451                   /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
    452              IF (zri(i)<0.0) THEN  ! situation instable
    453                 zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 &
    454                      / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
    455                 zscf = SQRT(-zri(i)*zscf)
    456                 zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
    457                 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
    458                 pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
    459                 pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
    460              ELSE ! situation stable
    461                 zscf=SQRT(1.+cd*zri(i))
    462                 pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
    463                 pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
    464              ENDIF
    465           ELSE
    466              zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) &
    467                   /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
    468              pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable))
    469              pcfm(i,k)= zl2(i)* pcfm(i,k)
    470              pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
    471           ENDIF
    472        ENDDO
    473     ENDDO
    474 
    475 ! Au-dela du sommet, pas de diffusion turbulente:
     462        ELSE
     463          zl2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, itop(i) + 1)) &
     464                  / (paprs(i, 2) - paprs(i, itop(i) + 1))))**2
     465          pcfm(i, k) = SQRT(MAX(zcdn * zcdn * (ric - zri(i)) / ric, kstable))
     466          pcfm(i, k) = zl2(i) * pcfm(i, k)
     467          pcfh(i, k) = pcfm(i, k) / prandtl ! h et m different
     468        ENDIF
     469      ENDDO
     470    ENDDO
     471
     472    ! Au-dela du sommet, pas de diffusion turbulente:
    476473
    477474    DO i = 1, knon
    478        IF (itop(i)+1 <= klev) THEN
    479           DO k = itop(i)+1, klev
    480              pcfh(i,k) = 0.0
    481              pcfm(i,k) = 0.0
    482           ENDDO
    483        ENDIF
    484     ENDDO
    485      
     475      IF (itop(i) + 1 <= klev) THEN
     476        DO k = itop(i) + 1, klev
     477          pcfh(i, k) = 0.0
     478          pcfm(i, k) = 0.0
     479        ENDDO
     480      ENDIF
     481    ENDDO
     482
    486483  END SUBROUTINE coefkz
    487484
    488 !****************************************************************************************
    489 
    490   SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
    491        pcfm, pcfh)
     485  !****************************************************************************************
     486
     487  SUBROUTINE coefkz2(nsrf, knon, paprs, pplay, t, &
     488          pcfm, pcfh)
    492489
    493490    USE dimphy
    494491    USE indice_sol_mod
    495492
    496 !======================================================================
    497 ! J'introduit un peu de diffusion sauf dans les endroits
    498 ! ou une forte inversion est presente
    499 ! On peut dire qu'il represente la convection peu profonde
    500 
    501 ! Arguments:
    502 ! nsrf-----input-I- indicateur de la nature du sol
    503 ! knon-----input-I- nombre de points a traiter
    504 ! paprs----input-R- pression a chaque intercouche (en Pa)
    505 ! pplay----input-R- pression au milieu de chaque couche (en Pa)
    506 ! t--------input-R- temperature (K)
    507 
    508 ! pcfm-----output-R- coefficients a calculer (vitesse)
    509 ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
    510 !======================================================================
    511 
    512 ! Arguments:
    513 
    514     INTEGER, INTENT(IN)                       :: knon, nsrf
    515     REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs
    516     REAL, DIMENSION(klon, klev), INTENT(IN)   :: pplay
    517     REAL, DIMENSION(klon, klev), INTENT(IN)   :: t(klon,klev)
    518    
    519     REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pcfm, pcfh
    520 
    521 ! Quelques constantes et options:
    522 
    523     REAL, PARAMETER :: prandtl=0.4
    524     REAL, PARAMETER :: kstable=0.002
    525 !   REAL, PARAMETER :: kstable=0.001
    526     REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange
    527     REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible
    528 !    PARAMETER (seuil=-0.04)
    529 !    PARAMETER (seuil=-0.06)
    530 !    PARAMETER (seuil=-0.09)
    531 
    532 ! Variables locales:
     493    !======================================================================
     494    ! J'introduit un peu de diffusion sauf dans les endroits
     495    ! ou une forte inversion est presente
     496    ! On peut dire qu'il represente la convection peu profonde
     497
     498    ! Arguments:
     499    ! nsrf-----input-I- indicateur de la nature du sol
     500    ! knon-----input-I- nombre de points a traiter
     501    ! paprs----input-R- pression a chaque intercouche (en Pa)
     502    ! pplay----input-R- pression au milieu de chaque couche (en Pa)
     503    ! t--------input-R- temperature (K)
     504
     505    ! pcfm-----output-R- coefficients a calculer (vitesse)
     506    ! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
     507    !======================================================================
     508
     509    ! Arguments:
     510
     511    INTEGER, INTENT(IN) :: knon, nsrf
     512    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
     513    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
     514    REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon, klev)
     515
     516    REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh
     517
     518    ! Quelques constantes et options:
     519
     520    REAL, PARAMETER :: prandtl = 0.4
     521    REAL, PARAMETER :: kstable = 0.002
     522    !   REAL, PARAMETER :: kstable=0.001
     523    REAL, PARAMETER :: mixlen = 35.0 ! constante controlant longueur de melange
     524    REAL, PARAMETER :: seuil = -0.02 ! au-dela l'inversion est consideree trop faible
     525    !    PARAMETER (seuil=-0.04)
     526    !    PARAMETER (seuil=-0.06)
     527    !    PARAMETER (seuil=-0.09)
     528
     529    ! Variables locales:
    533530
    534531    INTEGER i, k, invb(knon)
    … …  
    538535    INCLUDE "YOMCST.h"
    539536
    540 ! Initialiser les sorties
     537    ! Initialiser les sorties
    541538
    542539    DO k = 1, klev
    543        DO i = 1, knon
    544           pcfm(i,k) = 0.0
    545           pcfh(i,k) = 0.0
    546        ENDDO
    547     ENDDO
    548 
    549 ! Chercher la zone d'inversion forte
     540      DO i = 1, knon
     541        pcfm(i, k) = 0.0
     542        pcfh(i, k) = 0.0
     543      ENDDO
     544    ENDDO
     545
     546    ! Chercher la zone d'inversion forte
    550547
    551548    DO i = 1, knon
    552        invb(i) = klev
    553        zdthmin(i)=0.0
    554     ENDDO
    555     DO k = 2, klev/2-1
    556        DO i = 1, knon
    557           zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) &
    558                - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
    559           zdthdp = zdthdp * 100.0
    560           IF (pplay(i,k)>0.8*paprs(i,1) .AND. &
    561                zdthdp<zdthmin(i) ) THEN
    562              zdthmin(i) = zdthdp
    563              invb(i) = k
     549      invb(i) = klev
     550      zdthmin(i) = 0.0
     551    ENDDO
     552    DO k = 2, klev / 2 - 1
     553      DO i = 1, knon
     554        zdthdp = (t(i, k) - t(i, k + 1)) / (pplay(i, k) - pplay(i, k + 1)) &
     555                - RD * 0.5 * (t(i, k) + t(i, k + 1)) / RCPD / paprs(i, k + 1)
     556        zdthdp = zdthdp * 100.0
     557        IF (pplay(i, k)>0.8 * paprs(i, 1) .AND. &
     558                zdthdp<zdthmin(i)) THEN
     559          zdthmin(i) = zdthdp
     560          invb(i) = k
     561        ENDIF
     562      ENDDO
     563    ENDDO
     564
     565    ! Introduire une diffusion:
     566
     567    IF (nsrf==is_oce) THEN
     568      DO k = 2, klev
     569        DO i = 1, knon
     570          !IM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
     571          !IM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
     572          !IM cf JLD/ GKtest TERkz2
     573          ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
     574          ! fin GKtest
     575
     576
     577          ! s'il n'y a pas d'inversion ou si l'inversion est trop faible
     578          !          IF ( (nsrf.EQ.is_oce) .AND. &
     579          IF ((invb(i)==klev) .OR. (zdthmin(i)>seuil)) THEN
     580            zl2(i) = (mixlen * MAX(0.0, (paprs(i, k) - paprs(i, klev + 1)) &
     581                    / (paprs(i, 2) - paprs(i, klev + 1))))**2
     582            pcfm(i, k) = zl2(i) * kstable
     583            pcfh(i, k) = pcfm(i, k) / prandtl ! h et m different
    564584          ENDIF
    565        ENDDO
    566     ENDDO
    567 
    568 ! Introduire une diffusion:
    569 
    570     IF ( nsrf==is_oce ) THEN
    571        DO k = 2, klev
    572           DO i = 1, knon
    573 !IM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
    574 !IM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
    575       !IM cf JLD/ GKtest TERkz2
    576       ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
    577       ! fin GKtest
    578 
    579 
    580 ! s'il n'y a pas d'inversion ou si l'inversion est trop faible
    581 !          IF ( (nsrf.EQ.is_oce) .AND. &
    582              IF ( (invb(i)==klev) .OR. (zdthmin(i)>seuil) ) THEN
    583                 zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) &
    584                      /(paprs(i,2)-paprs(i,klev+1)) ))**2
    585                 pcfm(i,k)= zl2(i)* kstable
    586                 pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
    587              ENDIF
    588           ENDDO
    589        ENDDO
     585        ENDDO
     586      ENDDO
    590587    ENDIF
    591588
    592589  END SUBROUTINE coefkz2
    593590
    594 !****************************************************************************************
     591  !****************************************************************************************
    595592
    596593END MODULE coef_diff_turb_mod
Note: See TracChangeset for help on using the changeset viewer.