Changeset 1864 for LMDZ5/branches/testing/libf/phylmd/cv3_routines.F
- Timestamp:
- Sep 11, 2013, 11:45:01 AM (11 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 1797-1799,1801-1811,1813-1834,1836,1838-1840,1842-1860
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/cv3_routines.F
r1795 r1864 143 143 144 144 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) 146 146 implicit none 147 147 … … 157 157 158 158 c 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) 160 160 real gz(len,nd), h(len,nd), hm(len,nd) 161 161 real th(len,nd) … … 178 178 cdebug lv(i,k)= lv0-clmcpv*(t(i,k)-t0) 179 179 lv(i,k)= lv0-clmcpv*(t(i,k)-273.15) 180 lf(i,k)= lf0-clmci*(t(i,k)-273.15) 180 181 cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k) 181 182 cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k) … … 918 919 end 919 920 921 SUBROUTINE Icefrac(t,clw,qi,nl,len) 922 implicit none 923 924 925 cJAM-------------------------------------------------------------------- 926 C Calcul de la quantité d'eau sous forme de glace 927 C-------------------------------------------------------------------- 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 945 C print*,t(i,k),qi(i,k),'temp,testglace' 946 enddo 947 enddo 948 949 return 950 951 end 952 920 953 SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk 921 954 : ,tnk,qnk,gznk,hnk,t,q,qs,gz 922 : ,p,h,tv,lv, pbase,buoybase,plcl923 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) 924 957 implicit none 925 958 … … 945 978 #include "cv3param.h" 946 979 #include "conema3.h" 980 #include "cvflag.h" 947 981 948 982 c inputs: 949 integer ncum, nd, nloc 983 integer ncum, nd, nloc, j 950 984 integer icb(nloc), icbs(nloc), nk(nloc) 951 985 real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd) … … 953 987 real tnk(nloc), qnk(nloc), gznk(nloc) 954 988 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) 956 990 real pbase(nloc), buoybase(nloc), plcl(nloc) 957 991 … … 964 998 c local variables: 965 999 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 968 1004 real ah0(nloc), cape(nloc), capem(nloc), byp(nloc) 1005 real frac(nloc,nd) 969 1006 logical lcape(nloc) 970 1007 integer iposit(nloc) 1008 Real fracg 971 1009 972 1010 !===================================================================== … … 978 1016 ep(i,k)=0.0 979 1017 sigp(i,k)=spfac 1018 qi(i,k)=0. 980 1019 160 continue 981 1020 170 continue … … 1060 1099 1061 1100 c 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 1107 c 1064 1108 clw(i,k)=qnk(i)-qg 1065 1109 clw(i,k)=max(0.0,clw(i,k)) … … 1068 1112 c convect3: (qg utilise au lieu du vrai mixing ratio rg): 1069 1113 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 1070 1120 endif 1121 1122 if (cvflag_ice) then 1123 cCR:attention boucle en klon dans Icefrac 1124 c 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 1135 cCR: fin test 1136 if(t(i,k).lt.263.15) then 1137 cCR: on commente les calculs d'Arnaud car division par zero 1138 cnouveau calcul propose par JYG 1139 c alv=lv0-clmcpv*(t(i,k)-273.15) 1140 c alf=lf0-clmci*(t(i,k)-273.15) 1141 c tg=tp(i,k) 1142 c tc=tp(i,k)-273.15 1143 c denom=243.5+tc 1144 c do j=1,3 1145 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1146 c il faudra que esi vienne en argument de la convection 1147 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1148 c tbis=t(i,k)+(tp(i,k)-tg) 1149 c esi=exp(23.33086-(6111.72784/tbis) 1150 c : +0.15215*log(tbis)) 1151 c qsat_new=eps*esi/(p(i,k)-esi*(1.-eps)) 1152 c snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/ 1153 c : (rrv*tbis*tbis) 1154 c snew=1./snew 1155 c print*,esi,qsat_new,snew,'esi,qsat,snew' 1156 c tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew 1157 c print*,k,tp(i,k),qnk(i),'avec glace' 1158 c print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew 1159 c 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 1173 cc 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 1177 c print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), 1178 c : 'k,tp,q,qt,qi avec glace' 1179 enddo 1180 1181 cCR: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))) 1185 c 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 1071 1193 290 continue 1072 1194 300 continue … … 1303 1425 do 590 i=1,ncum 1304 1426 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 1305 1435 hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k) 1436 endif 1437 1306 1438 endif 1307 1439 590 continue … … 1555 1687 1556 1688 SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb 1557 : ,ph,t,rr,rs,u,v,tra,h,lv, qnk1689 : ,ph,t,rr,rs,u,v,tra,h,lv,lf,frac,qnk 1558 1690 : ,unk,vnk,hp,tv,tvp,ep,clw,m,sig 1559 1691 : ,ment,qent,uent,vent,nent,sij,elij,ments,qents,traent) … … 1567 1699 #include "cvthermo.h" 1568 1700 #include "cv3param.h" 1701 #include "cvflag.h" 1569 1702 1570 1703 c inputs: … … 1578 1711 real tra(nloc,nd,ntra) ! input of convect3 1579 1712 real lv(nloc,na), h(nloc,na), hp(nloc,na) 1713 real lf(nloc,na), frac(nloc,na) 1580 1714 real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na) 1581 1715 real m(nloc,na) ! input of convect3 … … 1659 1793 rti=qnk(il)-ep(il,i)*clw(il,i) 1660 1794 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 1798 c 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 1661 1805 anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j)) 1662 1806 denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j) … … 1671 1815 if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) 1672 1816 : .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 1673 1823 anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2) 1674 1824 denom=denom+lv(il,j)*(rr(il,i)-rti) 1825 endif 1826 1675 1827 if(abs(denom).lt.0.01)denom=0.01 1676 1828 sij(il,i,j)=anum/denom … … 1780 1932 lwork(il)=(nent(il,i).ne.0) 1781 1933 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 1782 1943 anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) 1783 1944 : +(cpv-cpd)*t(il,i)*(qp-rr(il,i)) 1784 1945 denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) 1785 1946 : +(cpd-cpv)*t(il,i)*(rr(il,i)-qp) 1947 endif 1948 1786 1949 if(abs(denom).lt.0.01)denom=0.01 1787 1950 scrit(il)=anum/denom … … 1948 2111 SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb,iflag 1949 2112 : ,t,rr,rs,gz,u,v,tra,p,ph 1950 : ,th,tv,lv, cpn,ep,sigp,clw2113 : ,th,tv,lv,lf,cpn,ep,sigp,clw 1951 2114 : ,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 1953 2117 o ,wdtrainA,wdtrainM) ! RomP 1954 2118 implicit none … … 1969 2133 real ep(nloc,na), sigp(nloc,na), clw(nloc,na) 1970 2134 real th(nloc,na),tv(nloc,na),lv(nloc,na),cpn(nloc,na) 2135 real lf(nloc,na) 1971 2136 real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na) 1972 2137 real coef_clos(nloc) … … 1978 2143 real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na) 1979 2144 real water(nloc,na), evap(nloc,na), wt(nloc,na) 2145 real ice(nloc,na), fondue(nloc,na),faci(nloc,na) 1980 2146 real trap(nloc,na,ntra) 1981 2147 real b(nloc,na), sigd(nloc) … … 1988 2154 c local variables 1989 2155 integer i,j,k,il,num1,ndp1 1990 real tinv, delti 2156 real tinv, delti, coef 1991 2157 real awat, afac, afac1, afac2, bfac 1992 real pr1, pr2, sigt, b6, c6, revap, delth2158 real pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth 1993 2159 real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf 1994 real ampmax 2160 real ampmax, thaw 1995 2161 real tevap(nloc) 1996 real lvcp(nloc,na) 2162 real lvcp(nloc,na),lfcp(nloc,na) 1997 2163 real h(nloc,na),hm(nloc,na) 2164 real frac(nloc,na) 2165 real fraci(nloc,na),prec(nloc,na) 1998 2166 real wdtrain(nloc) 1999 2167 logical lwork(nloc),mplus(nloc) … … 2015 2183 wt(il,i)=0.001 2016 2184 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 2017 2191 evap(il,i)=0.0 2018 2192 b(il,i)=0.0 2019 2193 lvcp(il,i)=lv(il,i)/cpn(il,i) 2194 lfcp(il,i)=lf(il,i)/cpn(il,i) 2020 2195 enddo 2021 2196 enddo … … 2115 2290 wt(il,i)=45.0 2116 2291 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 2117 2300 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 2118 2310 rp(il,i)=rp(il,i+1) 2119 2311 : +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i) 2120 2312 rp(il,i)=0.5*(rp(il,i)+rr(il,i)) 2121 2313 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) 2122 2316 rp(il,i)=max(rp(il,i),0.0) 2123 2317 rp(il,i)=amin1(rp(il,i),rs(il,i)) … … 2126 2320 if(i.eq.1)then 2127 2321 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 2128 2325 else 2129 2326 rp(il,i-1)=rp(il,i) … … 2161 2358 c b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac 2162 2359 c c6 = water(il,i+1) + wdtrain(il)*bfac 2360 c c6 = prec(il,i+1) + wdtrain(il)*bfac 2163 2361 c revap=0.5*(-b6+sqrt(b6*b6+4.*c6)) 2164 2362 c evap(il,i)=sigt*afac*revap 2165 2363 c water(il,i)=revap*revap 2364 c prec(il,i)=revap*revap 2166 2365 cc print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', 2167 2366 cc $ i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) … … 2169 2368 c 2170 2369 c--------retour à la formulation originale d'Emanuel. 2370 if (cvflag_ice) then 2371 2372 c b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac 2373 c c6=prec(il,i+1)+bfac*wdtrain(il) 2374 c : -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1) 2375 c if(c6.gt.0.0)then 2376 c revap=0.5*(-b6+sqrt(b6*b6+4.*c6)) 2377 c 2378 cJAM Attention: evap=sigt*E 2379 c Modification: evap devient l'évaporation en milieu de couche 2380 c car nécessaire dans cv3_yield 2381 c Du coup, il faut modifier pas mal d'équations... 2382 c et l'expression de afac qui devient afac1 2383 c 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) 2387 c print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1 2388 c print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il) 2389 c 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)) 2396 c print*,prec(il,i),'neige' 2397 2398 cjyg Dans sa formulation originale, Emanuel calcule l'evaporation par: 2399 cc evap(il,i)=sigt*afac*revap 2400 c ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee. 2401 c Ici,l'evaporation evap est simplement calculee par l'equation de 2402 c conservation. 2403 c prec(il,i)=revap*revap 2404 c else 2405 cjyg---- Correction : si c6 <= 0, water(il,i)=0. 2406 c prec(il,i)=0. 2407 c endif 2408 2409 cjyg--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i] 2410 c moins [tt ce qui sort de la couche i] 2411 c 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 2438 c water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6 2439 c water(il,i)=max(water(il,i),0.) 2440 c ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6 2441 c ice(il,i)=max(ice(il,i),0.) 2442 c fondue(il,i)=ice(il,i)*thaw 2443 c water(il,i)=water(il,i)+fondue(il,i) 2444 c ice(il,i)=ice(il,i)-fondue(il,i) 2445 c 2446 c if((water(il,i)+ice(il,i)).lt.1.e-30)then 2447 c faci(il,i)=0. 2448 c else 2449 c faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i)) 2450 c endif 2451 2452 else 2171 2453 b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac 2172 2454 c6=water(il,i+1)+bfac*wdtrain(il) … … 2174 2456 if(c6.gt.0.0)then 2175 2457 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*revap2178 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 de2180 c conservation.2181 2458 water(il,i)=revap*revap 2182 2459 else 2183 cjyg---- Correction : si c6 <= 0, water(il,i)=0. 2184 water(il,i) = 0. 2460 water(il,i)= 0. 2185 2461 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] 2462 c print *, 'evap sans ice' 2200 2463 evap(il,i)= 2201 2464 : (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i))) 2202 2465 : /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.) 2203 c 2466 2467 endif 2204 2468 endif !(i.le.inb(il) .and. lwork(il)) 2205 2469 995 Continue … … 2215 2479 tevap(il)=max(0.0,evap(il,i)) 2216 2480 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 2217 2499 if (cvflag_grav) then 2218 2500 mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il) … … 2223 2505 endif 2224 2506 c 2507 endif 2508 c 2225 2509 endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1) 2226 2510 996 Continue … … 2243 2527 : /(lvcp(il,i)*sigd(il)*th(il,i)) 2244 2528 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 2245 2538 bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf 2246 2539 : +50.*(p(il,i-1)-p(il,i))*xf*tevap(il) 2540 endif 2541 2247 2542 fac2=1.0 2248 2543 if(bf.lt.0.0)fac2=-1.0 … … 2262 2557 mp(il,i)=max(0.0,mp(il,i)) 2263 2558 2559 if (cvflag_ice) then 2264 2560 if (cvflag_grav) then 2265 2561 Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante: 2266 2562 C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1). 2267 2563 C 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) 2269 2583 2 /(mp(il,i)+sigd(il)*0.1) 2270 2584 3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i) … … 2276 2590 : *sigd(il)*th(il,i)) 2277 2591 endif 2592 endif 2278 2593 b(il,i-1)=max(b(il,i-1),0.0) 2279 2594 c … … 2391 2706 : ,icb,inb,delt 2392 2707 : ,t,rr,t_wake,rr_wake,s_wake,u,v,tra 2393 : ,gz,p,ph,h,hp,lv, cpn,th,th_wake2708 : ,gz,p,ph,h,hp,lv,lf,cpn,th,th_wake 2394 2709 : ,ep,clw,m,tp,mp,rp,up,vp,trap 2395 : ,wt,water, evap,b,sigd2710 : ,wt,water,ice,evap,fondue,faci,b,sigd 2396 2711 : ,ment,qent,hent,iflag_mix,uent,vent 2397 2712 : ,nent,elij,traent,sig … … 2422 2737 real th(nloc,na), p(nloc,nd), tp(nloc,na) 2423 2738 real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na) 2739 real lf(nloc,na) 2424 2740 real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na) 2425 2741 real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra) 2426 2742 real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc) 2743 real fondue(nloc,na),faci(nloc,na), ice(nloc,na) 2427 2744 real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na) 2428 2745 real hent(nloc,na,na) … … 2455 2772 real cpinv, rdcp, dpinv 2456 2773 real awat(nloc) 2457 real lvcp(nloc,na), mke(nloc,na)2774 real lvcp(nloc,na), lfcp(nloc,na), mke(nloc,na) 2458 2775 real am(nloc), work(nloc), ad(nloc), amp1(nloc) 2459 2776 c!! real up1(nloc), dn1(nloc) … … 2513 2830 do il=1,ncum 2514 2831 lvcp(il,i)=lv(il,i)/cpn(il,i) 2832 lfcp(il,i)=lf(il,i)/cpn(il,i) 2515 2833 enddo 2516 2834 enddo … … 2522 2840 do il=1,ncum 2523 2841 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 2524 2850 if (cvflag_grav) then 2525 2851 precip(il)=wt(il,1)*sigd(il)*water(il,1)*86400.*1000. … … 2527 2853 else 2528 2854 precip(il)=wt(il,1)*sigd(il)*water(il,1)*8640. 2855 endif 2529 2856 endif 2530 2857 endif … … 2539 2866 if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il) 2540 2867 : .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 2541 2877 if (cvflag_grav) then 2542 2878 VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/grav 2543 2879 else 2544 2880 VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/10. 2881 endif 2545 2882 endif 2546 2883 endif … … 2585 2922 cjyg Correction pour conserver l'eau 2586 2923 ccc 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 2588 2932 2589 2933 if (cvflag_grav) then … … 2595 2939 endif 2596 2940 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 2597 2947 ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2) 2598 2948 : *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1) 2949 endif 2599 2950 2600 2951 ftd(il,1) = ft(il,1) ! fin precip … … 2791 3142 ! precip 2792 3143 ccc ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1)) 3144 if (cvflag_ice) then 2793 3145 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 2794 3153 rat=cpn(il,i-1)*cpinv 2795 3154 c … … 2798 3157 : *(mp(il,i+1)*t_wake(il,i)*b(il,i) 2799 3158 : -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv 3159 if (cvflag_ice) then 2800 3160 ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1) 2801 3161 : *(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 2802 3169 ftd(il,i)=ft(il,i) 2803 3170 ! fin precip … … 2818 3185 : *(mp(il,i+1)*t_wake(il,i)*b(il,i) 2819 3186 : -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv 3187 3188 if (cvflag_ice) then 2820 3189 ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1) 2821 3190 : *(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 2822 3198 ftd(il,i)=ft(il,i) 2823 3199 ! fin precip … … 3744 4120 return 3745 4121 end 3746
Note: See TracChangeset
for help on using the changeset viewer.