Changeset 2680


Ignore:
Timestamp:
Oct 21, 2016, 7:58:08 PM (8 years ago)
Author:
fhourdin
Message:

Rustine sur ocean_albedo pour le 1D
Correction sur la version conservative de Mellor et Yamada

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

Legend:

Unmodified
Added
Removed

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

    r2677 r2680  
    8080!--initialisations-------------
    8181!
     82
     83IF (knon==0) RETURN ! A verifier pourquoi on en a besoin...
     84
    8285alb_dir_new(:,:) = 0.
    8386alb_dif_new(:,:) = 0.
    … …  
    182185  ZUE=0.676               ! equivalent u_unif for diffuse incidence
    183186  ZUE2=SQRT(1.0-(1.0-ZUE**2)/ZREFM**2)
     187print*,'knon',knon
     188stop
    184189  ZRR0(1:knon)=0.50*(((ZUE2-ZREFM*ZUE)/(ZUE2+ZREFM*ZUE))**2 +((ZUE-ZREFM*ZUE2)/(ZUE+ZREFM*ZUE2))**2)
    185190  ZRRR(1:knon)=0.50*(((ZUE2-1.34*ZUE)/(ZUE2+1.34*ZUE))**2 +((ZUE-1.34*ZUE2)/(ZUE+1.34*ZUE2))**2)

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

    r2670 r2680  
    849849 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
    850850 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
    851  zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
    852851 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
    853852 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
    … …  
    21032102        CALL yamada_c(knon,dtime,ypaprs,ypplay &
    21042103    &   ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
    2105     &   ,iflag_pbl,nsrf)
     2104    &   ,iflag_pbl)
    21062105     ENDIF
    21072106!     print*,'yamada_c OK'
    … …  
    21212120    &   ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &
    21222121        ,ycoefq_x,y_d_t_diss_x,yustar_x &
    2123     &   ,iflag_pbl,nsrf)
     2122    &   ,iflag_pbl)
    21242123     ENDIF
    21252124!     print*,'yamada_c OK'
    … …  
    21382137    &   ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &
    21392138        ,ycoefq_w,y_d_t_diss_w,yustar_w &
    2140     &   ,iflag_pbl,nsrf)
     2139    &   ,iflag_pbl)
    21412140     ENDIF
    21422141!     print*,'yamada_c OK'
    … …  
    23262325!!!
    23272326       IF (iflag_split .eq.0) THEN
    2328         DO k = 2, klev
     2327        DO k = 1, klev
    23292328           DO j = 1, knon
    23302329              i = ni(j)
    … …  
    23392338
    23402339       ELSE
    2341         DO k = 2, klev
     2340        DO k = 1, klev
    23422341          DO j = 1, knon
    23432342            i = ni(j)

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

    r2391 r2680  
    44      SUBROUTINE yamada_c(ngrid,timestep,plev,play &
    55     &   ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar &
    6      &   ,iflag_pbl,okiophys)
     6     &   ,iflag_pbl)
    77      USE dimphy, ONLY: klon, klev
    88      USE print_control_mod, ONLY: prt_level
     9      USE ioipsl_getin_p_mod, ONLY : getin_p
     10
    911      IMPLICIT NONE
    1012#include "YOMCST.h"
    … …  
    5052      REAL, DIMENSION(klon,klev) :: pu,pv,pt
    5153      REAL, DIMENSION(klon,klev) :: d_t_diss
    52       INTEGER okiophys
    5354
    5455      REAL timestep
    … …  
    6869      REAL unsdzdec(klon,klev+1)
    6970
    70       REAL km(klon,klev+1)
     71      REAL km(klon,klev)
    7172      REAL kmpre(klon,klev+1),tmp2
    7273      REAL mpre(klon,klev+1)
    73       REAL kn(klon,klev+1)
    74       REAL kq(klon,klev+1)
     74      REAL kn(klon,klev)
     75      REAL kq(klon,klev)
    7576      real ff(klon,klev+1),delta(klon,klev+1)
    7677      real aa(klon,klev+1),aa0,aa1
    … …  
    8485      data first,ipas/.false.,0/
    8586!$OMP THREADPRIVATE( first,ipas)
     87       INTEGER, SAVE :: iflag_tke_diff=0
     88!$OMP THREADPRIVATE(iflag_tke_diff)
     89
    8690
    8791      integer ig,k
    … …  
    119123REAL, DIMENSION(klon,klev) :: exner,masse
    120124REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
     125      LOGICAL okiophys
    121126
    122127      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    … …  
    128133
    129134
     135      okiophys=klon==1
    130136      if (firstcall) then
     137        CALL getin_p('iflag_tke_diff',iflag_tke_diff)
    131138        allocate(l0(klon))
    132 #ifdef IOPHYS
    133         call iophys_ini
     139#define IOPHYS
     140#ifdef IOPHYS
     141!        call iophys_ini
    134142#endif
    135143        firstcall=.false.
    136144      endif
    137145
    138 
    139 #ifdef IOPHYS
    140 if (okiophys==1) then
     146   IF (ngrid<=0) RETURN ! Bizarre : on n a pas ce probeleme pour coef_diff_turb
     147
     148#ifdef IOPHYS
     149if (okiophys) then
    141150call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
    142151call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
    … …  
    146155      nlay=klev
    147156      nlev=klev+1
     157
    148158
    149159!-------------------------------------------------------------------------
    … …  
    152162
    153163
    154    zu(:,:)=pu(:,:)+0.5*d_u(:,:)
    155    zv(:,:)=pv(:,:)+0.5*d_v(:,:)
    156    zt(:,:)=pt(:,:)+0.5*d_t(:,:)
     164   zalpha=0.5 ! Anciennement 0.5. Essayer de voir pourquoi ?
     165   zu(:,:)=pu(:,:)+zalpha*d_u(:,:)
     166   zv(:,:)=pv(:,:)+zalpha*d_v(:,:)
     167   zt(:,:)=pt(:,:)+zalpha*d_t(:,:)
     168
    157169   do k=1,klev
    158170      exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA
    159171      masse(:,k)=(plev(:,k)-plev(:,k+1))/RG
     172      teta(:,k)=zt(:,k)/exner(:,k)
    160173   enddo
    161    teta(:,:)=zt(:,:)/exner(:,:)
    162174
    163175! Atmospheric mass at layer interfaces, where the TKE is computed
    … …  
    168180    enddo
    169181    masseb(:,:)=0.5*masseb(:,:)
    170 
    171 
    172182
    173183   zlev(:,1)=0.
    … …  
    202212
    203213#ifdef IOPHYS
    204 if (okiophys==1) then
     214if (okiophys) then
    205215      call iophys_ecrit('zlay',klev,'Geop','m',zlay)
    206216      call iophys_ecrit('teta',klev,'teta','K',teta)
    207217      call iophys_ecrit('temp',klev,'temp','K',zt)
    208218      call iophys_ecrit('pt',klev,'temp','K',pt)
     219      call iophys_ecrit('pu',klev,'u','m/s',pu)
     220      call iophys_ecrit('pv',klev,'v','m/s',pv)
    209221      call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u)
    210222      call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v)
    … …  
    213225      call iophys_ecrit('masse',klev,'masse','',masse)
    214226      call iophys_ecrit('masseb',klev,'masseb','',masseb)
    215       call iophys_ecrit('Cm2',klev,'m2 conserv','m/s',(dddu(:,1:klev)+dddv(:,1:klev))/(masseb(:,1:klev)*timestep))
    216       call iophys_ecrit('Cn2',klev,'m2 conserv','m/s',(rcpd*dddt(:,1:klev)/masseb(:,1:klev))/timestep)
    217       call iophys_ecrit('rifc',klev,'rif conservative','',rcpd*dddt(:,1:klev)/min(dddu(:,1:klev)+dddv(:,1:klev),-1.e-20))
    218227endif
    219228#endif
    … …  
    273282            rif(ig,k)=rifc
    274283         endif
    275          if(rif(ig,k).lt.0.16) then
     284         if(rif(ig,k)<0.16) then
    276285            alpha(ig,k)=falpha(rif(ig,k))
    277286            sm(ig,k)=fsm(rif(ig,k))
    … …  
    338347
    339348#ifdef IOPHYS
    340 if (okiophys==1) then
     349if (okiophys) then
    341350call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
    342351call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
    343 call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
     352call iophys_ecrit('Km2app',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
    344353call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
    345354endif
    … …  
    357366! Evolution of TKE under source terms K M2 and K N2
    358367   leff(:,:)=max(l(:,:),1.)
    359    IF (iflag_pbl==29) THEN
    360       km2(:,:)=km(:,:)*m2(:,:)
    361       kn2(:,:)=kn2(:,:)*rif(:,:)
    362    ELSEIF (iflag_pbl==25) THEN
    363       DO k=1,klev
    364          km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
    365          &        /(masse(:,k)*timestep)
    366          kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
    367          leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
    368       ENDDO
    369       km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
    370    ELSE
     368
     369!##################################################################
     370!#  IF (iflag_pbl==29) THEN
     371!#     STOP'Ne pas utiliser iflag_pbl=29'
     372!#     km2(:,:)=km(:,:)*m2(:,:)
     373!#     kn2(:,:)=kn2(:,:)*rif(:,:)
     374!#  ELSEIF (iflag_pbl==25) THEN
     375! VERSION AVEC LA TKE EN MILIEU DE COUCHE
     376!#     STOP'Ne pas utiliser iflag_pbl=25'
     377!#     DO k=1,klev
     378!#        km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
     379!#        &        /(masse(:,k)*timestep)
     380!#        kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
     381!#        leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
     382!#     ENDDO
     383!#     km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
     384!#  ELSE
     385!#################################################################
     386
    371387      km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep)
    372388      kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep)
    373    ENDIF
     389!   ENDIF
    374390   q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:))
    375391   q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4)
     392
     393 
     394#ifdef IOPHYS
     395if (okiophys) then
     396      call iophys_ecrit('km2',klev,'m2 conserv','m/s',km2(:,1:klev))
     397      call iophys_ecrit('kn2',klev,'n2 conserv','m/s',kn2(:,1:klev))
     398endif
     399#endif
    376400
    377401! Dissipation of TKE
    … …  
    379403   q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1))
    380404   q2(:,:)=q2(:,:)*q2(:,:)
    381    IF (iflag_pbl<=24) THEN
     405!  IF (iflag_pbl<=24) THEN
    382406      DO k=1,klev
    383407         d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k))
    384408      ENDDO
    385    ELSE IF (iflag_pbl<=27) THEN
    386       DO k=1,klev
    387          d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
    388       ENDDO
    389    ENDIF
     409
     410!###################################################################
     411!  ELSE IF (iflag_pbl<=27) THEN
     412!     DO k=1,klev
     413!        d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
     414!     ENDDO
     415!  ENDIF
    390416!  print*,'iflag_pbl ',d_t_diss
     417!###################################################################
    391418
    392419
    393420! Compuation of stability functions
    394    IF (iflag_pbl/=29) THEN
     421!   IF (iflag_pbl/=29) THEN
    395422      DO k=1,klev
    396423      DO ig=1,ngrid
    … …  
    409436      ENDDO
    410437      ENDDO
    411     ENDIF
     438!    ENDIF
    412439
    413440! Computation of turbulent diffusivities
    414    IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
    415      DO k=2,klev
    416         sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
    417      ENDDO
    418    ELSE
    419      DO k=2,klev
    420         sqrtq(:,k)=sqrt(q2(:,k))
    421      ENDDO
    422    ENDIF
    423    DO k=2,klev
    424    DO ig=1,ngrid
    425       km(ig,k)=leff(ig,k)*sqrtq(ig,k)*sm(ig,k)
    426       kn(ig,k)=km(ig,k)*alpha(ig,k)
    427       kq(ig,k)=leff(ig,k)*zq*0.2
    428 !         print*,q2(ig,k),zq,km(ig,k)
     441!  IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
     442!    DO k=2,klev
     443!       sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
     444!    ENDDO
     445!  ELSE
     446   kq(:,:)=0.
     447   DO k=1,klev
     448      ! Coefficient au milieu des couches pour diffuser la TKE
     449      kq(:,k)=0.5*leff(:,k)*sqrt(q2(:,k))*0.2
    429450   ENDDO
     451
     452#ifdef IOPHYS
     453if (okiophys) then
     454call iophys_ecrit('q2b',klev,'KTE inter','m2/s',q2(:,1:klev))
     455endif
     456#endif
     457
     458  IF (iflag_tke_diff==1) THEN
     459    CALL vdif_q2(timestep, RG, RD, ngrid, plev, pt, kq, q2)
     460  ENDIF
     461
     462   km(:,:)=0.
     463   kn(:,:)=0.
     464   DO k=1,klev
     465      km(:,k)=leff(:,k)*sqrt(q2(:,k))*sm(:,k)
     466      kn(:,k)=km(:,k)*alpha(:,k)
    430467   ENDDO
    431468
    432469
    433 
    434 #ifdef IOPHYS
    435 if (okiophys==1) then
     470#ifdef IOPHYS
     471if (okiophys) then
    436472call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
    437473call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
    … …  
    447483#endif
    448484
     485
    449486ENDIF
    450487
    … …  
    452489!  print*,'OK2'
    453490      RETURN
    454 !====================================================================
    455 !   Yamada 2.0
    456 !====================================================================
    457       if (iflag_pbl.eq.6) then
    458 
    459       do k=2,klev
    460          q2(:,k)=l(:,k)**2*zz(:,k)
    461       enddo
    462 
    463 
    464       else if (iflag_pbl.eq.7) then
    465 !====================================================================
    466 !   Yamada 2.Fournier
    467 !====================================================================
    468 
    469 !  Calcul de l,  km, au pas precedent
    470       do k=2,klev
    471                                                           do ig=1,ngrid
    472 !        print*,'SMML=',sm(ig,k),l(ig,k)
    473          delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
    474          kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
    475          mpre(ig,k)=sqrt(m2(ig,k))
    476 !        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    477                                                           enddo
    478       enddo
    479 
    480       do k=2,klev-1
    481                                                           do ig=1,ngrid
    482         m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
    483         mcstat=sqrt(m2cstat)
    484 
    485 !        print*,'M2 L=',k,mpre(ig,k),mcstat
    486 !
    487 !  -----{puis on ecrit la valeur de q qui annule l'equation de m
    488 !        supposee en q3}
    489 !
    490         IF (k.eq.2) THEN
    491           kmcstat=1.E+0 / mcstat &
    492      &    *( unsdz(ig,k)*kmpre(ig,k+1) &
    493      &                        *mpre(ig,k+1) &
    494      &      +unsdz(ig,k-1) &
    495      &              *cd(ig) &
    496      &              *( sqrt(zu(ig,3)**2+zv(ig,3)**2) &
    497      &                -mcstat/unsdzdec(ig,k) &
    498      &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) &
    499      &      /( unsdz(ig,k)+unsdz(ig,k-1) )
    500         ELSE
    501           kmcstat=1.E+0 / mcstat &
    502      &    *( unsdz(ig,k)*kmpre(ig,k+1) &
    503      &                        *mpre(ig,k+1) &
    504      &      +unsdz(ig,k-1)*kmpre(ig,k-1) &
    505      &                          *mpre(ig,k-1) ) &
    506      &      /( unsdz(ig,k)+unsdz(ig,k-1) )
    507         ENDIF
    508 !       print*,'T2 L=',k,tmp2
    509         tmp2=kmcstat &
    510      &      /( sm(ig,k)/q2(ig,k) ) &
    511      &      /l(ig,k)
    512         q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
    513 !       print*,'Q2 L=',k,q2(ig,k)
    514 !
    515                                                           enddo
    516       enddo
    517 
    518       else if (iflag_pbl==8.or.iflag_pbl==9) then
    519 !====================================================================
    520 !   Yamada 2.5 a la Didi
    521 !====================================================================
    522 
    523 
    524 !  Calcul de l,  km, au pas precedent
    525       do k=2,klev
    526                                                           do ig=1,ngrid
    527 !        print*,'SMML=',sm(ig,k),l(ig,k)
    528          delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
    529          if (delta(ig,k).lt.1.e-20) then
    530 !     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
    531             delta(ig,k)=1.e-20
    532          endif
    533          km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
    534          aa0=(m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
    535          aa1=(m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
    536 ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
    537          aa(ig,k)=aa1*timestep/(delta(ig,k)*l(ig,k))
    538 !     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    539          qpre=sqrt(q2(ig,k))
    540 !        if (iflag_pbl.eq.8 ) then
    541             if (aa(ig,k).gt.0.) then
    542                q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
    543             else
    544                q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
    545             endif
    546 !        else ! iflag_pbl=9
    547 !           if (aa(ig,k)*qpre.gt.0.9) then
    548 !              q2(ig,k)=(qpre*10.)**2
    549 !           else
    550 !              q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
    551 !           endif
    552 !        endif
    553          q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
    554 !     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    555                                                           enddo
    556       enddo
    557 
    558       else if (iflag_pbl>=10) then
    559 
    560 !        print*,'Schema mixte D'
    561 !        print*,'Longueur ',l(:,:)
    562          do k=2,klev-1
    563             l(:,k)=max(l(:,k),1.)
    564             km(:,k)=l(:,k)*sqrt(q2(:,k))*sm(:,k)
    565             q2(:,k)=q2(:,k)+timestep*km(:,k)*m2(:,k)*(1.-rif(:,k))
    566             q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
    567             q2(:,k)=1./(1./sqrt(q2(:,k))+timestep/(2*l(:,k)*b1))
    568             q2(:,k)=q2(:,k)*q2(:,k)
    569          enddo
    570 
    571 
    572       else
    573          CALL abort_physic(modname,'Cas nom prevu dans yamada4',1)
    574 
    575       endif ! Fin du cas 8
    576 
    577 !     print*,'OK8'
    578 
    579 !====================================================================
    580 !   Calcul des coefficients de m�ange
    581 !====================================================================
    582       do k=2,klev
    583 !     print*,'k=',k
    584                                                           do ig=1,ngrid
    585 !abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
    586          zq=sqrt(q2(ig,k))
    587          km(ig,k)=l(ig,k)*zq*sm(ig,k)
    588          kn(ig,k)=km(ig,k)*alpha(ig,k)
    589          kq(ig,k)=l(ig,k)*zq*0.2
    590 !     print*,'KML=',km(ig,k),kn(ig,k)
    591                                                           enddo
    592       enddo
    593 
    594 ! Transport diffusif vertical de la TKE.
    595       if (iflag_pbl.ge.12) then
    596 !       print*,'YAMADA VDIF'
    597         q2(:,1)=q2(:,2)
    598         call vdif_q2(timestep,RG,RD,ngrid,plev,zt,kq,q2)
    599       endif
    600 
    601 !   Traitement des cas noctrunes avec l'introduction d'une longueur
    602 !   minilale.
    603 
    604 !====================================================================
    605 !   Traitement particulier pour les cas tres stables.
    606 !   D'apres Holtslag Boville.
    607 
    608       if (prt_level>1) THEN
    609        print*,'YAMADA4 0'
    610       endif !(prt_level>1) THEN
    611                                                           do ig=1,ngrid
    612       coriol(ig)=1.e-4
    613       pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
    614                                                           enddo
    615 
    616 !     print*,'pblhmin ',pblhmin
    617 !Test a remettre 21 11 02
    618 ! test abd 13 05 02      if(0.eq.1) then
    619       if(1==1) then
    620       if(iflag_pbl==8.or.iflag_pbl==10) then
    621 
    622       do k=2,klev
    623          do ig=1,ngrid
    624             if (teta(ig,2).gt.teta(ig,1)) then
    625                qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
    626                kmin=kap*zlev(ig,k)*qmin
    627             else
    628                kmin=-1. ! kmin n'est utilise que pour les SL stables.
    629             endif
    630             if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
    631 !               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    632 !     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    633                kn(ig,k)=kmin
    634                km(ig,k)=kmin
    635                kq(ig,k)=kmin
    636 !   la longueur de melange est suposee etre l= kap z
    637 !   K=l q Sm d'ou q2=(K/l Sm)**2
    638                q2(ig,k)=(qmin/sm(ig,k))**2
    639             endif
    640          enddo
    641       enddo
    642 
    643       else
    644 
    645       do k=2,klev
    646          do ig=1,ngrid
    647             if (teta(ig,2).gt.teta(ig,1)) then
    648                qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
    649                kmin=kap*zlev(ig,k)*qmin
    650             else
    651                kmin=-1. ! kmin n'est utilise que pour les SL stables.
    652             endif
    653             if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
    654 !               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    655 !     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    656                kn(ig,k)=kmin
    657                km(ig,k)=kmin
    658                kq(ig,k)=kmin
    659 !   la longueur de melange est suposee etre l= kap z
    660 !   K=l q Sm d'ou q2=(K/l Sm)**2
    661                sm(ig,k)=1.
    662                alpha(ig,k)=1.
    663                q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
    664                zq=sqrt(q2(ig,k))
    665                km(ig,k)=l(ig,k)*zq*sm(ig,k)
    666                kn(ig,k)=km(ig,k)*alpha(ig,k)
    667                kq(ig,k)=l(ig,k)*zq*0.2
    668             endif
    669          enddo
    670       enddo
    671       endif
    672 
    673       endif
    674 
    675       if (prt_level>1) THEN
    676        print*,'YAMADA4 1'
    677       endif !(prt_level>1) THEN
    678 !   Diagnostique pour stokage
    679 
    680       if(1.eq.0)then
    681       rino=rif
    682       smyam(1:ngrid,1)=0.
    683       styam(1:ngrid,1)=0.
    684       lyam(1:ngrid,1)=0.
    685       knyam(1:ngrid,1)=0.
    686       w2yam(1:ngrid,1)=0.
    687       t2yam(1:ngrid,1)=0.
    688 
    689       smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
    690       styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
    691       lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
    692       knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
    693 
    694 !   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
    695 
    696       w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 &
    697      &    +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev) &
    698      &    *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
    699 
    700       t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) &
    701      &    *dtetadz(1:ngrid,2:klev)**2 &
    702      &    /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
    703       endif
    704 
    705 !     print*,'OKFIN'
    706       first=.false.
    707       return
    708       end
     491      END
Note: See TracChangeset for help on using the changeset viewer.