Ignore:
Timestamp:
Aug 27, 2013, 12:55:18 PM (11 years ago)
Author:
Ehouarn Millour
Message:

Including the thermodynamics of ice in the convection scheme (iactive only if iflag_ice_thermo==1).
CR+JYG

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cv3_routines.F

    r1774 r1849  
    143143
    144144      SUBROUTINE cv3_prelim(len,nd,ndp1,t,q,p,ph
    145      :                    ,lv,cpn,tv,gz,h,hm,th)
     145     :                    ,lv,lf,cpn,tv,gz,h,hm,th)
    146146      implicit none
    147147
    … …  
    157157
    158158c outputs:
    159       real lv(len,nd), cpn(len,nd), tv(len,nd)
     159      real lv(len,nd), lf(len,nd), cpn(len,nd), tv(len,nd)
    160160      real gz(len,nd), h(len,nd), hm(len,nd)
    161161      real th(len,nd)
    … …  
    178178cdebug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    179179          lv(i,k)= lv0-clmcpv*(t(i,k)-273.15)
     180          lf(i,k)= lf0-clmci*(t(i,k)-273.15)
    180181          cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
    181182          cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
    … …  
    918919      end
    919920
     921      SUBROUTINE Icefrac(t,clw,qi,nl,len)
     922      implicit none
     923
     924
     925cJAM--------------------------------------------------------------------
     926C    Calcul de la quantité d'eau sous forme de glace
     927C--------------------------------------------------------------------
     928      Real qi(len,nl)
     929      Real t(len,nl), clw(len,nl)
     930      Real fracg
     931      Integer nl, len, k, i
     932
     933      do k=3, nl
     934           do i = 1, len
     935           if (t(i,k).gt.263.15) then
     936           qi(i,k)=0.
     937           else
     938           if (t(i,k).lt.243.15) then
     939           qi(i,k)=clw(i,k)
     940           else
     941           fracg=(263.15-t(i,k))/20
     942           qi(i,k)=clw(i,k)*fracg
     943           endif
     944           endif
     945C           print*,t(i,k),qi(i,k),'temp,testglace'
     946           enddo
     947       enddo
     948
     949      return
     950
     951      end
     952
    920953      SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk
    921954     :                       ,tnk,qnk,gznk,hnk,t,q,qs,gz
    922      :                       ,p,h,tv,lv,pbase,buoybase,plcl
    923      o                       ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
     955     :                       ,p,h,tv,lv,lf,pbase,buoybase,plcl
     956     o                       ,inb,tp,tvp,clw,hp,ep,sigp,buoy,frac)
    924957      implicit none
    925958
    … …  
    945978#include "cv3param.h"
    946979#include "conema3.h"
     980#include "cvflag.h"
    947981
    948982c inputs:
    949       integer ncum, nd, nloc
     983      integer ncum, nd, nloc, j
    950984      integer icb(nloc), icbs(nloc), nk(nloc)
    951985      real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
    … …  
    953987      real tnk(nloc), qnk(nloc), gznk(nloc)
    954988      real hnk(nloc)
    955       real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
     989      real lv(nloc,nd), lf(nloc,nd), tv(nloc,nd), h(nloc,nd)
    956990      real pbase(nloc), buoybase(nloc), plcl(nloc)
    957991
    … …  
    964998c local variables:
    965999      integer i, k
    966       real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
    967       real by, defrac, pden
     1000       real tg,qg,ahg,alv,alf,s,tc,es,esi,denom,rg,tca,elacrit
     1001      real als
     1002      real qsat_new,snew, qi(nloc,nd)
     1003      real by, defrac, pden, tbis
    9681004      real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
     1005      real frac(nloc,nd)
    9691006      logical lcape(nloc)
    9701007      integer iposit(nloc)
     1008      Real fracg
    9711009
    9721010!=====================================================================
    … …  
    9781016       ep(i,k)=0.0
    9791017       sigp(i,k)=spfac
     1018       qi(i,k)=0.
    9801019 160  continue
    9811020 170  continue
    … …  
    10601099
    10611100c convect3: no approximation:
    1062            tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
    1063 
     1101           if (cvflag_ice) then
     1102              tp(i,k)=Max(0.,(ah0(i)-gz(i,k)-alv*qg)
     1103     &                /(cpd+(cl-cpd)*qnk(i)))
     1104           else
     1105              tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
     1106           endif
     1107c         
    10641108               clw(i,k)=qnk(i)-qg
    10651109               clw(i,k)=max(0.0,clw(i,k))
    … …  
    10681112c convect3: (qg utilise au lieu du vrai mixing ratio rg):
    10691113               tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing
     1114               if (cvflag_ice) then
     1115               if(clw(i,k).lt.1.e-11) then
     1116               tp(i,k)=tv(i,k)
     1117               tvp(i,k)=tv(i,k)
     1118               endif
     1119               endif
    10701120            endif
     1121
     1122        if (cvflag_ice) then
     1123cCR:attention boucle en klon dans Icefrac         
     1124c          Call Icefrac(t,clw,qi,nl,nloc)
     1125           if (t(i,k).gt.263.15) then
     1126           qi(i,k)=0.
     1127           else
     1128           if (t(i,k).lt.243.15) then
     1129           qi(i,k)=clw(i,k)
     1130           else
     1131           fracg=(263.15-t(i,k))/20
     1132           qi(i,k)=clw(i,k)*fracg
     1133           endif
     1134           endif     
     1135cCR: fin test     
     1136          if(t(i,k).lt.263.15) then
     1137cCR: on commente les calculs d'Arnaud car division par zero
     1138cnouveau calcul propose par JYG
     1139c           alv=lv0-clmcpv*(t(i,k)-273.15)
     1140c           alf=lf0-clmci*(t(i,k)-273.15)
     1141c           tg=tp(i,k)
     1142c           tc=tp(i,k)-273.15
     1143c           denom=243.5+tc
     1144c             do j=1,3
     1145ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1146c         il faudra que esi vienne en argument de la convection
     1147ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1148c             tbis=t(i,k)+(tp(i,k)-tg)
     1149c             esi=exp(23.33086-(6111.72784/tbis)
     1150c     :               +0.15215*log(tbis))
     1151c             qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
     1152c             snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/
     1153c     :                               (rrv*tbis*tbis)
     1154c             snew=1./snew
     1155c             print*,esi,qsat_new,snew,'esi,qsat,snew'
     1156c             tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
     1157c             print*,k,tp(i,k),qnk(i),'avec glace'
     1158c             print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
     1159c             enddo
     1160
     1161           alv=lv0-clmcpv*(t(i,k)-273.15)
     1162           alf=lf0+clmci*(t(i,k)-273.15)
     1163           als=alf+alv
     1164           tg=tp(i,k)
     1165           tp(i,k) = t(i,k)
     1166           do j=1,3
     1167             esi=exp(23.33086-(6111.72784/tp(i,k))
     1168     :               +0.15215*log(tp(i,k)))
     1169             qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
     1170             snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*als*qsat_new/
     1171     :                               (rrv*tp(i,k)*tp(i,k))
     1172             snew=1./snew
     1173cc             print*,esi,qsat_new,snew,'esi,qsat,snew'
     1174             tp(i,k)=tp(i,k)+
     1175     :         ( (cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) +
     1176     :         alv*(qg-qsat_new) + alf*qi(i,k) )*snew
     1177c             print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k),
     1178c     :             'k,tp,q,qt,qi avec glace'
     1179           enddo
     1180
     1181cCR:reprise du code AJ
     1182           clw(i,k)=qnk(i)-qsat_new
     1183           clw(i,k)=max(0.0,clw(i,k))
     1184           tvp(i,k)=Max(0.,tp(i,k)*(1.+qsat_new/eps-qnk(i)))
     1185c           print*,tvp(i,k),'tvp'
     1186          endif
     1187          if(clw(i,k).lt.1.e-11) then
     1188           tp(i,k)=tv(i,k)
     1189           tvp(i,k)=tv(i,k)
     1190          endif
     1191        endif  ! (cvflag_ice)
     1192
    10711193  290     continue
    10721194  300   continue
    … …  
    13031425        do 590 i=1,ncum
    13041426        if((k.ge.icb(i)).and.(k.le.inb(i)))then
     1427
     1428          if (cvflag_ice) then
     1429          frac(i,k)=1.-(t(i,k)-243.15)/(263.15-243.15)
     1430          frac(i,k)=min(max(frac(i,k),0.0),1.0)
     1431          hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))
     1432     :            *ep(i,k)*clw(i,k)
     1433
     1434          else
    13051435          hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
     1436          endif
     1437
    13061438        endif
    13071439 590    continue
    … …  
    15551687
    15561688      SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
    1557      :                    ,ph,t,rr,rs,u,v,tra,h,lv,qnk
     1689     :                    ,ph,t,rr,rs,u,v,tra,h,lv,lf,frac,qnk
    15581690     :                    ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
    15591691     :   ,ment,qent,uent,vent,nent,sij,elij,ments,qents,traent)
    … …  
    15671699#include "cvthermo.h"
    15681700#include "cv3param.h"
     1701#include "cvflag.h"
    15691702
    15701703c inputs:
    … …  
    15781711      real tra(nloc,nd,ntra) ! input of convect3
    15791712      real lv(nloc,na), h(nloc,na), hp(nloc,na)
     1713      real lf(nloc,na), frac(nloc,na)
    15801714      real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
    15811715      real m(nloc,na)        ! input of convect3
    … …  
    16591793          rti=qnk(il)-ep(il,i)*clw(il,i)
    16601794          bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
     1795
     1796
     1797          if (cvflag_ice) then
     1798c          print*,cvflag_ice,'cvflag_ice dans do 700'
     1799          if (t(il,j).le.263.15) then
     1800          bf2=1.+(lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)
     1801     :          *lf(il,j))*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
     1802          endif
     1803          endif
     1804
    16611805          anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
    16621806          denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
    … …  
    16711815          if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
    16721816     :                 .and.j.gt.i)then
     1817
     1818          if (cvflag_ice) then
     1819          anum=anum-(lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)
     1820     :          -cwat*bf2)
     1821          denom=denom+(lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
     1822          else
    16731823           anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
    16741824           denom=denom+lv(il,j)*(rr(il,i)-rti)
     1825          endif
     1826
    16751827           if(abs(denom).lt.0.01)denom=0.01
    16761828           sij(il,i,j)=anum/denom
    … …  
    17801932        lwork(il)=(nent(il,i).ne.0)
    17811933        qp=qnk(il)-ep(il,i)*clw(il,i)
     1934
     1935        if(cvflag_ice) then
     1936
     1937        anum=h(il,i)-hp(il,i)-(lv(il,i)+frac(il,i)*lf(il,i))*
     1938     :       (qp-rs(il,i))+(cpv-cpd)*t(il,i)*(qp-rr(il,i))
     1939        denom=h(il,i)-hp(il,i)+(lv(il,i)+frac(il,i)*lf(il,i))*
     1940     :        (rr(il,i)-qp)+(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
     1941        else
     1942
    17821943        anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
    17831944     :           +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
    17841945        denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
    17851946     :           +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
     1947        endif
     1948
    17861949        if(abs(denom).lt.0.01)denom=0.01
    17871950        scrit(il)=anum/denom
    … …  
    19482111      SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb,iflag
    19492112     :              ,t,rr,rs,gz,u,v,tra,p,ph
    1950      :              ,th,tv,lv,cpn,ep,sigp,clw
     2113     :              ,th,tv,lv,lf,cpn,ep,sigp,clw
    19512114     :              ,m,ment,elij,delt,plcl,coef_clos
    1952      o              ,mp,rp,up,vp,trap,wt,water,evap,b,sigd
     2115     o              ,mp,rp,up,vp,trap,wt,water,evap,fondue,ice
     2116     :              ,faci,b,sigd
    19532117     o              ,wdtrainA,wdtrainM)                                ! RomP
    19542118      implicit none
    … …  
    19692133      real ep(nloc,na), sigp(nloc,na), clw(nloc,na)
    19702134      real th(nloc,na),tv(nloc,na),lv(nloc,na),cpn(nloc,na)
     2135      real lf(nloc,na)
    19712136      real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
    19722137      real coef_clos(nloc)
    … …  
    19782143      real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
    19792144      real water(nloc,na), evap(nloc,na), wt(nloc,na)
     2145      real ice(nloc,na), fondue(nloc,na),faci(nloc,na)
    19802146      real trap(nloc,na,ntra)
    19812147      real b(nloc,na), sigd(nloc)
    … …  
    19882154c local variables
    19892155      integer i,j,k,il,num1,ndp1
    1990       real tinv, delti
     2156      real tinv, delti, coef
    19912157      real awat, afac, afac1, afac2, bfac
    1992       real pr1, pr2, sigt, b6, c6, revap, delth
     2158      real pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
    19932159      real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
    1994       real ampmax
     2160      real ampmax, thaw
    19952161      real tevap(nloc)
    1996       real lvcp(nloc,na)
     2162      real lvcp(nloc,na),lfcp(nloc,na)
    19972163      real h(nloc,na),hm(nloc,na)
     2164      real frac(nloc,na)
     2165      real fraci(nloc,na),prec(nloc,na)
    19982166      real wdtrain(nloc)
    19992167      logical lwork(nloc),mplus(nloc)
    … …  
    20152183          wt(il,i)=0.001
    20162184          water(il,i)=0.0
     2185          frac(il,i)=0.0
     2186          faci(il,i)=0.0
     2187          fraci(il,i)=0.0
     2188          ice(il,i)=0.0
     2189          prec(il,i)=0.0
     2190          fondue(il,i)=0.0
    20172191          evap(il,i)=0.0
    20182192          b(il,i)=0.0
    20192193          lvcp(il,i)=lv(il,i)/cpn(il,i)
     2194          lfcp(il,i)=lf(il,i)/cpn(il,i)
    20202195         enddo
    20212196        enddo
    … …  
    21152290      wt(il,i)=45.0
    21162291
     2292      if (cvflag_ice) then
     2293      frac(il,inb(il)) = 1. -(t(il,inb(il))-243.15)/(263.15-243.15)
     2294      frac(il,inb(il)) = min(max(frac(il,inb(il)),0.),1.)
     2295      fraci(il,inb(il)) = frac(il,inb(il))
     2296      else
     2297      continue
     2298      endif
     2299
    21172300      if(i.lt.inb(il))then
     2301
     2302      if (cvflag_ice) then
     2303      thaw = (t(il,i)-273.15)/(275.15-273.15)
     2304      thaw = min(max(thaw,0.0),1.0)
     2305      frac(il,i)=frac(il,i)*(1.-thaw)
     2306      else
     2307      continue
     2308      endif
     2309
    21182310       rp(il,i)=rp(il,i+1)
    21192311     :       +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
    21202312       rp(il,i)=0.5*(rp(il,i)+rr(il,i))
    21212313      endif
     2314      fraci(il,i)=1.-(t(il,i)-243.15)/(263.15-243.15)
     2315      fraci(il,i)=min(max(fraci(il,i),0.0),1.0)
    21222316      rp(il,i)=max(rp(il,i),0.0)
    21232317      rp(il,i)=amin1(rp(il,i),rs(il,i))
    … …  
    21262320      if(i.eq.1)then
    21272321       afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))
     2322       if (cvflag_ice) then
     2323       afac1=p(il,i)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))
     2324       endif
    21282325      else
    21292326       rp(il,i-1)=rp(il,i)
    … …  
    21612358c       b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    21622359c       c6 = water(il,i+1) + wdtrain(il)*bfac
     2360c       c6 = prec(il,i+1) + wdtrain(il)*bfac
    21632361c        revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    21642362c        evap(il,i)=sigt*afac*revap
    21652363c        water(il,i)=revap*revap
     2364c        prec(il,i)=revap*revap
    21662365cc        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ',
    21672366cc     $            i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
    … …  
    21692368c
    21702369c--------retour à la formulation originale d'Emanuel.
     2370       if (cvflag_ice) then
     2371
     2372c      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2373c      c6=prec(il,i+1)+bfac*wdtrain(il)
     2374c     :    -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
     2375c      if(c6.gt.0.0)then
     2376c       revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2377c
     2378cJAM  Attention: evap=sigt*E
     2379c     Modification: evap devient l'évaporation en milieu de couche
     2380c     car nécessaire dans cv3_yield
     2381c     Du coup, il faut modifier pas mal d'équations...
     2382c     et l'expression de afac qui devient afac1
     2383c     revap=sqrt((prec(i+1)+prec(i))/2)
     2384
     2385      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
     2386      c6=prec(il,i+1)+0.5*bfac*wdtrain(il)
     2387c      print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
     2388c      print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
     2389c      print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
     2390      if (c6 .gt. b6*b6+1.e-20) then
     2391       revap=2.*c6/(b6+sqrt(b6*b6+4.*c6))
     2392      else
     2393       revap=(-b6+sqrt(b6*b6+4.*c6))/2.
     2394      endif
     2395      prec(il,i)=Max(0.,2.*revap*revap-prec(il,i+1))
     2396c      print*,prec(il,i),'neige'
     2397
     2398cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
     2399cc             evap(il,i)=sigt*afac*revap
     2400c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
     2401c     Ici,l'evaporation evap est simplement calculee par l'equation de
     2402c     conservation.
     2403c       prec(il,i)=revap*revap
     2404c      else
     2405cjyg----   Correction : si c6 <= 0, water(il,i)=0.
     2406c       prec(il,i)=0.
     2407c      endif
     2408
     2409cjyg---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
     2410c                                    moins [tt ce qui sort de la couche i]
     2411c       print *, 'evap avec ice'
     2412       evap(il,i)=
     2413     :       (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i)))
     2414     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
     2415
     2416      d6=bfac*wdtrain(il)-100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))
     2417     :   *evap(il,i)
     2418      e6=bfac*wdtrain(il)
     2419      f6=-100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*
     2420     :   evap(il,i)
     2421
     2422      thaw=(t(il,i)-273.15)/(275.15-273.15)
     2423      thaw=min(max(thaw,0.0),1.0)
     2424      water(il,i)=water(il,i+1)+(1-fraci(il,i))*d6
     2425      water(il,i)=max(water(il,i),0.)
     2426      ice(il,i)=ice(il,i+1)+fraci(il,i)*d6
     2427      ice(il,i)=max(ice(il,i),0.)
     2428      fondue(il,i)=ice(il,i)*thaw
     2429      water(il,i)=water(il,i)+fondue(il,i)
     2430      ice(il,i)=ice(il,i)-fondue(il,i)
     2431
     2432      if(water(il,i)+ice(il,i).lt.1.e-30)then
     2433      faci(il,i)=0.
     2434      else
     2435      faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
     2436      endif
     2437
     2438c      water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
     2439c      water(il,i)=max(water(il,i),0.)
     2440c      ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
     2441c      ice(il,i)=max(ice(il,i),0.)
     2442c      fondue(il,i)=ice(il,i)*thaw
     2443c      water(il,i)=water(il,i)+fondue(il,i)
     2444c      ice(il,i)=ice(il,i)-fondue(il,i)
     2445c
     2446c      if((water(il,i)+ice(il,i)).lt.1.e-30)then
     2447c       faci(il,i)=0.
     2448c      else
     2449c       faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
     2450c      endif
     2451
     2452      else
    21712453      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    21722454      c6=water(il,i+1)+bfac*wdtrain(il)
    … …  
    21742456      if(c6.gt.0.0)then
    21752457       revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    2176 cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
    2177 cc             evap(il,i)=sigt*afac*revap
    2178 c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
    2179 c     Ici,l'evaporation evap est simplement calculee par l'equation de
    2180 c     conservation.
    21812458       water(il,i)=revap*revap
    21822459      else
    2183 cjyg----   Correction : si c6 <= 0, water(il,i)=0.
    2184        water(il,i) = 0.
     2460       water(il,i)= 0.
    21852461      endif
    2186 cJYG/IM : ci-dessous formulation originale de KE
    2187 c      evap(il,i)=-evap(il,i+1)
    2188 c    :            +(wdtrain(il)+sigd(il)*wt(il,i)*water(il,i+1))
    2189 c    :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*50.)
    2190 c
    2191 cJYG/IM : ci-dessous modification formulation originale de KE
    2192 c        pour eliminer oscillations verticales de pluie se produisant
    2193 c        lorsqu'il y a evaporation totale de la pluie
    2194 c
    2195 c       evap(il,i)= +(wdtrain(il)+sigd(il)*wt(il,i)*water(il,i+1)) !itlmd(jyg)
    2196 c     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    2197 c      end if  !itlmd(jyg)
    2198 cjyg---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
    2199 c                                    moins [tt ce qui sort de la couche i]
     2462c       print *, 'evap sans ice'
    22002463       evap(il,i)=
    22012464     :       (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))
    22022465     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    2203 c
     2466
     2467       endif
    22042468       endif !(i.le.inb(il) .and. lwork(il))
    22052469995   Continue
    … …  
    22152479      tevap(il)=max(0.0,evap(il,i))
    22162480      delth=max(0.001,(th(il,i)-th(il,i-1)))
     2481      if (cvflag_ice) then
     2482      if (cvflag_grav) then
     2483       mp(il,i)=100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)
     2484     :          *(p(il,i-1)-p(il,i))/delth
     2485     :          +lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)
     2486     :          *(p(il,i-1)-p(il,i))/delth
     2487     :          +lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)
     2488     :          *(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
     2489      else
     2490       mp(il,i)=10.*(lvcp(il,i)*sigd(il)*tevap(il)
     2491     :          *(p(il,i-1)-p(il,i))/delth
     2492     :          +lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)
     2493     :          *(p(il,i-1)-p(il,i))/delth
     2494     :          +lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)
     2495     :          *(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
     2496
     2497      endif
     2498      else
    22172499      if (cvflag_grav) then
    22182500       mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il)
    … …  
    22232505      endif
    22242506c
     2507      endif
     2508c
    22252509       endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    22262510996   Continue
    … …  
    22432527     :               /(lvcp(il,i)*sigd(il)*th(il,i))
    22442528       af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
     2529
     2530       if (cvflag_ice) then
     2531       bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
     2532     :            +50.*(p(il,i-1)-p(il,i))*xf*
     2533     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
     2534     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/
     2535     :            (ph(il,i)-ph(il,i+1)))
     2536       else
     2537
    22452538       bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
    22462539     :            +50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
     2540       endif
     2541
    22472542       fac2=1.0
    22482543       if(bf.lt.0.0)fac2=-1.0
    … …  
    22622557       mp(il,i)=max(0.0,mp(il,i))
    22632558
     2559       if (cvflag_ice) then
    22642560       if (cvflag_grav) then
    22652561Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
    22662562C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
    22672563C Et il faut bien revoir les facteurs 100.
    2268         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
     2564         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*
     2565     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
     2566     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)
     2567     :            /(ph(il,i)-ph(il,i+1)))
     2568     2            /(mp(il,i)+sigd(il)*0.1)
     2569     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
     2570     : *sigd(il)*th(il,i))
     2571       else
     2572        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*
     2573     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
     2574     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)
     2575     :            /(ph(il,i)-ph(il,i+1)))
     2576     2            /(mp(il,i)+sigd(il)*0.1)
     2577     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
     2578     : *sigd(il)*th(il,i))
     2579       endif
     2580       else
     2581       if (cvflag_grav) then
     2582         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
    22692583     2   /(mp(il,i)+sigd(il)*0.1)
    22702584     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
    … …  
    22762590     : *sigd(il)*th(il,i))
    22772591       endif
     2592      endif
    22782593       b(il,i-1)=max(b(il,i-1),0.0)
    22792594c
    … …  
    23912706     :                    ,icb,inb,delt
    23922707     :                    ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
    2393      :                    ,gz,p,ph,h,hp,lv,cpn,th,th_wake
     2708     :                    ,gz,p,ph,h,hp,lv,lf,cpn,th,th_wake
    23942709     :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
    2395      :                    ,wt,water,evap,b,sigd
     2710     :                    ,wt,water,ice,evap,fondue,faci,b,sigd
    23962711     :                    ,ment,qent,hent,iflag_mix,uent,vent
    23972712     :                    ,nent,elij,traent,sig
    … …  
    24222737      real th(nloc,na), p(nloc,nd), tp(nloc,na)
    24232738      real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
     2739      real lf(nloc,na)
    24242740      real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
    24252741      real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
    24262742      real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc)
     2743      real fondue(nloc,na),faci(nloc,na), ice(nloc,na)
    24272744      real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
    24282745      real hent(nloc,na,na)
    … …  
    24552772      real cpinv, rdcp, dpinv
    24562773      real awat(nloc)
    2457       real lvcp(nloc,na), mke(nloc,na)
     2774      real lvcp(nloc,na), lfcp(nloc,na), mke(nloc,na)
    24582775      real am(nloc), work(nloc), ad(nloc), amp1(nloc)
    24592776c!!      real up1(nloc), dn1(nloc)
    … …  
    25132830       do il=1,ncum
    25142831         lvcp(il,i)=lv(il,i)/cpn(il,i)
     2832         lfcp(il,i)=lf(il,i)/cpn(il,i)
    25152833       enddo
    25162834      enddo
    … …  
    25222840      do il=1,ncum
    25232841       if(ep(il,inb(il)).ge.0.0001 .and. iflag(il) .le. 1)then
     2842        if (cvflag_ice) then
     2843        if (cvflag_grav) then
     2844           precip(il)=wt(il,1)*sigd(il)*(water(il,1)+ice(il,1))*86400.
     2845     :               *1000./(rowl*grav)
     2846        else
     2847         precip(il)=wt(il,1)*sigd(il)*(water(il,1)+ice(il,1))*8640.
     2848        endif
     2849        else
    25242850        if (cvflag_grav) then
    25252851           precip(il)=wt(il,1)*sigd(il)*water(il,1)*86400.*1000.
    … …  
    25272853        else
    25282854         precip(il)=wt(il,1)*sigd(il)*water(il,1)*8640.
     2855        endif
    25292856        endif
    25302857       endif
    … …  
    25392866       if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il)
    25402867     :    .and. iflag(il) .le. 1)then
     2868        if (cvflag_ice) then
     2869        if (cvflag_grav) then
     2870           VPrecip(il,i) = wt(il,i)*sigd(il)*(water(il,i)+ice(il,i))
     2871     :                     /grav
     2872        else
     2873           VPrecip(il,i) = wt(il,i)*sigd(il)*(water(il,i)+ice(il,i))
     2874     :                     /10.
     2875        endif
     2876        else
    25412877        if (cvflag_grav) then
    25422878           VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/grav
    25432879        else
    25442880           VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/10.
     2881        endif
    25452882        endif
    25462883       endif
    … …  
    25852922cjyg  Correction pour conserver l'eau
    25862923ccc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2)) !precip
    2587        ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)                  !precip
     2924       if (cvflag_ice) then
     2925       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)
     2926     :           -lfcp(il,1)*sigd(il)*evap(il,1)*faci(il,1)
     2927     :           -lfcp(il,1)*sigd(il)*(fondue(il,1)*wt(il,1))/
     2928     :           (100.*(ph(il,1)-ph(il,2)))                               !precip
     2929       else
     2930       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)
     2931       endif
    25882932
    25892933      if (cvflag_grav) then
    … …  
    25952939      endif
    25962940
     2941      if (cvflag_ice) then
     2942        ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2)
     2943     :     *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)+
     2944     :     0.01*sigd(il)*wt(il,1)*(ci-cpd)*ice(il,2)*
     2945     :     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)
     2946      else
    25972947      ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2)
    25982948     :     *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)
     2949      endif
    25992950
    26002951      ftd(il,1) = ft(il,1)                        ! fin precip
    … …  
    27913142       ! precip
    27923143ccc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
     3144       if (cvflag_ice) then
    27933145       ft(il,i)= -sigd(il)*lvcp(il,i)*evap(il,i)
     3146     :           -sigd(il)*lfcp(il,i)*evap(il,i)*faci(il,i)
     3147     :           -sigd(il)*lfcp(il,i)*fondue(il,i)*wt(il,i)/
     3148     :           (100.*(p(il,i-1)-p(il,i)))
     3149       else
     3150       ft(il,i)= -sigd(il)*lvcp(il,i)*evap(il,i)
     3151       endif
     3152
    27943153        rat=cpn(il,i-1)*cpinv
    27953154c
    … …  
    27983157     :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
    27993158     :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
     3159       if (cvflag_ice) then
    28003160       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
    28013161     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3162     :           +0.01*sigd(il)*wt(il,i)*(ci-cpd)*ice(il,i+1)*
     3163     :           (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3164      else
     3165       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
     3166     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3167      endif
     3168
    28023169       ftd(il,i)=ft(il,i)
    28033170        ! fin precip
    … …  
    28183185     :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
    28193186     :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
     3187
     3188       if (cvflag_ice) then
    28203189       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
    28213190     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3191     :           +0.01*sigd(il)*wt(il,i)*(ci-cpd)*ice(il,i+1)*
     3192     :           (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3193       else
     3194       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
     3195     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3196       endif
     3197
    28223198       ftd(il,i)=ft(il,i)
    28233199        ! fin precip
    … …  
    37444120        return
    37454121        end
    3746 
Note: See TracChangeset for help on using the changeset viewer.