Changeset 1849 for LMDZ5/trunk


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

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

Legend:

Unmodified
Added
Removed

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

    r1828 r1849  
    7474       REAL freq_COSP
    7575       LOGICAL :: ok_cosp,ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP
    76        INTEGER :: ip_ebil_phy, iflag_rrtm
     76       INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo
    7777       LOGICAL :: ok_strato
    7878       LOGICAL :: ok_hines
    … …  
    102102     &     , ok_lic_melt, cvl_corr, aer_type                            &
    103103     &     , qsol0, iflag_rrtm, ok_strato,ok_hines,ecrit_LES            &
    104      &     , co2_ppm0
     104     &     , co2_ppm0, iflag_ice_thermo
    105105     
    106106       save /clesphys/

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

    r1836 r1849  
    1       SUBROUTINE concvl (iflag_con,iflag_clos,
     1      SUBROUTINE concvl (iflag_clos,
    22     .             dtime,paprs,pplay,
    33     .             t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
    … …  
    7070c======================================================================
    7171c
     72
     73#include "clesphys.h"
    7274#include "dimensions.h"
    7375c
    74        INTEGER iflag_con,iflag_clos
     76       INTEGER iflag_clos
    7577c
    7678       REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
    … …  
    376378      nloc=klon
    377379      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
    378      $              iflag_con,iflag_mix,iflag_clos,dtime,
     380     $              iflag_con,iflag_mix,iflag_ice_thermo,
     381     $              iflag_clos,dtime,
    379382     :              t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
    380383     $              em_p,em_ph,

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

    r1843 r1849  
    141141  LOGICAL,SAVE :: reevap_ice_omp
    142142  INTEGER,SAVE :: iflag_pdf_omp
     143  INTEGER,SAVE :: iflag_ice_thermo_omp
    143144  REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
    144145  REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
    … …  
    966967  t_glace_max_omp = 273.13
    967968  call getin('t_glace_max',t_glace_max_omp)
     969
     970!
     971!Config Key  = iflag_ice_thermo
     972!Config Desc = 
     973!Config Def  = 0
     974!Config Help =
     975!
     976  iflag_ice_thermo_omp = 0
     977  call getin('iflag_ice_thermo',iflag_ice_thermo_omp)
    968978
    969979!Config Key  = rei_min
    … …  
    16941704    t_glace_min = t_glace_min_omp
    16951705    t_glace_max = t_glace_max_omp
     1706    iflag_ice_thermo = iflag_ice_thermo_omp
    16961707    rei_min = rei_min_omp
    16971708    rei_max = rei_max_omp
    … …  
    19111922  write(lunout,*)' t_glace_min = ',t_glace_min
    19121923  write(lunout,*)' t_glace_max = ',t_glace_max
     1924  write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
    19131925  write(lunout,*)' rei_min = ',rei_min
    19141926  write(lunout,*)' rei_max = ',rei_max

  • 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 

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

    r1650 r1849  
    66     :    ,u1,v1,gz1,th1,th1_wake
    77     :    ,tra1
    8      :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
    9      :    ,h1_wake,lv1_wake,cpn1_wake     ,tv1_wake
     8     :    ,h1     ,lv1, lf1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
     9     :    ,h1_wake,lv1_wake,lf1_wake,cpn1_wake     ,tv1_wake
    1010     :    ,sig1,w01,ptop21
    1111     :    ,Ale1,Alp1
    … …  
    1616     o    ,u,v,gz,th,th_wake
    1717     o    ,tra
    18      o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
    19      o    ,h_wake,lv_wake,cpn_wake    ,tv_wake
     18     o    ,h     ,lv, lf     ,cpn    ,p,ph,tv    ,tp,tvp,clw
     19     o    ,h_wake,lv_wake,lf_wake,cpn_wake    ,tv_wake
    2020     o    ,sig,w0,ptop2
    2121     o    ,Ale,Alp  )
    … …  
    4545      real gz1(len,nd),th1(len,nd),th1_wake(len,nd)
    4646      real tra1(len,nd,ntra)
    47       real h1(len,nd),lv1(len,nd),cpn1(len,nd)
     47      real h1(len,nd),lv1(len,nd),lf1(len,nd),cpn1(len,nd)
    4848      real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
    4949      real tvp1(len,nd),clw1(len,nd)
    5050      real h1_wake(len,nd),lv1_wake(len,nd),cpn1_wake(len,nd)
    51       real tv1_wake(len,nd)
     51      real tv1_wake(len,nd),lf1_wake(len,nd)
    5252      real sig1(len,nd), w01(len,nd), ptop21(len)
    5353      real Ale1(len),Alp1(len)
    … …  
    6565      real gz(len,nd),th(len,nd),th_wake(len,nd)
    6666      real tra(len,nd,ntra)
    67       real h(len,nd),lv(len,nd),cpn(len,nd)
     67      real h(len,nd),lv(len,nd),lf(len,nd),cpn(len,nd)
    6868      real p(len,nd),ph(len,nd+1),tv(len,nd),tp(len,nd)
    6969      real tvp(len,nd),clw(len,nd)
    7070      real h_wake(len,nd),lv_wake(len,nd),cpn_wake(len,nd)
    71       real tv_wake(len,nd)
     71      real tv_wake(len,nd),lf_wake(len,nd)
    7272      real sig(len,nd), w0(len,nd), ptop2(len)
    7373      real Ale(len),Alp(len)
    … …  
    9999        h(nn,k)=h1(i,k)
    100100        lv(nn,k)=lv1(i,k)
     101        lf(nn,k)=lf1(i,k)
    101102        cpn(nn,k)=cpn1(i,k)
    102103        p(nn,k)=p1(i,k)
    … …  
    108109        h_wake(nn,k)=h1_wake(i,k)
    109110        lv_wake(nn,k)=lv1_wake(i,k)
     111        lf_wake(nn,k)=lf1_wake(i,k)
    110112        cpn_wake(nn,k)=cpn1_wake(i,k)
    111113        tv_wake(nn,k)=tv1_wake(i,k)

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

    r1836 r1849  
    345345c   (common cvflag)
    346346
    347        CALL cv_flag
     347       CALL cv_flag(0)
    348348
    349349c -- set thermodynamical constants:
    … …  
    701701
    702702!==================================================================
    703       SUBROUTINE cv_flag
     703      SUBROUTINE cv_flag(iflag_ice_thermo)
    704704      implicit none
     705
     706c Argument : iflag_ice_thermo : ice thermodynamics is taken into account if
     707c                               iflag_ice_thermo >=1
     708      INTEGER iflag_ice_thermo
    705709
    706710#include "cvflag.h"
    … …  
    709713c differente de 10.0 dans convect3:
    710714      cvflag_grav = .TRUE.
     715      cvflag_ice = iflag_ice_thermo .GE. 1
    711716
    712717      return
    … …  
    744749       cpv = RCPV
    745750       cl  = RCW
     751       ci  = RCS
    746752       rrv = RV
    747753       rrd = RD
    748754       lv0 = RLVTT
     755       lf0 = RLSTT-RLVTT
    749756       g   = RG     ! not used in convect3
    750757c ori      t0  = RTT
    … …  
    758765      clmcpv=cl-cpv
    759766      clmcpd=cl-cpd
     767      clmci=cl-ci
    760768      cpdmcp=cpd-cpv
    761769      cpvmcpd=cpv-cpd

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

    r1774 r1849  
    33!
    44      SUBROUTINE cva_driver(len,nd,ndp1,ntra,nloc,
    5      &                   iflag_con,iflag_mix,
     5     &                   iflag_con,iflag_mix,iflag_ice_thermo,
    66     &                   iflag_clos,delt,
    77     &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake,
    … …  
    5252C      iflag_con     Integer        Input        version of convect (3/4)
    5353C      iflag_mix     Integer        Input        version of mixing  (0/1/2)
     54C   iflag_ice_thermo Integer        Input        accounting for ice thermodynamics (0/1)
    5455C      iflag_clos    Integer        Input        version of closure (0/1)
    5556C      delt          Real           Input        time step
    … …  
    155156      integer iflag_con
    156157      integer iflag_mix
     158      integer iflag_ice_thermo
    157159      integer iflag_clos
    158160      real delt
    … …  
    381383      real buoybase1(klon)
    382384
     385      real lf1(klon,klev) ,lf1_wake(klon,klev)
    383386      real lv1(klon,klev) ,lv1_wake(klon,klev)
    384387      real cpn1(klon,klev),cpn1_wake(klon,klev)
    … …  
    419422      real gz(nloc,klev),h(nloc,klev)     ,hm(nloc,klev)
    420423      real               h_wake(nloc,klev),hm_wake(nloc,klev)
    421       real lv(nloc,klev)     ,cpn(nloc,klev)
    422       real lv_wake(nloc,klev),cpn_wake(nloc,klev)
     424      real lv(nloc,klev), lf(nloc,klev), cpn(nloc,klev)
     425      real lv_wake(nloc,klev), lf_wake(nloc,klev), cpn_wake(nloc,klev)
    423426      real p(nloc,klev),ph(nloc,klev+1),tv(nloc,klev)    ,tp(nloc,klev)
    424427      real                              tv_wake(nloc,klev)
    … …  
    430433      real sig(nloc,klev), w0(nloc,klev), ptop2(nloc)
    431434      real hp(nloc,klev), ep(nloc,klev), sigp(nloc,klev)
    432       real frac(nloc), buoy(nloc,klev)
     435      real buoy(nloc,klev)
    433436      real cape(nloc)
    434437      real cin(nloc)
    … …  
    444447      real sigd(nloc)
    445448!      real mp(nloc,klev), qp(nloc,klev), up(nloc,klev), vp(nloc,klev)
    446 !      real wt(nloc,klev), water(nloc,klev), evap(nloc,klev)
     449!      real wt(nloc,klev), water(nloc,klev), evap(nloc,klev), ice(nloc,klev)
    447450!      real b(nloc,klev), sigd(nloc)
    448451!      save mp,qp,up,vp,wt,water,evap,b
    449452      real, save, allocatable :: mp(:,:),qp(:,:),up(:,:),vp(:,:)
    450       real, save, allocatable :: wt(:,:),water(:,:),evap(:,:), b(:,:)
    451 c$OMP THREADPRIVATE(mp,qp,up,vp,wt,water,evap,b)
     453      real, save, allocatable :: wt(:,:),water(:,:),evap(:,:)
     454      real, save, allocatable :: ice(:,:),fondue(:,:), b(:,:)
     455      real, save, allocatable :: frac(:,:), faci(:,:)
     456c$OMP THREADPRIVATE(mp,qp,up,vp,wt,water,evap,ice,fondue,b,frac,faci)
    452457      real  ft(nloc,klev), fq(nloc,klev)
    453458      real ftd(nloc,klev), fqd(nloc,klev)
    … …  
    496501         allocate(mp(nloc,klev), qp(nloc,klev), up(nloc,klev))
    497502         allocate(vp(nloc,klev), wt(nloc,klev), water(nloc,klev))
     503         allocate(ice(nloc,klev), fondue(nloc,klev))
    498504         allocate(evap(nloc,klev), b(nloc,klev))
     505         allocate(frac(nloc,klev), faci(nloc,klev))
    499506         first=.false.
    500507       endif
    … …  
    502509c   (common cvflag)
    503510
    504        CALL cv_flag
     511       CALL cv_flag(iflag_ice_thermo)
    505512
    506513c -- set thermodynamical constants:
    … …  
    610617!       print*,'t1, q1 ',t1,q1
    611618       CALL cv3_prelim(len,nd,ndp1,t1,q1,p1,ph1      ! nd->na
    612      o               ,lv1,cpn1,tv1,gz1,h1,hm1,th1)
     619     o               ,lv1,lf1,cpn1,tv1,gz1,h1,hm1,th1)
    613620   
    614621c
    615622       CALL cv3_prelim(len,nd,ndp1,t1_wake,q1_wake,p1,ph1 ! nd->na
    616      o               ,lv1_wake,cpn1_wake,tv1_wake,gz1_wake
     623     o               ,lv1_wake,lf1_wake,cpn1_wake,tv1_wake,gz1_wake
    617624     o               ,h1_wake,bid,th1_wake)
    618625   
    … …  
    755762     :    ,u1,v1,gz1,th1,th1_wake
    756763     :    ,tra1
    757      :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
    758      :    ,h1_wake,lv1_wake,cpn1_wake     ,tv1_wake
     764     :    ,h1     ,lv1, lf1,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
     765     :    ,h1_wake,lv1_wake,lf1_wake,cpn1_wake     ,tv1_wake
    759766     :    ,sig1,w01,ptop21
    760767     :    ,Ale1,Alp1
    … …  
    765772     o    ,u,v,gz,th,th_wake
    766773     o    ,tra
    767      o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
    768      o    ,h_wake,lv_wake,cpn_wake    ,tv_wake
     774     o    ,h     ,lv, lf     ,cpn    ,p,ph,tv    ,tp,tvp,clw
     775     o    ,h_wake,lv_wake,lf_wake,cpn_wake    ,tv_wake
    769776     o    ,sig,w0,ptop2
    770777     o    ,Ale,Alp  )
    … …  
    800807       CALL cv3_undilute2(nloc,ncum,nd,icb,icbs,nk        !na->nd
    801808     :                        ,tnk,qnk,gznk,hnk,t,q,qs,gz
    802      :                        ,p,h,tv,lv,pbase,buoybase,plcl
    803      o                        ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
     809     :                        ,p,h,tv,lv,lf,pbase,buoybase,plcl
     810     o                        ,inb,tp,tvp,clw,hp,ep,sigp,buoy
     811     :                        ,frac)
    804812   
    805813      endif
    … …  
    816824!-------------------------------------------------------------------
    817825      IF (iflag_con .eq. 3) THEN
     826       if ((iflag_ice_thermo.eq.1).and.(iflag_mix.ne.0)) then
     827         write(*,*) " iflag_ice_thermo==1 requires iflag_mix==0",
     828     &   " but iflag_mix=",iflag_mix, ". Might as well stop here."
     829         stop
     830       endif
    818831       IF (iflag_mix .ge. 1 ) THEN
    819832         CALL zilch(supmax,nloc*klev)   
    … …  
    877890        IF (iflag_mix.eq.0) THEN
    878891         CALL cv3_mixing(nloc,ncum,nd,nd,ntra,icb,nk,inb    ! na->nd
    879      :                       ,ph,t,q,qs,u,v,tra,h,lv,qnk
     892     :                       ,ph,t,q,qs,u,v,tra,h,lv,lf,frac,qnk
    880893     :                       ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
    881894     o   ,ment,qent,uent,vent,nent,sigij,elij,ments,qents,traent)
    … …  
    913926       CALL cv3_unsat(nloc,ncum,nd,nd,ntra,icb,inb,iflag    ! na->nd
    914927     :               ,t_wake,q_wake,qs_wake,gz,u,v,tra,p,ph
    915      :               ,th_wake,tv_wake,lv_wake,cpn_wake
     928     :               ,th_wake,tv_wake,lv_wake,lf_wake,cpn_wake
    916929     :               ,ep,sigp,clw
    917930     :               ,m,ment,elij,delt,plcl,coef_clos
    918      o          ,mp,qp,up,vp,trap,wt,water,evap,b,sigd
     931     o          ,mp,qp,up,vp,trap,wt,water,evap,fondue,ice
     932     :          ,faci,b,sigd
    919933     o          ,wdtrainA,wdtrainM)   ! RomP
    920934      endif
    … …  
    943957     :                     ,icb,inb,delt
    944958     :                     ,t,q,t_wake,q_wake,s_wake,u,v,tra
    945      :                     ,gz,p,ph,h,hp,lv,cpn,th,th_wake
     959     :                     ,gz,p,ph,h,hp,lv,lf,cpn,th,th_wake
    946960     :                     ,ep,clw,m,tp,mp,qp,up,vp,trap
    947      :                     ,wt,water,evap,b,sigd
     961     :                     ,wt,water,ice,evap,fondue,faci,b,sigd
    948962     :                    ,ment,qent,hent,iflag_mix,uent,vent
    949963     :                    ,nent,elij,traent,sig
    … …  
    10431057      return
    10441058      end
    1045 

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

    r766 r1849  
    33!
    44      logical cvflag_grav
     5      logical cvflag_ice
    56
    6       COMMON /cvflag/ cvflag_grav
     7      COMMON /cvflag/ cvflag_grav, cvflag_ice
    78c$OMP THREADPRIVATE(/cvflag/)

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

    r766 r1849  
    44c Thermodynamical constants for convectL:
    55
    6       real cpd, cpv, cl, rrv, rrd, lv0, g, rowl, t0
    7       real clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl
     6      real cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl, t0
     7      real clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl, clmci
    88      real eps, epsi, epsim1
    99      real ginv, hrd
    1010      real grav
    1111
    12       COMMON /cvthermo/ cpd, cpv, cl, rrv, rrd, lv0, g, rowl, t0
    13      :                 ,clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl
    14      :                 ,eps, epsi, epsim1, ginv, hrd, grav
     12      COMMON /cvthermo/ cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl
     13     :                 ,t0, clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl
     14     :                 ,clmci, eps, epsi, epsim1, ginv, hrd, grav
    1515
    1616c$OMP THREADPRIVATE(/cvthermo/)

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

    r1746 r1849  
    88     frac_impa, frac_nucl, beta,                        &
    99     prfl, psfl, rhcl, zqta, fraca,                     &
    10      ztv, zpspsk, ztla, zthl, iflag_cldcon)
     10     ztv, zpspsk, ztla, zthl, iflag_cldcon,             &
     11     iflag_ice_thermo)
    1112
    1213  !
    … …  
    5152  REAL zpspsk(klon,klev),ztla(klon,klev)
    5253  REAL zthl(klon,klev)
     54  REAL ztfondue, qsl, qsi
    5355
    5456  logical lognormale(klon)
     57  logical ice_thermo
    5558
    5659  !AA
    … …  
    7780  INTEGER ncoreczq 
    7881  INTEGER iflag_cldcon
     82  INTEGER iflag_ice_thermo
    7983  PARAMETER (ninter=5)
    8084  LOGICAL evap_prec ! evaporation de la pluie
    … …  
    97101  INTEGER i, k, n, kk
    98102  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
    99   REAL zrfl(klon), zrfln(klon), zqev, zqevt
    100   REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
     103  REAL zrfl(klon), zrfln(klon), zqev, zqevt
     104  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
     105  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
     106  REAL zoliqp(klon), zoliqi(klon)
    101107  REAL ztglace, zt(klon)
    102108  INTEGER nexpo ! exponentiel pour glace/eau
    103109  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
    104110  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
     111  REAL zmelt, zpluie, zice, zcondold
     112  PARAMETER (ztfondue=278.15)
    105113  !
    106114  LOGICAL appel1er
    … …  
    145153  DATA appel1er /.TRUE./
    146154  !ym
     155  ice_thermo = iflag_ice_thermo .GE. 1
    147156  zdelq=0.0
    148157
    … …  
    188197  !
    189198  ztglace = RTT - 15.0
    190   nexpo = 6
     199!AJ<
     200  IF (ice_thermo) THEN
     201    nexpo = 2
     202  ELSE
     203    nexpo = 6
     204  ENDIF
     205!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
     206!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
     207!>AJ
    191208  !cc      nexpo = 1
    192209  !
    … …  
    223240     !     DO i = 1, klon
    224241     zrfl(i) = 0.0
     242     zifl(i) = 0.0
    225243     zneb(i) = seuil_neb
    226244  ENDDO
    … …  
    272290
    273291
     292     ! Calculer l'evaporation de la precipitation
     293     !
     294
     295
    274296     IF (evap_prec) THEN
    275297        DO i = 1, klon
    276            IF (zrfl(i) .GT.0.) THEN
     298!AJ<
     299!!           IF (zrfl(i) .GT.0.) THEN
     300           IF (zrfl(i)+zifl(i).GT.0.) THEN
     301!>AJ
    277302              IF (thermcep) THEN
    278303                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
    … …  
    288313                 ENDIF
    289314              ENDIF
    290               zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    291               zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
    292                    * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    293               zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
    294                    * RG*dtime/(paprs(i,k)-paprs(i,k+1))
    295               zqev = MIN (zqev, zqevt)
    296               zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
    297                    /RG/dtime
    298 
    299               ! pour la glace, on ré-évapore toute la précip dans la
    300               ! couche du dessous
    301               ! la glace venant de la couche du dessus est simplement
    302               ! dans la couche du dessous.
    303 
    304               IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
    305 
    306               zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
    307                    * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    308               zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
    309                    * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    310                    * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    311               zrfl(i) = zrfln(i)
    312            ENDIF
    313         ENDDO
    314      ENDIF
     315           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
     316        ENDDO
     317!AJ<
     318       IF (.NOT. ice_thermo) THEN
     319        DO i = 1, klon
     320!AJ<
     321!!           IF (zrfl(i) .GT.0.) THEN
     322           IF (zrfl(i)+zifl(i).GT.0.) THEN
     323!>AJ
     324                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     325                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
     326                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
     327                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
     328                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
     329                zqev = MIN (zqev, zqevt)
     330                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
     331                     /RG/dtime
     332       
     333                ! pour la glace, on ré-évapore toute la précip dans la
     334                ! couche du dessous
     335                ! la glace venant de la couche du dessus est simplement
     336                ! dans la couche du dessous.
     337       
     338                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
     339       
     340                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
     341                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     342                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     343                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     344                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     345                zrfl(i) = zrfln(i)
     346                zifl(i) = 0.
     347           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
     348        ENDDO
     349!
     350       ELSE ! (.NOT. ice_thermo)
     351!
     352        DO i = 1, klon
     353!AJ<
     354!!           IF (zrfl(i) .GT.0.) THEN
     355           IF (zrfl(i)+zifl(i).GT.0.) THEN
     356!>AJ
     357     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     358     ! Modification de l'évaporation avec la glace
     359     ! Différentiation entre précipitation liquide et solide
     360     ! On suppose que coef_evai=2*coef_eva
     361     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     362     
     363         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     364     !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
     365
     366     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     367     ! On différencie qsat pour l'eau et la glace
     368     ! Si zdelta=1. --> glace
     369     ! Si zdelta=0. --> eau liquide
     370     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     371         
     372         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
     373         qsl= MIN(0.5,qsl)
     374         zcor= 1./(1.-RETV*qsl)
     375         qsl= qsl*zcor
     376         
     377         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
     378              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
     379         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
     380              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
     381
     382         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
     383         qsi= MIN(0.5,qsi)
     384         zcor= 1./(1.-RETV*qsi)
     385         qsi= qsi*zcor
     386
     387         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
     388              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
     389         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
     390              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
     391
     392             
     393     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     394     ! Vérification sur l'évaporation
     395     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     396     
     397         IF (zqevt+zqevti.GT.zqev0) THEN
     398                  zqev=zqev0*zqevt/(zqevt+zqevti)
     399                  zqevi=zqev0*zqevti/(zqevt+zqevti)
     400             
     401         ELSE
     402             IF (zqevt+zqevti.GT.0.) THEN
     403                  zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
     404                  zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
     405             ELSE
     406             zqev=0.
     407             zqevi=0.
     408             ENDIF
     409         ENDIF
     410     
     411         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
     412                                 /RG/dtime)
     413         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
     414                                 /RG/dtime)
     415         
     416     ! Pour la glace, on révapore toute la précip dans la couche du dessous
     417     ! la glace venant de la couche du dessus est simplement dans la couche
     418     ! du dessous.
     419     
     420     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
     421!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
     422         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
     423                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     424         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     425                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     426                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
     427                  + (zifln(i)-zifl(i)) &
     428                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     429                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
     430     
     431         zrfl(i) = zrfln(i)
     432         zifl(i) = zifln(i)
     433
     434           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
     435        ENDDO
     436
     437      ENDIF ! (.NOT. ice_thermo)
     438     
     439     ENDIF ! (evap_prec)
    315440     !
    316441     ! Calculer Qs et L/Cp*dQs/dT:
    … …  
    478603        zq(i) = zq(i) - zcond(i)
    479604        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
    480         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    481      ENDDO
     605     ENDDO
     606!AJ<
     607     IF (.NOT. ice_thermo) THEN
     608       DO i = 1, klon
     609         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     610       ENDDO
     611     ELSE
     612       DO i = 1, klon
     613         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
     614         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     615         zfice(i) = zfice(i)**nexpo
     616         zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
     617                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
     618!         print*,zt(i),zrfl(i),zifl(i),'temp1'
     619       ENDDO
     620     ENDIF
     621!>AJ
    482622     !
    483623     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
    … …  
    488628           zrho(i) = pplay(i,k) / zt(i) / RD
    489629           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
     630        ENDIF
     631     ENDDO
     632!AJ<
     633     IF (.NOT. ice_thermo) THEN
     634     DO i = 1, klon
     635        IF (rneb(i,k).GT.0.0) THEN
    490636           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
    491637           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    492638           zfice(i) = zfice(i)**nexpo
     639     !!      zfice(i)=0.
     640        ENDIF
     641     ENDDO
     642     ENDIF
     643     DO i = 1, klon
     644        IF (rneb(i,k).GT.0.0) THEN
    493645           zneb(i) = MAX(rneb(i,k), seuil_neb)
     646     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
     647!           print*,zt(i),'fractionglace'
     648!>AJ
    494649           radliq(i,k) = zoliq(i)/REAL(ninter+1)
    495650        ENDIF
    … …  
    502657
    503658              IF (zneb(i).EQ.seuil_neb) THEN
     659                 zpluie= 0.0
     660                 zice = 0.0 
    504661                 ztot = 0.0
    505662              ELSE
    … …  
    519676                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
    520677                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
    521                  ztot    = zchau    + zfroi
     678!AJ<
     679                 IF (.NOT. ice_thermo) THEN
     680                   ztot    = zchau + zfroi
     681                 ELSE
     682                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
     683                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
     684                   ztot    = zpluie    + zice
     685                 ENDIF
     686!>AJ
    522687                 ztot    = MAX(ztot   ,0.0)
    523688              ENDIF
    524689              ztot    = MIN(ztot,zoliq(i))
     690!AJ<
     691     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
     692     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
     693              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
     694              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
    525695              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
     696!>AJ
    526697              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
    527698           ENDIF
    … …  
    529700     ENDDO
    530701     !
    531      DO i = 1, klon
    532         IF (rneb(i,k).GT.0.0) THEN
     702         IF (.NOT. ice_thermo) THEN
     703       DO i = 1, klon
     704         IF (rneb(i,k).GT.0.0) THEN
    533705           d_ql(i,k) = zoliq(i)
    534706           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
    535707                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    536         ENDIF
    537         IF (zt(i).LT.RTT) THEN
     708         ENDIF
     709       ENDDO
     710     ELSE
     711       DO i = 1, klon
     712         IF (rneb(i,k).GT.0.0) THEN
     713           d_ql(i,k) = zoliq(i)
     714!AJ<
     715           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
     716               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     717           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
     718                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
     719     !      zrfl(i) = zrfl(i)+  zpluie                         &
     720     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     721     !      zifl(i) = zifl(i)+  zice                    &
     722     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
     723
     724         ENDIF                     
     725       ENDDO
     726     ENDIF
     727
     728     IF (ice_thermo) THEN
     729       DO i = 1, klon
     730           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
     731           zmelt = MIN(MAX(zmelt,0.),1.)
     732           zrfl(i)=zrfl(i)+zmelt*zifl(i)
     733           zifl(i)=zifl(i)*(1.-zmelt)
     734!           print*,zt(i),'octavio1'
     735           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     736                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
     737!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
     738       ENDDO
     739     ENDIF
     740
     741       
     742     IF (.NOT. ice_thermo) THEN
     743       DO i = 1, klon
     744         IF (zt(i).LT.RTT) THEN
    538745           psfl(i,k)=zrfl(i)
    539         ELSE
     746         ELSE
    540747           prfl(i,k)=zrfl(i)
    541         ENDIF
    542      ENDDO
     748         ENDIF
     749       ENDDO
     750     ELSE
     751     ! JAM*************************************************
     752     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
     753     ! *****************************************************
     754
     755       DO i = 1, klon
     756     !   IF (zt(i).LT.RTT) THEN
     757           psfl(i,k)=zifl(i)
     758     !   ELSE
     759           prfl(i,k)=zrfl(i)
     760     !   ENDIF
     761!>AJ
     762       ENDDO
     763     ENDIF
     764     !
    543765     !
    544766     ! Calculer les tendances de q et de t:
    … …  
    604826  DO i = 1, klon
    605827     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
    606         snow(i) = zrfl(i)
     828!AJ<
     829!!        snow(i) = zrfl(i)
     830        snow(i) = zrfl(i)+zifl(i)
     831!>AJ
    607832        zlh_solid(i) = RLSTT-RLVTT
    608833     ELSE

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

    r1828 r1849  
    18781878      DO i = 1, klon
    18791879         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    1880 c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    1881          zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
     1880cjyg<
     1881c      Attention : Arnaud a propose des formules completement differentes
     1882c                  A verifier !!!
     1883         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
     1884         IF (iflag_ice_thermo .EQ. 0) THEN
     1885           zlsdcp=zlvdcp
     1886         ENDIF
     1887c>jyg
     1888
    18821889         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
    18831890         zb = MAX(0.0,ql_seri(i,k))
    … …  
    22632270                nbtr_tmp=nbtr
    22642271             END IF
    2265           CALL concvl (iflag_con,iflag_clos,
     2272cjyg   iflag_con est dans clesphys
     2273cc          CALL concvl (iflag_con,iflag_clos,
     2274          CALL concvl (iflag_clos,
    22662275     .        dtime,paprs,pplay,t_undi,q_undi,
    22672276     .        t_wake,q_wake,wake_s,
    … …  
    27652774     .           frac_impa, frac_nucl, beta_prec_fisrt,
    27662775     .           prfl, psfl, rhcl,
    2767      .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
     2776     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon,
     2777     .           iflag_ice_thermo)
    27682778
    27692779      WHERE (rain_lsc < 0) rain_lsc = 0.
Note: See TracChangeset for help on using the changeset viewer.