/[lmdze]/trunk/libf/phylmd/Thermcell/thermcell.f90
ViewVC logotype

Diff of /trunk/libf/phylmd/Thermcell/thermcell.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC
# Line 8  c    s                  ,pu_therm,pv_the Line 8  c    s                  ,pu_therm,pv_the
8    
9        use dimens_m        use dimens_m
10        use dimphy        use dimphy
11        use YOMCST        use SUPHEC_M
12        IMPLICIT NONE        IMPLICIT NONE
13    
14  c=======================================================================  c=======================================================================
# Line 45  c   ---------- Line 45  c   ----------
45        REAL pu(ngrid,nlay),pduadj(ngrid,nlay)        REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
46        REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)        REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
47        REAL po(ngrid,nlay),pdoadj(ngrid,nlay)        REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
48        REAL pplay(ngrid,nlay)        REAL, intent(in):: pplay(ngrid,nlay)
49        real, intent(in):: pplev(ngrid,nlay+1)        real, intent(in):: pplev(ngrid,nlay+1)
50        real pphi(ngrid,nlay)        real pphi(ngrid,nlay)
51    
# Line 722  c    s         ,'   FM=',fm(ig,l) Line 722  c    s         ,'   FM=',fm(ig,l)
722              endif              endif
723              if(.not.masse(ig,l).ge.1.e-10              if(.not.masse(ig,l).ge.1.e-10
724       s         .or..not.masse(ig,l).le.1.e4) then       s         .or..not.masse(ig,l).le.1.e4) then
 c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l  
 c    s         ,'   M=',masse(ig,l)  
 c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',  
 c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)  
 c     print*,'zlev(ig,l+1),zlev(ig,l)'  
 c    s                ,zlev(ig,l+1),zlev(ig,l)  
 c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'  
 c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)  
725              endif              endif
726              if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then              if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
727  c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l  c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
# Line 933  c             wqd(ig,k)=fm(ig,k)*0.5*(q( Line 925  c             wqd(ig,k)=fm(ig,k)*0.5*(q(
925    
926        return        return
927        end        end
       subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse  
      .    ,fraca,larga  
      .    ,u,v,du,dv,ua,va)  
       use dimens_m  
       use dimphy  
       implicit none  
   
 c=======================================================================  
 c  
 c   Calcul du transport verticale dans la couche limite en presence  
 c   de "thermiques" explicitement representes  
 c   calcul du dq/dt une fois qu'on connait les ascendances  
 c  
 c=======================================================================  
   
   
       integer ngrid,nlay  
928    
       real ptimestep  
       real masse(ngrid,nlay),fm(ngrid,nlay+1)  
       real fraca(ngrid,nlay+1)  
       real larga(ngrid)  
       real entr(ngrid,nlay)  
       real u(ngrid,nlay)  
       real ua(ngrid,nlay)  
       real du(ngrid,nlay)  
       real v(ngrid,nlay)  
       real va(ngrid,nlay)  
       real dv(ngrid,nlay)  
   
       real qa(klon,klev),detr(klon,klev)  
       real wvd(klon,klev+1),wud(klon,klev+1)  
       real gamma0,gamma(klon,klev+1)  
       real dua,dva  
       integer iter  
   
       integer ig,k  
   
 c   calcul du detrainement  
   
       do k=1,nlay  
          do ig=1,ngrid  
             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)  
          enddo  
       enddo  
   
 c   calcul de la valeur dans les ascendances  
       do ig=1,ngrid  
          ua(ig,1)=u(ig,1)  
          va(ig,1)=v(ig,1)  
       enddo  
   
       do k=2,nlay  
          do ig=1,ngrid  
             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  
      s         1.e-5*masse(ig,k)) then  
 c   On itère sur la valeur du coeff de freinage.  
 c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))  
                gamma0=masse(ig,k)  
      s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  
      s         *0.5/larga(ig)  
 c              gamma0=0.  
 c   la première fois on multiplie le coefficient de freinage  
 c   par le module du vent dans la couche en dessous.  
                dua=ua(ig,k-1)-u(ig,k-1)  
                dva=va(ig,k-1)-v(ig,k-1)  
                do iter=1,5  
                   gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)  
                   ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  
      s               +(entr(ig,k)+gamma(ig,k))*u(ig,k))  
      s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))  
                   va(ig,k)=(fm(ig,k)*va(ig,k-1)  
      s               +(entr(ig,k)+gamma(ig,k))*v(ig,k))  
      s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))  
 c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva  
                   dua=ua(ig,k)-u(ig,k)  
                   dva=va(ig,k)-v(ig,k)  
                enddo  
             else  
                ua(ig,k)=u(ig,k)  
                va(ig,k)=v(ig,k)  
                gamma(ig,k)=0.  
             endif  
          enddo  
       enddo  
   
       do k=2,nlay  
          do ig=1,ngrid  
             wud(ig,k)=fm(ig,k)*u(ig,k)  
             wvd(ig,k)=fm(ig,k)*v(ig,k)  
          enddo  
       enddo  
       do ig=1,ngrid  
          wud(ig,1)=0.  
          wud(ig,nlay+1)=0.  
          wvd(ig,1)=0.  
          wvd(ig,nlay+1)=0.  
       enddo  
   
       do k=1,nlay  
          do ig=1,ngrid  
             du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  
      s               -(entr(ig,k)+gamma(ig,k))*u(ig,k)  
      s               -wud(ig,k)+wud(ig,k+1))  
      s               /masse(ig,k)  
             dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)  
      s               -(entr(ig,k)+gamma(ig,k))*v(ig,k)  
      s               -wvd(ig,k)+wvd(ig,k+1))  
      s               /masse(ig,k)  
          enddo  
       enddo  
   
       return  
       end  
929        subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac        subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
930       .    ,q,dq,qa)       .    ,q,dq,qa)
931        use dimens_m        use dimens_m
Legend:
Removed from v.3  
changed lines
  Added in v.38
  ViewVC Help
Powered by ViewVC 1.1.21