Ignore:
Timestamp:
Jun 10, 2016, 2:04:19 PM (8 years ago)
Author:
fhourdin
Message:

Parametrisation d'une longueur de melange verticale minimum associee
aux circulations meso-echelle introduites par le relief sous maille.
D'apres Etienne Vignon et Frédéric Hourdin

Location:
LMDZ5/trunk/libf/phylmd
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/clesphys.h

    r2524 r2561  
    4040!IM, MAFo fmagic, pmagic : parametres - additionnel et multiplicatif -
    4141!                          pour regler l albedo sur ocean
     42       REAL pbl_lmixmin_alpha
    4243       REAL fmagic, pmagic
    4344! Hauteur (imposee) du contenu en eau du sol
    … …  
    9495     &     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt                     &
    9596     &     , CH4_ppb_per, N2O_ppb_per, CFC11_ppt_per, CFC12_ppt_per     &
    96      &     , cdmmax, cdhmax, ksta, ksta_ter, f_ri_cd_min                &
     97     &     , cdmmax,cdhmax,ksta,ksta_ter,f_ri_cd_min,pbl_lmixmin_alpha  &
    9798     &     , fmagic, pmagic                                             &
    9899     &     , f_cdrag_ter,f_cdrag_oce,f_rugoro,z0min                     &

  • LMDZ5/trunk/libf/phylmd/coef_diff_turb_mod.F90

    r2311 r2561  
    164164               iflag_pbl)
    165165       ELSE IF (iflag_pbl<20) THEN
    166           CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
     166          CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &
    167167               yzlev,yzlay,yu,yv,yteta, &
    168168               ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &

  • LMDZ5/trunk/libf/phylmd/conf_phys_m.F90

    r2547 r2561  
    179179    REAL,SAVE :: cdmmax_omp,cdhmax_omp,ksta_omp,ksta_ter_omp,f_ri_cd_min_omp
    180180    LOGICAL,SAVE :: ok_kzmin_omp
     181    REAL, SAVE   :: pbl_lmixmin_alpha_omp
    181182    REAL, SAVE ::  fmagic_omp, pmagic_omp
    182183    INTEGER,SAVE :: iflag_pbl_omp,lev_histhf_omp,lev_histday_omp,lev_histmth_omp
    … …  
    12691270    ok_kzmin_omp = .true.
    12701271    call getin('ok_kzmin',ok_kzmin_omp)
     1272
     1273    pbl_lmixmin_alpha_omp=0.0
     1274    call getin('pbl_lmixmin_alpha',pbl_lmixmin_alpha_omp)
     1275
    12711276
    12721277    !
    … …  
    20662071    f_ri_cd_min = f_ri_cd_min_omp
    20672072    ok_kzmin = ok_kzmin_omp
     2073    pbl_lmixmin_alpha=pbl_lmixmin_alpha_omp
    20682074    fmagic = fmagic_omp
    20692075    pmagic = pmagic_omp
    … …  
    23832389    write(lunout,*)' f_ri_cd_min = ',f_ri_cd_min
    23842390    write(lunout,*)' ok_kzmin = ',ok_kzmin
     2391    write(lunout,*)' pbl_lmixmin_alpha = ',pbl_lmixmin_alpha
    23852392    write(lunout,*)' fmagic = ',fmagic
    23862393    write(lunout,*)' pmagic = ',pmagic

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

    r2536 r2561  
    1616      REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
    1717      !$OMP THREADPRIVATE(u_seri, v_seri)
     18      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:)
     19      !$OMP THREADPRIVATE(l_mixmin, l_mix)
    1820
    1921      REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
    … …  
    406408      allocate(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
    407409      allocate(u_seri(klon,klev),v_seri(klon,klev))
     410      allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbtr))
    408411
    409412      allocate(tr_seri(klon,klev,nbtr))
    … …  
    625628      deallocate(t_seri,q_seri,ql_seri,qs_seri)
    626629      deallocate(u_seri,v_seri)
     630      deallocate(l_mixmin,l_mix)
    627631
    628632      deallocate(tr_seri)

  • LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90

    r2553 r2561  
    745745  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_srf      = (/             &
    746746      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_ter',       &
    747       "Max Turb. Kinetic Energy "//clnsurf(1),"-", (/ ('', i=1, 9) /)), &
     747      "Max Turb. Kinetic Energy "//clnsurf(1),"m2/s2", (/ ('', i=1, 9) /)), &
    748748      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_lic',       &
    749       "Max Turb. Kinetic Energy "//clnsurf(2),"-", (/ ('', i=1, 9) /)), &
     749      "Max Turb. Kinetic Energy "//clnsurf(2),"m2/s2", (/ ('', i=1, 9) /)), &
    750750      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_oce',       &
    751       "Max Turb. Kinetic Energy "//clnsurf(3),"-", (/ ('', i=1, 9) /)), &
     751      "Max Turb. Kinetic Energy "//clnsurf(3),"m2/s2", (/ ('', i=1, 9) /)), &
    752752      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_sic',       &
    753       "Max Turb. Kinetic Energy "//clnsurf(4),"-", (/ ('', i=1, 9) /)) /)
     753      "Max Turb. Kinetic Energy "//clnsurf(4),"m2/s2", (/ ('', i=1, 9) /)) /)
     754
     755  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_l_mixmin      = (/             &
     756      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_ter',       &
     757      "PBL mixing length "//clnsurf(1),"m", (/ ('', i=1, 9) /)), &
     758      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_lic',       &
     759      "PBL mixing length "//clnsurf(2),"m", (/ ('', i=1, 9) /)), &
     760      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_oce',       &
     761      "PBL mixing length "//clnsurf(3),"m", (/ ('', i=1, 9) /)), &
     762      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_sic',       &
     763      "PBL mixing length "//clnsurf(4),"m", (/ ('', i=1, 9) /)) /)
     764
     765  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_l_mix      = (/             &
     766      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_ter',       &
     767      "min PBL mixing length "//clnsurf(1),"m", (/ ('', i=1, 9) /)), &
     768      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_lic',       &
     769      "min PBL mixing length "//clnsurf(2),"m", (/ ('', i=1, 9) /)), &
     770      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_oce',       &
     771      "min PBL mixing length "//clnsurf(3),"m", (/ ('', i=1, 9) /)), &
     772      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_sic',       &
     773      "min PBL mixing length "//clnsurf(4),"m", (/ ('', i=1, 9) /)) /)
    754774
    755775  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_max_srf  = (/                          &

  • LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90

    r2553 r2561  
    6060         o_fsw_srf, o_wbils_srf, o_wbilo_srf, &
    6161         o_tke_srf, o_tke_max_srf,o_dltpbltke_srf, o_wstar, &
     62         o_l_mixmin,o_l_mix, &
    6263         o_cdrm, o_cdrh, o_cldl, o_cldm, o_cldh, &
    6364         o_cldt, o_JrNt, o_cldljn, o_cldmjn, &
    … …  
    198199    USE phys_local_var_mod, only: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
    199200         t2m_min_mon, t2m_max_mon, evap, &
     201         l_mixmin,l_mix, &
    200202         zu10m, zv10m, zq2m, zustar, zxqsurf, &
    201203         rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
    … …  
    619621          IF (iflag_pbl > 1) THEN
    620622             CALL histwrite_phy(o_tke_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
     623             CALL histwrite_phy(o_l_mix(nsrf),  l_mix(:,1:klev,nsrf))
     624             CALL histwrite_phy(o_l_mixmin(nsrf),  l_mixmin(:,1:klev,nsrf))
    621625             CALL histwrite_phy(o_tke_max_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
    622626          ENDIF

  • LMDZ5/trunk/libf/phylmd/yamada4.F90

    r2441 r2561  
    1 
    2 ! $Header$
    3 
    4 SUBROUTINE yamada4(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
     1!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2
     3SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
    54    cd, q2, km, kn, kq, ustar, iflag_pbl)
     5
    66  USE dimphy
    7   USE print_control_mod, ONLY: prt_level
    8   USE ioipsl_getin_p_mod, ONLY : getin_p
    9 
    107  IMPLICIT NONE
    11 
     8  include "iniprint.h"
     9  ! .......................................................................
     10  ! ym#include "dimensions.h"
     11  ! ym#include "dimphy.h"
     12  ! ************************************************************************************************
     13  !
     14  ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
     15  !
     16  ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
     17  !            Yamada T.
     18  !            J. Atmos. Sci, 40, 91-106, 1983
     19  !
     20  !************************************************************************************************
     21  ! Input :
     22  !'======
     23  ! ni: indice horizontal sur la grille de base, non restreinte
     24  ! nsrf: type de surface
     25  ! ngrid: nombre de mailles concern??es sur l'horizontal
    1226  ! dt : pas de temps
    1327  ! g  : g
     28  ! rconst: constante de l'air sec
    1429  ! zlev : altitude a chaque niveau (interface inferieure de la couche
    1530  ! de meme indice)
    … …  
    1732  ! u,v : vitesse au centre de chaque couche
    1833  ! (en entree : la valeur au debut du pas de temps)
    19   ! teta : temperature potentielle au centre de chaque couche
     34  ! teta : temperature potentielle virtuelle au centre de chaque couche
    2035  ! (en entree : la valeur au debut du pas de temps)
    21   ! cd : cdrag
     36  ! cd : cdrag pour la quantit?? de mouvement
    2237  ! (en entree : la valeur au debut du pas de temps)
     38  ! ustar: vitesse de friction calcul??e par une formule diagnostique
     39  ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
     40  !
     41  !             iflag_pbl doit valoir entre 6 et 9
     42  !             l=6, on prend  systematiquement une longueur d'equilibre
     43  !             iflag_pbl=6 : MY 2.0
     44  !             iflag_pbl=7 : MY 2.0.Fournier
     45  !             iflag_pbl=8/9 : MY 2.5
     46  !             iflag_pbl=8 with special obsolete treatments for convergence
     47  !             with Cmpi5 NPv3.1 simulations
     48  !             iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
     49  !             iflag_pbl=12 = 11 with vertical diffusion off q2
     50  !
     51  !             2013/04/01 (FH hourdin@lmd.jussieu.fr)
     52  !             Correction for very stable PBLs (iflag_pbl=10 and 11)
     53  !             iflag_pbl=8 converges numerically with NPv3.1
     54  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
     55  !                          -> the model can run with longer time-steps.
     56  !
     57  ! Inpout/Output :
     58  !==============
    2359  ! q2 : $q^2$ au bas de chaque couche
    2460  ! (en entree : la valeur au debut du pas de temps)
    2561  ! (en sortie : la valeur a la fin du pas de temps)
     62 
     63  ! Outputs:
     64  !==========
    2665  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
    2766  ! couche)
    … …  
    2968  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
    3069  ! (en sortie : la valeur a la fin du pas de temps)
    31 
    32   ! iflag_pbl doit valoir entre 6 et 9
    33   ! l=6, on prend  systematiquement une longueur d'equilibre
    34   ! iflag_pbl=6 : MY 2.0
    35   ! iflag_pbl=7 : MY 2.0.Fournier
    36   ! iflag_pbl=8/9 : MY 2.5
    37   ! iflag_pbl=8 with special obsolete treatments for convergence
    38   ! with Cmpi5 NPv3.1 simulations
    39   ! iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
    40   ! iflag_pbl=12 = 11 with vertical diffusion off q2
    41 
    42   ! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
    43   ! Correction for very stable PBLs (iflag_pbl=10 and 11)
    44   ! iflag_pbl=8 converges numerically with NPv3.1
    45   ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
    46   ! -> the model can run with longer time-steps.
    47   ! .......................................................................
    48 
     70  !
     71  !.......................................................................
     72
     73  !=======================================================================
     74  ! Declarations:
     75  !=======================================================================
     76
     77
     78  ! Inputs/Outputs
     79  !----------------
    4980  REAL dt, g, rconst
    5081  REAL plev(klon, klev+1), temp(klon, klev)
    … …  
    5788  REAL teta(klon, klev)
    5889  REAL cd(klon)
    59   REAL q2(klon, klev+1), qpre
     90  REAL q2(klon, klev+1)
    6091  REAL unsdz(klon, klev)
    6192  REAL unsdzdec(klon, klev+1)
    62 
     93  REAL kn(klon, klev+1)
    6394  REAL km(klon, klev+1)
    64   REAL kmpre(klon, klev+1), tmp2
     95  INTEGER iflag_pbl, ngrid, nsrf
     96  INTEGER ni(klon)
     97
     98
     99  ! Local
     100  !-------
     101
     102  INCLUDE "clesphys.h"
     103
     104
     105  REAL kmpre(klon, klev+1), tmp2, qpre
    65106  REAL mpre(klon, klev+1)
    66   REAL kn(klon, klev+1)
    67107  REAL kq(klon, klev+1)
    68108  REAL ff(klon, klev+1), delta(klon, klev+1)
    69109  REAL aa(klon, klev+1), aa0, aa1
    70   INTEGER iflag_pbl, ngrid
    71110  INTEGER nlay, nlev
    72 
    73111  LOGICAL first
    74112  INTEGER ipas
    … …  
    77115  DATA first, ipas/.FALSE., 0/
    78116  !$OMP THREADPRIVATE( first,ipas)
    79   REAL,SAVE :: lmixmin=1.
    80   !$OMP THREADPRIVATE(lmixmin)
    81 
    82 
     117  LOGICAL hboville
    83118  INTEGER ig, k
    84 
    85 
    86119  REAL ri, zrif, zalpha, zsm, zsn
    87120  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
    88 
    89121  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
    90122  REAL dtetadz(klon, klev+1)
    91123  REAL m2cstat, mcstat, kmcstat
    92124  REAL l(klon, klev+1)
    93   REAL, ALLOCATABLE, SAVE :: l0(:)
    94   !$OMP THREADPRIVATE(l0)
    95   REAL sq(klon), sqz(klon), zz(klon, klev+1)
     125  REAL zz(klon, klev+1)
    96126  INTEGER iter
    97 
    98127  REAL ric, rifc, b1, kap
    99128  SAVE ric, rifc, b1, kap
    100129  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
    101130  !$OMP THREADPRIVATE(ric,rifc,b1,kap)
     131  REAL seuilsm, seuilalpha, lmixmin
    102132  REAL frif, falpha, fsm
    103   REAL fl, zzz, zl0, zq2, zn2
    104 
    105133  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
    106134    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
    107135  LOGICAL, SAVE :: firstcall = .TRUE.
    108136  !$OMP THREADPRIVATE(firstcall)
     137
     138
     139  ! Fonctions utilis??es
     140  !--------------------
     141
    109142  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    110143  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
    111144  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
    112   fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
    113     k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
    114     max(n2(ig,k),1.E-10))), lmixmin)
    115 
    116 
    117   nlay = klev
    118   nlev = klev + 1
    119 
    120   IF (firstcall) THEN
    121     ALLOCATE (l0(klon))
    122     firstcall = .FALSE.
    123     CALL getin_p('lmixmin',lmixmin)
    124   END IF
     145 
     146
     147!===============================================================================
     148! Flags, tests et ??valuations de constantes
     149!===============================================================================
     150
     151! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
     152   hboville=.TRUE.
    125153
    126154
    … …  
    129157  END IF
    130158
     159! Seuil dans le code de turbulence
     160 seuilalpha=1.12
     161 seuilsm=0.085
     162 lmixmin=1.
     163
     164  nlay = klev
     165  nlev = klev + 1
    131166  ipas = ipas + 1
    132167
    133168
    134   ! .......................................................................
    135   ! les increments verticaux
    136   ! .......................................................................
    137 
    138   ! !!!!! allerte !!!!!c
    139   ! !!!!! zlev n'est pas declare a nlev !!!!!c
    140   ! !!!!! ---->
     169!========================================================================
     170! Calcul des increments verticaux
     171!=========================================================================
     172
     173 
     174! Attention: zlev n'est pas declare a nlev
    141175  DO ig = 1, ngrid
    142176    zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
    143177  END DO
    144   ! !!!!! <----
    145   ! !!!!! allerte !!!!!c
     178
    146179
    147180  DO k = 1, nlay
    … …  
    162195  END DO
    163196
    164   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    165   ! Computing M^2, N^2, Richardson numbers, stability functions
    166   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    167     ! initialize arrays:
     197!=========================================================================
     198! Richardson number and stability functions
     199!=========================================================================
     200 
     201! initialize arrays:
     202
    168203  m2(:, :) = 0.0
    169204  sm(:, :) = 0.0
    170205  rif(:, :) = 0.0
    171206
     207!------------------------------------------------------------
    172208  DO k = 2, klev
     209
    173210    DO ig = 1, ngrid
    174211      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
    … …  
    177214      dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
    178215      n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
    179       ! n2(ig,k)=0.
    180216      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
    181217      IF (ri<ric) THEN
    … …  
    184220        rif(ig, k) = rifc
    185221      END IF
     222
    186223      IF (rif(ig,k)<0.16) THEN
    187224        alpha(ig, k) = falpha(rif(ig,k))
    188225        sm(ig, k) = fsm(rif(ig,k))
    189226      ELSE
    190         alpha(ig, k) = 1.12
    191         sm(ig, k) = 0.085
     227        alpha(ig, k) = seuilalpha
     228        sm(ig, k) = seuilsm
    192229      END IF
    193230      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
    … …  
    196233
    197234
    198   ! ====================================================================
    199   ! Computing the mixing length
    200   ! ====================================================================
    201 
    202   ! Mise a jour de l0
    203   IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
    204 
    205     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    206     ! Iterative computation of l0
    207     ! This version is kept for iflag_pbl only for convergence
    208     ! with NPv3.1 Cmip5 simulations
    209     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    210 
    211     DO ig = 1, ngrid
    212       sq(ig) = 1.E-10
    213       sqz(ig) = 1.E-10
    214     END DO
    215     DO k = 2, klev - 1
    216       DO ig = 1, ngrid
    217         zq = sqrt(q2(ig,k))
    218         sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
    219         sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
    220       END DO
    221     END DO
    222     DO ig = 1, ngrid
    223       l0(ig) = 0.2*sqz(ig)/sq(ig)
    224     END DO
     235
     236! ====================================================================
     237! Computing the mixing length
     238! ====================================================================
     239
     240 
     241  CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, firstcall, l)
     242
     243
     244
     245  !=======================================================================
     246  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
     247  !=======================================================================
     248
     249  !--------------
     250  ! Yamada 2.0
     251  !--------------
     252  IF (iflag_pbl==6) THEN
     253 
     254    DO k = 2, klev
     255      q2(:, k) = l(:, k)**2*zz(:, k)
     256    END DO
     257
     258
     259  !------------------
     260  ! Yamada 2.Fournier
     261  !------------------
     262
     263  ELSE IF (iflag_pbl==7) THEN
     264
     265
     266    ! Calcul de l,  km, au pas precedent
     267    !....................................
    225268    DO k = 2, klev
    226269      DO ig = 1, ngrid
    227         l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
    228       END DO
    229     END DO
    230     ! print*,'L0 cas 8 ou 10 ',l0
    231 
    232   ELSE
    233 
    234     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    235     ! In all other case, the assymptotic mixing length l0 is imposed (100m)
    236     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    237 
    238     l0(:) = 150.
    239     DO k = 2, klev
    240       DO ig = 1, ngrid
    241         l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
    242       END DO
    243     END DO
    244     ! print*,'L0 cas autres ',l0
    245 
    246   END IF
    247 
    248 
    249   ! ====================================================================
    250   ! Yamada 2.0
    251   ! ====================================================================
    252   IF (iflag_pbl==6) THEN
    253 
    254     DO k = 2, klev
    255       q2(:, k) = l(:, k)**2*zz(:, k)
    256     END DO
    257 
    258 
    259   ELSE IF (iflag_pbl==7) THEN
    260     ! ====================================================================
    261     ! Yamada 2.Fournier
    262     ! ====================================================================
    263 
    264     ! Calcul de l,  km, au pas precedent
    265     DO k = 2, klev
    266       DO ig = 1, ngrid
    267         ! print*,'SMML=',sm(ig,k),l(ig,k)
    268270        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
    269271        kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
    270272        mpre(ig, k) = sqrt(m2(ig,k))
    271         ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    272273      END DO
    273274    END DO
    … …  
    278279        mcstat = sqrt(m2cstat)
    279280
    280         ! print*,'M2 L=',k,mpre(ig,k),mcstat
    281 
    282         ! -----{puis on ecrit la valeur de q qui annule l'equation de m
    283         ! supposee en q3}
     281     ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
     282     !.........................................................................
    284283
    285284        IF (k==2) THEN
    … …  
    293292            (unsdz(ig,k)+unsdz(ig,k-1))
    294293        END IF
    295         ! print*,'T2 L=',k,tmp2
     294
    296295        tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
    297296        q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
    298         ! print*,'Q2 L=',k,q2(ig,k)
    299297
    300298      END DO
    301299    END DO
    302300
     301
     302    ! ------------------------
     303    ! Yamada 2.5 a la Didi
     304    !-------------------------
     305
    303306  ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
    304     ! ====================================================================
    305     ! Yamada 2.5 a la Didi
    306     ! ====================================================================
    307 
    308 
    309     ! Calcul de l,  km, au pas precedent
     307
     308    ! Calcul de l, km, au pas precedent
     309    !....................................
    310310    DO k = 2, klev
    311311      DO ig = 1, ngrid
    312         ! print*,'SMML=',sm(ig,k),l(ig,k)
    313312        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
    314313        IF (delta(ig,k)<1.E-20) THEN
    315           ! print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
    316314          delta(ig, k) = 1.E-20
    317315        END IF
    … …  
    319317        aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
    320318        aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
    321         ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
    322319        aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
    323         ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    324320        qpre = sqrt(q2(ig,k))
    325         ! if (iflag_pbl.eq.8 ) then
    326321        IF (aa(ig,k)>0.) THEN
    327322          q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
    … …  
    337332        ! endif
    338333        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
    339         ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    340334      END DO
    341335    END DO
    … …  
    343337  ELSE IF (iflag_pbl>=10) THEN
    344338
    345     ! print*,'Schema mixte D'
    346     ! print*,'Longueur ',l(:,:)
    347339    DO k = 2, klev - 1
    348       DO ig = 1, ngrid
    349         l(ig, k) = max(l(ig,k), lmixmin)
    350         km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
    351         q2(ig, k) = q2(ig, k) + dt*km(ig, k)*m2(ig, k)*(1.-rif(ig,k))
    352         q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
    353         q2(ig, k) = 1./(1./sqrt(q2(ig,k))+dt/(2*l(ig,k)*b1))
    354         q2(ig, k) = q2(ig, k)*q2(ig, k)
    355       END DO
     340      km(:, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k)
     341      q2(:, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k))
     342      q2(:, k) = min(max(q2(:,k),1.E-10), 1.E4)
     343      q2(:, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))
     344      q2(:, k) = q2(:, k)*q2(:, k)
    356345    END DO
    357346
    … …  
    362351  END IF ! Fin du cas 8
    363352
    364   ! print*,'OK8'
    365353
    366354  ! ====================================================================
    367   ! Calcul des coefficients de m�ange
     355  ! Calcul des coefficients de melange
    368356  ! ====================================================================
     357
    369358  DO k = 2, klev
    370     ! print*,'k=',k
    371359    DO ig = 1, ngrid
    372       ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
    373360      zq = sqrt(q2(ig,k))
    374       km(ig, k) = l(ig, k)*zq*sm(ig, k)
    375       kn(ig, k) = km(ig, k)*alpha(ig, k)
    376       kq(ig, k) = l(ig, k)*zq*0.2
    377       ! print*,'KML=',km(ig,k),kn(ig,k)
    378     END DO
    379   END DO
     361      km(ig, k) = l(ig, k)*zq*sm(ig, k)     ! For momentum
     362      kn(ig, k) = km(ig, k)*alpha(ig, k)    ! For scalars
     363      kq(ig, k) = l(ig, k)*zq*0.2           ! For TKE
     364    END DO
     365  END DO
     366
     367
     368  !====================================================================
     369  ! Transport diffusif vertical de la TKE par la TKE
     370  !====================================================================
     371
     372
    380373    ! initialize near-surface and top-layer mixing coefficients
    381   kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
    382   kq(1:ngrid, klev+1) = 0 ! zero at the top
    383 
    384   ! Transport diffusif vertical de la TKE.
     374    !...........................................................
     375
     376  kq(1:ngrid, 1) = kq(1:ngrid, 2)    ! constant (ie no gradient) near the surface
     377  kq(1:ngrid, klev+1) = 0            ! zero at the top
     378
     379    ! Transport diffusif vertical de la TKE.
     380    !.......................................
     381
    385382  IF (iflag_pbl>=12) THEN
    386     ! print*,'YAMADA VDIF'
    387383    q2(:, 1) = q2(:, 2)
    388384    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
    389385  END IF
    390386
    391   ! Traitement des cas noctrunes avec l'introduction d'une longueur
    392   ! minilale.
    393 
    394   ! ====================================================================
    395   ! Traitement particulier pour les cas tres stables.
    396   ! D'apres Holtslag Boville.
     387
     388  !====================================================================
     389  ! Traitement particulier pour les cas tres stables, introduction d'une
     390  ! longueur de m??lange minimale
     391  !====================================================================
     392  !
     393  ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
     394  !            Holtslag A.A.M. and Boville B.A.
     395  !            J. Clim., 6, 1825-1842, 1993
     396
     397
     398 IF (hboville) THEN
     399
    397400
    398401  IF (prt_level>1) THEN
    399402    PRINT *, 'YAMADA4 0'
    400   END IF !(prt_level>1) THEN
     403  END IF
     404
    401405  DO ig = 1, ngrid
    402406    coriol(ig) = 1.E-4
    … …  
    404408  END DO
    405409
    406   ! print*,'pblhmin ',pblhmin
    407   ! Test a remettre 21 11 02
    408   ! test abd 13 05 02      if(0.eq.1) then
    409410  IF (1==1) THEN
    410411    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
    … …  
    419420          END IF
    420421          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
    421             ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    422             ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
     422
    423423            kn(ig, k) = kmin
    424424            km(ig, k) = kmin
    425425            kq(ig, k) = kmin
    426             ! la longueur de melange est suposee etre l= kap z
    427             ! K=l q Sm d'ou q2=(K/l Sm)**2
     426
     427 ! la longueur de melange est suposee etre l= kap z
     428 ! K=l q Sm d'ou q2=(K/l Sm)**2
     429
    428430            q2(ig, k) = (qmin/sm(ig,k))**2
    429431          END IF
    … …  
    432434
    433435    ELSE
    434 
    435436      DO k = 2, klev
    436437        DO ig = 1, ngrid
    … …  
    442443          END IF
    443444          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
    444             ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    445             ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    446445            kn(ig, k) = kmin
    447446            km(ig, k) = kmin
    448447            kq(ig, k) = kmin
    449             ! la longueur de melange est suposee etre l= kap z
    450             ! K=l q Sm d'ou q2=(K/l Sm)**2
     448 ! la longueur de melange est suposee etre l= kap z
     449 ! K=l q Sm d'ou q2=(K/l Sm)**2
    451450            sm(ig, k) = 1.
    452451            alpha(ig, k) = 1.
    … …  
    463462  END IF
    464463
     464 END IF ! hboville
     465
    465466  IF (prt_level>1) THEN
    466467    PRINT *, 'YAMADA4 1'
    467468  END IF !(prt_level>1) THEN
     469
     470
     471 !======================================================
     472 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
     473 !======================================================
     474 !
     475 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
     476 !            Abdella K and McFarlane N
     477 !            J. Atmos. Sci., 54, 1850-1867, 1997
     478
    468479  ! Diagnostique pour stokage
     480  !..........................
    469481
    470482  IF (1==0) THEN
    … …  
    482494    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
    483495
    484     ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
     496
     497  ! Calcul de w'2 et T'2
     498  !.......................
    485499
    486500    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
    … …  
    493507  END IF
    494508
    495   ! print*,'OKFIN'
     509!============================================================================
     510
    496511  first = .FALSE.
    497512  RETURN
     513
     514
    498515END SUBROUTINE yamada4
     516
     517!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     518
     519!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    499520SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
     521
    500522  USE dimphy
    501523  IMPLICIT NONE
    502 
    503   ! dt : pas de temps
     524 
     525  include "dimensions.h"
     526
     527!    vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
     528!             avec un schema implicite en temps avec
     529!             inversion d'un syst??me tridiagonal
     530!
     531!     Reference: Description of the interface with the surface and
     532!                the computation of the turbulet diffusion in LMDZ
     533!                Technical note on LMDZ
     534!                Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
     535!
     536!============================================================================
     537! Declarations
     538!============================================================================
    504539
    505540  REAL plev(klon, klev+1)
    … …  
    516551  INTEGER i, k
    517552
    518   ! print*,'RD=',rconst
     553
     554!=========================================================================
     555! Calcul
     556!=========================================================================
     557
    519558  DO k = 1, klev
    520559    DO i = 1, ngrid
    521       ! test
    522       ! print*,'i,k',i,k
    523       ! print*,'temp(i,k)=',temp(i,k)
    524       ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
    525560      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
    526561      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
    … …  
    546581      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
    547582        kstar(i, k-1)
    548       ! correction d'un bug 10 01 2001
    549583      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
    550584      beta(i, k) = kstar(i, k-1)/denom(i, k)
    … …  
    553587
    554588  ! Si on recalcule q2(1)
     589  !.......................
    555590  IF (1==0) THEN
    556591    DO i = 1, ngrid
    … …  
    559594    END DO
    560595  END IF
    561   ! sinon, on peut sauter cette boucle...
     596
    562597
    563598  DO k = 2, klev + 1
    … …  
    569604  RETURN
    570605END SUBROUTINE vdif_q2
    571 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
    572   USE dimphy
     606!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     607
     608
     609
     610!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     611 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
     612 
     613   USE dimphy
    573614  IMPLICIT NONE
    574615
    575   ! dt : pas de temps
     616  include "dimensions.h"
     617!
     618! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
     619!           avec un schema explicite en temps
     620
     621
     622!====================================================
     623! Declarations
     624!====================================================
    576625
    577626  REAL plev(klon, klev+1)
    … …  
    585634  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
    586635  INTEGER ngrid
    587 
    588636  INTEGER i, k
     637
     638
     639!==================================================
     640! Calcul
     641!==================================================
    589642
    590643  DO k = 1, klev
    … …  
    621674  RETURN
    622675END SUBROUTINE vdif_q2e
     676
     677!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     678
     679
     680!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     681
     682SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, firstcall, lmix)
     683
     684
     685
     686  USE dimphy
     687  USE phys_state_var_mod, only: zstd, zsig, zmea
     688  USE phys_local_var_mod, only: l_mixmin, l_mix
     689
     690 ! zstd: ecart type de la'altitud e sous-maille
     691 ! zmea: altitude moyenne sous maille
     692 ! zsig: pente moyenne de le maille
     693
     694  USE geometry_mod, only: cell_area
     695  ! aire_cell: aire de la maille
     696
     697  IMPLICIT NONE
     698!*************************************************************************
     699! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
     700! avec la formule de Blackadar
     701! Calcul d'un  minimum en fonction de l'orographie sous-maille:
     702! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
     703! induit par les circulations meso et submeso ??chelles.
     704!
     705! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
     706!               Blackadar A.K.
     707!               J. Geophys. Res., 64, No 8, 1962
     708!
     709!             * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
     710!               to large eddy simulations
     711!               Ayotte K et al
     712!               Boundary Layer Meteorology, 79, 131-175, 1996
     713!
     714!
     715!             * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
     716!               Van de Wiel B.J.H et al
     717!               Boundary-Lay Meteorol, 128, 103-166, 2008
     718!
     719!
     720! Histoire:
     721!----------
     722! * premi??re r??daction, Etienne et Frederic, 09/06/2016
     723!
     724! ***********************************************************************
     725
     726!==================================================================
     727! Declarations
     728!==================================================================
     729
     730! Inputs
     731!-------
     732 INTEGER            ni(klon)           ! indice sur la grille original (non restreinte)
     733 INTEGER            nsrf               ! Type de surface
     734 INTEGER            ngrid              ! Nombre de points concern??s sur l'horizontal
     735 INTEGER            iflag_pbl          ! Choix du sch??ma de turbulence
     736 REAL            pbl_lmixmin_alpha  ! on active ou non le calcul de la longueur de melange minimum
     737 REAL               lmixmin            ! Minimum absolu de la longueur de m??lange
     738 REAL               zlay(klon, klev)   ! altitude du centre de la couche
     739 REAL               zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
     740 REAL               u(klon, klev)      ! vitesse du vent zonal
     741 REAL               v(klon, klev)      ! vitesse du vent meridional
     742 REAL               q2(klon, klev+1)   ! energie cin??tique turbulente
     743 REAL               n2(klon, klev+1)   ! frequence de Brunt-Vaisala
     744
     745!In/out
     746!-------
     747
     748 LOGICAL            firstcall          ! premier appel au code
     749
     750! Outputs
     751!---------
     752
     753 REAL               lmix(klon, klev+1)    ! Longueur de melange 
     754
     755
     756! Local
     757!-------
     758 
     759 INTEGER  ig,jg, k
     760 REAL     h_oro(klon)
     761 REAL     hlim(klon)
     762 REAL, SAVE :: kap=0.4,kapb=0.4
     763 REAL zq
     764 REAL sq(klon), sqz(klon)
     765 REAL, ALLOCATABLE, SAVE :: l0(:)
     766  !$OMP THREADPRIVATE(l0)
     767 REAL fl, zzz, zl0, zq2, zn2
     768 REAL famorti, zzzz, zh_oro, zhlim
     769 REAL l1(klon, klev+1), l2(klon,klev+1)
     770 REAL winds(klon, klev)
     771 REAL xcell
     772 REAL zstdslope(klon) 
     773 REAL lmax
     774 REAL l2strat, l2neutre, extent 
     775 REAL l2limit(klon)
     776!===============================================================
     777! Fonctions utiles
     778!===============================================================
     779
     780! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
     781!..........................................................................
     782
     783 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
     784    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
     785    max(n2(ig,k),1.E-10))), 1.E-5)
     786 
     787! Fonction d'amortissement de la turbulence au dessus de la montagne
     788! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
     789!.....................................................................
     790
     791 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1.   
     792
     793
     794 
     795
     796
     797  IF (firstcall) THEN
     798    ALLOCATE (l0(klon))
     799    firstcall = .FALSE.
     800  END IF
     801
     802
     803
     804!=====================================================================
     805!         CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
     806!=====================================================================
     807
     808
     809  IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
     810
     811   
     812    ! Iterative computation of l0
     813    ! This version is kept for iflag_pbl only for convergence
     814    ! with NPv3.1 Cmip5 simulations
     815    !...................................................................
     816
     817    DO ig = 1, ngrid
     818      sq(ig) = 1.E-10
     819      sqz(ig) = 1.E-10
     820    END DO
     821    DO k = 2, klev - 1
     822      DO ig = 1, ngrid
     823        zq = sqrt(q2(ig,k))
     824        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
     825        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
     826      END DO
     827    END DO
     828    DO ig = 1, ngrid
     829      l0(ig) = 0.2*sqz(ig)/sq(ig)
     830    END DO
     831    DO k = 2, klev
     832      DO ig = 1, ngrid
     833        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
     834      END DO
     835    END DO
     836
     837  ELSE
     838
     839   
     840    ! In all other case, the assymptotic mixing length l0 is imposed (150m)
     841    !......................................................................
     842
     843    l0(:) = 150.
     844    DO k = 2, klev
     845      DO ig = 1, ngrid
     846        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
     847      END DO
     848    END DO
     849
     850  END IF
     851
     852!=================================================================================
     853!  CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
     854! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
     855! glacier, la glace de mer et les oc??ans)
     856!=================================================================================
     857
     858   l2(:,:)=0.0
     859   l_mixmin(:,:,nsrf)=0.
     860   l_mix(:,:,nsrf)=0.
     861
     862   IF (nsrf .EQ. 1) THEN
     863
     864! coefficients
     865!--------------
     866
     867  extent=2.                                                         ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. 
     868  lmax=150.                                                         ! Longueur de m??lange max dans l'absolu
     869
     870! calculs
     871!---------
     872
     873  DO ig=1,ngrid
     874
     875      ! On calcule la hauteur du relief
     876      !.................................
     877
     878      ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
     879      ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
     880      ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
     881      ! (en gros, une maille de taille xcell avec une pente constante zstdslope)
     882 
     883
     884      jg=ni(ig)
     885 
     886!     IF (zsig(jg) .EQ. 0.) THEN
     887!          zstdslope(ig)=0.         
     888!     ELSE
     889!     xcell=sqrt(cell_area(jg))
     890!     zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.)
     891!     zstdslope(ig)=sqrt(zstdslope(ig))
     892!     END IF
     893     
     894!     h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.)   ! Hauteur du relief
     895      h_oro(ig)=zstd(jg)
     896      hlim(ig)=extent*h_oro(ig)     
     897
     898     
     899
     900
     901   END DO
     902
     903  l2limit(1:ngrid)=0.
     904
     905  DO k=2,klev
     906   DO ig=1,ngrid
     907 
     908      winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
     909
     910      IF (zlev(ig,k) .LE. h_oro(ig)) THEN                                     ! Si on est sous l'orographie
     911
     912
     913      l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10))              ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
     914                                               
     915      l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig))            ! Dans le cas neutre, formule de blackadar qui tend asymptotiquement vers h
     916      l2neutre=MIN(l2neutre,lmax)                                             ! On majore par lmax
     917      l2limit(ig)=MIN(l2neutre,l2strat)                                       ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
     918     
     919      l2(ig,k)=l2limit(ig)
     920                                     
     921      ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
     922
     923      ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
     924      ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
     925      ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
     926
     927         l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
     928
     929     
     930      ELSE                                                                    ! Au dessus de extent*h, on prend l2=l0
     931
     932      l2(ig,k)=0.
     933
     934
     935      END IF
     936       
     937    END DO
     938  END DO
     939
     940
     941 END IF                                                                        ! pbl_lmixmin_alpha
     942
     943!==================================================================================
     944! On prend le max entre la longueur de melange de blackadar et celle calcul??e
     945! en fonction de la topographie
     946!===================================================================================
     947
     948
     949 DO k=2,klev
     950    DO ig=1,ngrid
     951
     952   lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
     953
     954   END DO
     955 END DO
     956
     957 DO k=2,klev
     958  DO ig=1,ngrid
     959   jg=ni(ig)
     960   l_mix(jg,k,nsrf)=lmix(ig,k)
     961   l_mixmin(jg,k,nsrf)=l2(ig,k)
     962  ENDDO
     963 ENDDO
     964  DO ig=1,ngrid
     965   jg=ni(ig)
     966      l_mix(jg,1,nsrf)=hlim(ig)
     967  ENDDO
     968
     969
     970
     971END SUBROUTINE mixinglength
     972
     973!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Note: See TracChangeset for help on using the changeset viewer.