! ! $Id: cv30_routines.F 1299 2010-01-20 14:27:21Z jyg $ ! c c SUBROUTINE cv30_param(nd,delt) implicit none c------------------------------------------------------------ c Set parameters for convectL for iflag_con = 3 c------------------------------------------------------------ C C *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** C *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** C *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** C *** EFFICIENCY IS ASSUMED TO BE UNITY *** C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** C *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** C *** OF CLOUD *** C C [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA] C *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** C *** APPROACH TO QUASI-EQUILIBRIUM *** C *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** C *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** C C *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** C *** APPROACH TO QUASI-EQUILIBRIUM *** C *** IT MUST BE LESS THAN 0 *** #include "cv30param.h" #include "conema3.h" integer nd real delt ! timestep (seconds) c noff: integer limit for convection (nd-noff) c minorig: First level of convection c -- limit levels for convection: noff = 1 minorig = 1 nl=nd-noff nlp=nl+1 nlm=nl-1 c -- "microphysical" parameters: sigd = 0.01 spfac = 0.15 pbcrit = 150.0 ptcrit = 500.0 cIM cf. FH epmax = 0.993 omtrain = 45.0 ! used also for snow (no disctinction rain/snow) c -- misc: dtovsh = -0.2 ! dT for overshoot dpbase = -40. ! definition cloud base (400m above LCL) dttrig = 5. ! (loose) condition for triggering c -- rate of approach to quasi-equilibrium: dtcrit = -2.0 tau = 8000. beta = 1.0 - delt/tau alpha = 1.5E-3 * delt/tau c increase alpha to compensate W decrease: alpha = alpha*1.5 c -- interface cloud parameterization: delta=0.01 ! cld c -- interface with boundary-layer (gust factor): (sb) betad=10.0 ! original value (from convect 4.3) return end SUBROUTINE cv30_prelim(len,nd,ndp1,t,q,p,ph : ,lv,cpn,tv,gz,h,hm,th) implicit none !===================================================================== ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY ! "ori": from convect4.3 (vectorized) ! "convect3": to be exactly consistent with convect3 !===================================================================== c inputs: integer len, nd, ndp1 real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1) c outputs: real lv(len,nd), cpn(len,nd), tv(len,nd) real gz(len,nd), h(len,nd), hm(len,nd) real th(len,nd) c local variables: integer k, i real rdcp real tvx,tvy ! convect3 real cpx(len,nd) #include "cvthermo.h" #include "cv30param.h" c ori do 110 k=1,nlp do 110 k=1,nl ! convect3 do 100 i=1,len cdebug lv(i,k)= lv0-clmcpv*(t(i,k)-t0) lv(i,k)= lv0-clmcpv*(t(i,k)-273.15) cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k) cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k) c ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) tv(i,k)=t(i,k)*(1.0+q(i,k)/eps-q(i,k)) rdcp=(rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i,k) th(i,k)=t(i,k)*(1000.0/p(i,k))**rdcp 100 continue 110 continue c c gz = phi at the full levels (same as p). c do 120 i=1,len gz(i,1)=0.0 120 continue c ori do 140 k=2,nlp do 140 k=2,nl ! convect3 do 130 i=1,len tvx=t(i,k)*(1.+q(i,k)/eps-q(i,k)) !convect3 tvy=t(i,k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3 gz(i,k)=gz(i,k-1)+0.5*rrd*(tvx+tvy) !convect3 & *(p(i,k-1)-p(i,k))/ph(i,k) !convect3 c ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) c ori & *(p(i,k-1)-p(i,k))/ph(i,k) 130 continue 140 continue c c h = phi + cpT (dry static energy). c hm = phi + cp(T-Tbase)+Lq c c ori do 170 k=1,nlp do 170 k=1,nl ! convect3 do 160 i=1,len h(i,k)=gz(i,k)+cpn(i,k)*t(i,k) hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k) 160 continue 170 continue return end SUBROUTINE cv30_feed(len,nd,t,q,qs,p,ph,hm,gz : ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl) implicit none C================================================================ C Purpose: CONVECTIVE FEED C C Main differences with cv_feed: C - ph added in input C - here, nk(i)=minorig C - icb defined differently (plcl compared with ph instead of p) C C Main differences with convect3: C - we do not compute dplcldt and dplcldr of CLIFT anymore C - values iflag different (but tests identical) C - A,B explicitely defined (!...) C================================================================ #include "cv30param.h" c inputs: integer len, nd real t(len,nd), q(len,nd), qs(len,nd), p(len,nd) real hm(len,nd), gz(len,nd) real ph(len,nd+1) c outputs: integer iflag(len), nk(len), icb(len), icbmax real tnk(len), qnk(len), gznk(len), plcl(len) c local variables: integer i, k integer ihmin(len) real work(len) real pnk(len), qsnk(len), rh(len), chi(len) real A, B ! convect3 cym plcl=0.0 c@ !------------------------------------------------------------------- c@ ! --- Find level of minimum moist static energy c@ ! --- If level of minimum moist static energy coincides with c@ ! --- or is lower than minimum allowable parcel origin level, c@ ! --- set iflag to 6. c@ !------------------------------------------------------------------- c@ c@ do 180 i=1,len c@ work(i)=1.0e12 c@ ihmin(i)=nl c@ 180 continue c@ do 200 k=2,nlp c@ do 190 i=1,len c@ if((hm(i,k).lt.work(i)).and. c@ & (hm(i,k).lt.hm(i,k-1)))then c@ work(i)=hm(i,k) c@ ihmin(i)=k c@ endif c@ 190 continue c@ 200 continue c@ do 210 i=1,len c@ ihmin(i)=min(ihmin(i),nlm) c@ if(ihmin(i).le.minorig)then c@ iflag(i)=6 c@ endif c@ 210 continue c@ c c@ !------------------------------------------------------------------- c@ ! --- Find that model level below the level of minimum moist static c@ ! --- energy that has the maximum value of moist static energy c@ !------------------------------------------------------------------- c@ c@ do 220 i=1,len c@ work(i)=hm(i,minorig) c@ nk(i)=minorig c@ 220 continue c@ do 240 k=minorig+1,nl c@ do 230 i=1,len c@ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then c@ work(i)=hm(i,k) c@ nk(i)=k c@ endif c@ 230 continue c@ 240 continue !------------------------------------------------------------------- ! --- Origin level of ascending parcels for convect3: !------------------------------------------------------------------- do 220 i=1,len nk(i)=minorig 220 continue !------------------------------------------------------------------- ! --- Check whether parcel level temperature and specific humidity ! --- are reasonable !------------------------------------------------------------------- do 250 i=1,len if( ( ( t(i,nk(i)).lt.250.0 ) & .or.( q(i,nk(i)).le.0.0 ) ) c@ & .or.( p(i,ihmin(i)).lt.400.0 ) ) & .and. & ( iflag(i).eq.0) ) iflag(i)=7 250 continue !------------------------------------------------------------------- ! --- Calculate lifted condensation level of air at parcel origin level ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) !------------------------------------------------------------------- A = 1669.0 ! convect3 B = 122.0 ! convect3 do 260 i=1,len if (iflag(i).ne.7) then ! modif sb Jun7th 2002 tnk(i)=t(i,nk(i)) qnk(i)=q(i,nk(i)) gznk(i)=gz(i,nk(i)) pnk(i)=p(i,nk(i)) qsnk(i)=qs(i,nk(i)) c rh(i)=qnk(i)/qsnk(i) c ori rh(i)=min(1.0,rh(i)) ! removed for convect3 c ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) chi(i)=tnk(i)/(A-B*rh(i)-tnk(i)) ! convect3 plcl(i)=pnk(i)*(rh(i)**chi(i)) if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0)) & .and.(iflag(i).eq.0))iflag(i)=8 endif ! iflag=7 260 continue !------------------------------------------------------------------- ! --- Calculate first level above lcl (=icb) !------------------------------------------------------------------- c@ do 270 i=1,len c@ icb(i)=nlm c@ 270 continue c@c c@ do 290 k=minorig,nl c@ do 280 i=1,len c@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i))) c@ & icb(i)=min(icb(i),k) c@ 280 continue c@ 290 continue c@c c@ do 300 i=1,len c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9 c@ 300 continue do 270 i=1,len icb(i)=nlm 270 continue c c la modification consiste a comparer plcl a ph et non a p: c icb est defini par : ph(icb)p(icb), then icbs=icb c c * the routine above computes tvp from minorig to icbs (included). c c * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1) c must be known. This is the case if icbs=icb+1, but not if icbs=icb. c c * therefore, in the case icbs=icb, we compute tvp at level icb+1 c (tvp at other levels will be computed in cv3_undilute2.F) c do i=1,len ticb(i)=t(i,icb(i)+1) gzicb(i)=gz(i,icb(i)+1) qsicb(i)=qs(i,icb(i)+1) enddo do 460 i=1,len tg=ticb(i) qg=qsicb(i) ! convect3 cdebug alv=lv0-clmcpv*(ticb(i)-t0) alv=lv0-clmcpv*(ticb(i)-273.15) c c First iteration. c c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3 : +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 s=1./s c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 tg=tg+s*(ah0(i)-ahg) c ori tg=max(tg,35.0) cdebug tc=tg-t0 tc=tg-273.15 denom=243.5+tc denom=MAX(denom,1.0) ! convect3 c ori if(tc.ge.0.0)then es=6.112*exp(17.67*tc/denom) c ori else c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) c ori endif c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps)) c c Second iteration. c c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) c ori s=1./s c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 tg=tg+s*(ah0(i)-ahg) c ori tg=max(tg,35.0) cdebug tc=tg-t0 tc=tg-273.15 denom=243.5+tc denom=MAX(denom,1.0) ! convect3 c ori if(tc.ge.0.0)then es=6.112*exp(17.67*tc/denom) c ori else c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) c ori end if c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps)) alv=lv0-clmcpv*(ticb(i)-273.15) c ori c approximation here: c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) c ori & -gz(i,icb(i))-alv*qg)/cpd c convect3: no approximation: tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg) : /(cpd+(cl-cpd)*qnk(i)) c ori clw(i,icb(i))=qnk(i)-qg c ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) clw(i,icb(i)+1)=qnk(i)-qg clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1)) rg=qg/(1.-qnk(i)) c ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) c convect3: (qg utilise au lieu du vrai mixing ratio rg) tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing 460 continue return end SUBROUTINE cv30_trigger(len,nd,icb,plcl,p,th,tv,tvp o ,pbase,buoybase,iflag,sig,w0) implicit none !------------------------------------------------------------------- ! --- TRIGGERING ! ! - computes the cloud base ! - triggering (crude in this version) ! - relaxation of sig and w0 when no convection ! ! Caution1: if no convection, we set iflag=4 ! (it used to be 0 in convect3) ! ! Caution2: at this stage, tvp (and thus buoy) are know up ! through icb only! ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy !------------------------------------------------------------------- #include "cv30param.h" c input: integer len, nd integer icb(len) real plcl(len), p(len,nd) real th(len,nd), tv(len,nd), tvp(len,nd) c output: real pbase(len), buoybase(len) c input AND output: integer iflag(len) real sig(len,nd), w0(len,nd) c local variables: integer i,k real tvpbase, tvbase, tdif, ath, ath1 c c *** set cloud base buoyancy at (plcl+dpbase) level buoyancy c do 100 i=1,len pbase(i) = plcl(i) + dpbase tvpbase = tvp(i,icb(i))*(pbase(i)-p(i,icb(i)+1)) : /(p(i,icb(i))-p(i,icb(i)+1)) : + tvp(i,icb(i)+1)*(p(i,icb(i))-pbase(i)) : /(p(i,icb(i))-p(i,icb(i)+1)) tvbase = tv(i,icb(i))*(pbase(i)-p(i,icb(i)+1)) : /(p(i,icb(i))-p(i,icb(i)+1)) : + tv(i,icb(i)+1)*(p(i,icb(i))-pbase(i)) : /(p(i,icb(i))-p(i,icb(i)+1)) buoybase(i) = tvpbase - tvbase 100 continue c c *** make sure that column is dry adiabatic between the surface *** c *** and cloud base, and that lifted air is positively buoyant *** c *** at cloud base *** c *** if not, return to calling program after resetting *** c *** sig(i) and w0(i) *** c c oct3 do 200 i=1,len c oct3 c oct3 tdif = buoybase(i) c oct3 ath1 = th(i,1) c oct3 ath = th(i,icb(i)-1) - dttrig c oct3 c oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then c oct3 do 60 k=1,nl c oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif c oct3 sig(i,k) = AMAX1(sig(i,k),0.0) c oct3 w0(i,k) = beta*w0(i,k) c oct3 60 continue c oct3 iflag(i)=4 ! pour version vectorisee c oct3c convect3 iflag(i)=0 c oct3cccc return c oct3 endif c oct3 c oct3200 continue c -- oct3: on reecrit la boucle 200 (pour la vectorisation) do 60 k=1,nl do 200 i=1,len tdif = buoybase(i) ath1 = th(i,1) ath = th(i,icb(i)-1) - dttrig if (tdif.lt.dtcrit .or. ath.gt.ath1) then sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif sig(i,k) = AMAX1(sig(i,k),0.0) w0(i,k) = beta*w0(i,k) iflag(i)=4 ! pour version vectorisee c convect3 iflag(i)=0 endif 200 continue 60 continue c fin oct3 -- return end SUBROUTINE cv30_compress( len,nloc,ncum,nd,ntra : ,iflag1,nk1,icb1,icbs1 : ,plcl1,tnk1,qnk1,gznk1,pbase1,buoybase1 : ,t1,q1,qs1,u1,v1,gz1,th1 : ,tra1 : ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1 : ,sig1,w01 o ,iflag,nk,icb,icbs o ,plcl,tnk,qnk,gznk,pbase,buoybase o ,t,q,qs,u,v,gz,th o ,tra o ,h,lv,cpn,p,ph,tv,tp,tvp,clw o ,sig,w0 ) implicit none #include "cv30param.h" include 'iniprint.h' c inputs: integer len,ncum,nd,ntra,nloc integer iflag1(len),nk1(len),icb1(len),icbs1(len) real plcl1(len),tnk1(len),qnk1(len),gznk1(len) real pbase1(len),buoybase1(len) real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd) real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd) real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd) real tvp1(len,nd),clw1(len,nd) real th1(len,nd) real sig1(len,nd), w01(len,nd) real tra1(len,nd,ntra) c outputs: c en fait, on a nloc=len pour l'instant (cf cv_driver) integer iflag(nloc),nk(nloc),icb(nloc),icbs(nloc) real plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc) real pbase(nloc),buoybase(nloc) real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd) real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd) real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd) real tvp(nloc,nd),clw(nloc,nd) real th(nloc,nd) real sig(nloc,nd), w0(nloc,nd) real tra(nloc,nd,ntra) c local variables: integer i,k,nn,j CHARACTER (LEN=20) :: modname='cv30_compress' CHARACTER (LEN=80) :: abort_message do 110 k=1,nl+1 nn=0 do 100 i=1,len if(iflag1(i).eq.0)then nn=nn+1 sig(nn,k)=sig1(i,k) w0(nn,k)=w01(i,k) t(nn,k)=t1(i,k) q(nn,k)=q1(i,k) qs(nn,k)=qs1(i,k) u(nn,k)=u1(i,k) v(nn,k)=v1(i,k) gz(nn,k)=gz1(i,k) h(nn,k)=h1(i,k) lv(nn,k)=lv1(i,k) cpn(nn,k)=cpn1(i,k) p(nn,k)=p1(i,k) ph(nn,k)=ph1(i,k) tv(nn,k)=tv1(i,k) tp(nn,k)=tp1(i,k) tvp(nn,k)=tvp1(i,k) clw(nn,k)=clw1(i,k) th(nn,k)=th1(i,k) endif 100 continue 110 continue c do 121 j=1,ntra c do 111 k=1,nd c nn=0 c do 101 i=1,len c if(iflag1(i).eq.0)then c nn=nn+1 c tra(nn,k,j)=tra1(i,k,j) c endif c 101 continue c 111 continue c 121 continue if (nn.ne.ncum) then write(lunout,*)'strange! nn not equal to ncum: ',nn,ncum abort_message = '' CALL abort_gcm (modname,abort_message,1) endif nn=0 do 150 i=1,len if(iflag1(i).eq.0)then nn=nn+1 pbase(nn)=pbase1(i) buoybase(nn)=buoybase1(i) plcl(nn)=plcl1(i) tnk(nn)=tnk1(i) qnk(nn)=qnk1(i) gznk(nn)=gznk1(i) nk(nn)=nk1(i) icb(nn)=icb1(i) icbs(nn)=icbs1(i) iflag(nn)=iflag1(i) endif 150 continue return end SUBROUTINE cv30_undilute2(nloc,ncum,nd,icb,icbs,nk : ,tnk,qnk,gznk,t,q,qs,gz : ,p,h,tv,lv,pbase,buoybase,plcl o ,inb,tp,tvp,clw,hp,ep,sigp,buoy) implicit none C--------------------------------------------------------------------- C Purpose: C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES C & C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD C & C FIND THE LEVEL OF NEUTRAL BUOYANCY C C Main differences convect3/convect4: C - icbs (input) is the first level above LCL (may differ from icb) C - many minor differences in the iterations C - condensed water not removed from tvp in convect3 C - vertical profile of buoyancy computed here (use of buoybase) C - the determination of inb is different C - no inb1, only inb in output C--------------------------------------------------------------------- #include "cvthermo.h" #include "cv30param.h" #include "conema3.h" c inputs: integer ncum, nd, nloc integer icb(nloc), icbs(nloc), nk(nloc) real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd) real p(nloc,nd) real tnk(nloc), qnk(nloc), gznk(nloc) real lv(nloc,nd), tv(nloc,nd), h(nloc,nd) real pbase(nloc), buoybase(nloc), plcl(nloc) c outputs: integer inb(nloc) real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd) real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd) real buoy(nloc,nd) c local variables: integer i, k real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit real by, defrac, pden real ah0(nloc), cape(nloc), capem(nloc), byp(nloc) logical lcape(nloc) !===================================================================== ! --- SOME INITIALIZATIONS !===================================================================== do 170 k=1,nl do 160 i=1,ncum ep(i,k)=0.0 sigp(i,k)=spfac 160 continue 170 continue !===================================================================== ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES !===================================================================== c c --- The procedure is to solve the equation. c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. c c *** Calculate certain parcel quantities, including static energy *** c c do 240 i=1,ncum ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) cdebug & +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i) & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) 240 continue c c c *** Find lifted parcel quantities above cloud base *** c c do 300 k=minorig+1,nl do 290 i=1,ncum c ori if(k.ge.(icb(i)+1))then if(k.ge.(icbs(i)+1))then ! convect3 tg=t(i,k) qg=qs(i,k) cdebug alv=lv0-clmcpv*(t(i,k)-t0) alv=lv0-clmcpv*(t(i,k)-273.15) c c First iteration. c c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3 : +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3 s=1./s c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3 tg=tg+s*(ah0(i)-ahg) c ori tg=max(tg,35.0) cdebug tc=tg-t0 tc=tg-273.15 denom=243.5+tc denom=MAX(denom,1.0) ! convect3 c ori if(tc.ge.0.0)then es=6.112*exp(17.67*tc/denom) c ori else c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) c ori endif qg=eps*es/(p(i,k)-es*(1.-eps)) c c Second iteration. c c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) c ori s=1./s c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3 tg=tg+s*(ah0(i)-ahg) c ori tg=max(tg,35.0) cdebug tc=tg-t0 tc=tg-273.15 denom=243.5+tc denom=MAX(denom,1.0) ! convect3 c ori if(tc.ge.0.0)then es=6.112*exp(17.67*tc/denom) c ori else c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) c ori endif qg=eps*es/(p(i,k)-es*(1.-eps)) c cdebug alv=lv0-clmcpv*(t(i,k)-t0) alv=lv0-clmcpv*(t(i,k)-273.15) c print*,'cpd dans convect2 ',cpd c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd c ori c approximation here: c ori tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd c convect3: no approximation: tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)) clw(i,k)=qnk(i)-qg clw(i,k)=max(0.0,clw(i,k)) rg=qg/(1.-qnk(i)) c ori tvp(i,k)=tp(i,k)*(1.+rg*epsi) c convect3: (qg utilise au lieu du vrai mixing ratio rg): tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing endif 290 continue 300 continue c !===================================================================== ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) !===================================================================== c c ori do 320 k=minorig+1,nl do 320 k=1,nl ! convect3 do 310 i=1,ncum pden=ptcrit-pbcrit ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax ep(i,k)=amax1(ep(i,k),0.0) ep(i,k)=amin1(ep(i,k),epmax) sigp(i,k)=spfac c ori if(k.ge.(nk(i)+1))then c ori tca=tp(i,k)-t0 c ori if(tca.ge.0.0)then c ori elacrit=elcrit c ori else c ori elacrit=elcrit*(1.0-tca/tlcrit) c ori endif c ori elacrit=max(elacrit,0.0) c ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8) c ori ep(i,k)=max(ep(i,k),0.0 ) c ori ep(i,k)=min(ep(i,k),1.0 ) c ori sigp(i,k)=sigs c ori endif 310 continue 320 continue c !===================================================================== ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL ! --- VIRTUAL TEMPERATURE !===================================================================== c c dans convect3, tvp est calcule en une seule fois, et sans retirer c l'eau condensee (~> reversible CAPE) c c ori do 340 k=minorig+1,nl c ori do 330 i=1,ncum c ori if(k.ge.(icb(i)+1))then c ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) c oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' c oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) c ori endif c ori 330 continue c ori 340 continue c ori do 350 i=1,ncum c ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd c ori 350 continue do 350 i=1,ncum ! convect3 tp(i,nlp)=tp(i,nl) ! convect3 350 continue ! convect3 c c===================================================================== c --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only): c===================================================================== c-- this is for convect3 only: c first estimate of buoyancy: do 500 i=1,ncum do 501 k=1,nl buoy(i,k)=tvp(i,k)-tv(i,k) 501 continue 500 continue c set buoyancy=buoybase for all levels below base c for safety, set buoy(icb)=buoybase do 505 i=1,ncum do 506 k=1,nl if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then buoy(i,k)=buoybase(i) endif 506 continue cIM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i) buoy(i,icb(i))=buoybase(i) 505 continue c-- end convect3 c===================================================================== c --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S c --- LEVEL OF NEUTRAL BUOYANCY c===================================================================== c c-- this is for convect3 only: do 510 i=1,ncum inb(i)=nl-1 510 continue do 530 i=1,ncum do 535 k=1,nl-1 if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then inb(i)=MIN(inb(i),k) endif 535 continue 530 continue c-- end convect3 c ori do 510 i=1,ncum c ori cape(i)=0.0 c ori capem(i)=0.0 c ori inb(i)=icb(i)+1 c ori inb1(i)=inb(i) c ori 510 continue c c Originial Code c c do 530 k=minorig+1,nl-1 c do 520 i=1,ncum c if(k.ge.(icb(i)+1))then c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) c cape(i)=cape(i)+by c if(by.ge.0.0)inb1(i)=k+1 c if(cape(i).gt.0.0)then c inb(i)=k+1 c capem(i)=cape(i) c endif c endif c520 continue c530 continue c do 540 i=1,ncum c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) c cape(i)=capem(i)+byp c defrac=capem(i)-cape(i) c defrac=max(defrac,0.001) c frac(i)=-cape(i)/defrac c frac(i)=min(frac(i),1.0) c frac(i)=max(frac(i),0.0) c540 continue c c K Emanuel fix c c call zilch(byp,ncum) c do 530 k=minorig+1,nl-1 c do 520 i=1,ncum c if(k.ge.(icb(i)+1))then c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) c cape(i)=cape(i)+by c if(by.ge.0.0)inb1(i)=k+1 c if(cape(i).gt.0.0)then c inb(i)=k+1 c capem(i)=cape(i) c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) c endif c endif c520 continue c530 continue c do 540 i=1,ncum c inb(i)=max(inb(i),inb1(i)) c cape(i)=capem(i)+byp(i) c defrac=capem(i)-cape(i) c defrac=max(defrac,0.001) c frac(i)=-cape(i)/defrac c frac(i)=min(frac(i),1.0) c frac(i)=max(frac(i),0.0) c540 continue c c J Teixeira fix c c ori call zilch(byp,ncum) c ori do 515 i=1,ncum c ori lcape(i)=.true. c ori 515 continue c ori do 530 k=minorig+1,nl-1 c ori do 520 i=1,ncum c ori if(cape(i).lt.0.0)lcape(i)=.false. c ori if((k.ge.(icb(i)+1)).and.lcape(i))then c ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) c ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) c ori cape(i)=cape(i)+by c ori if(by.ge.0.0)inb1(i)=k+1 c ori if(cape(i).gt.0.0)then c ori inb(i)=k+1 c ori capem(i)=cape(i) c ori endif c ori endif c ori 520 continue c ori 530 continue c ori do 540 i=1,ncum c ori cape(i)=capem(i)+byp(i) c ori defrac=capem(i)-cape(i) c ori defrac=max(defrac,0.001) c ori frac(i)=-cape(i)/defrac c ori frac(i)=min(frac(i),1.0) c ori frac(i)=max(frac(i),0.0) c ori 540 continue c c===================================================================== c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL c===================================================================== c cym do i=1,ncum*nlp cym hp(i,1)=h(i,1) cym enddo do k=1,nlp do i=1,ncum hp(i,k)=h(i,k) enddo enddo do 600 k=minorig+1,nl do 590 i=1,ncum if((k.ge.icb(i)).and.(k.le.inb(i)))then hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k) endif 590 continue 600 continue return end SUBROUTINE cv30_closure(nloc,ncum,nd,icb,inb : ,pbase,p,ph,tv,buoy o ,sig,w0,cape,m) implicit none !=================================================================== ! --- CLOSURE OF CONVECT3 ! ! vectorization: S. Bony !=================================================================== #include "cvthermo.h" #include "cv30param.h" c input: integer ncum, nd, nloc integer icb(nloc), inb(nloc) real pbase(nloc) real p(nloc,nd), ph(nloc,nd+1) real tv(nloc,nd), buoy(nloc,nd) c input/output: real sig(nloc,nd), w0(nloc,nd) c output: real cape(nloc) real m(nloc,nd) c local variables: integer i, j, k, icbmax real deltap, fac, w, amu real dtmin(nloc,nd), sigold(nloc,nd) c ------------------------------------------------------- c -- Initialization c ------------------------------------------------------- do k=1,nl do i=1,ncum m(i,k)=0.0 enddo enddo c ------------------------------------------------------- c -- Reset sig(i) and w0(i) for i>inb and i qnk(il) ! - vectorisation de la partie normalisation des flux (do 789...) !--------------------------------------------------------------------- #include "cvthermo.h" #include "cv30param.h" c inputs: integer ncum, nd, na, ntra, nloc integer icb(nloc), inb(nloc), nk(nloc) real sig(nloc,nd) real qnk(nloc) real ph(nloc,nd+1) real t(nloc,nd), rr(nloc,nd), rs(nloc,nd) real u(nloc,nd), v(nloc,nd) real tra(nloc,nd,ntra) ! input of convect3 real lv(nloc,na), h(nloc,na), hp(nloc,na) real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na) real m(nloc,na) ! input of convect3 c outputs: real ment(nloc,na,na), qent(nloc,na,na) real uent(nloc,na,na), vent(nloc,na,na) real sij(nloc,na,na), elij(nloc,na,na) real traent(nloc,nd,nd,ntra) real ments(nloc,nd,nd), qents(nloc,nd,nd) real sigij(nloc,nd,nd) c local variables: integer i, j, k, il, im, jm integer num1, num2 integer nent(nloc,na) real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp real alt, smid, sjmin, sjmax, delp, delm real asij(nloc), smax(nloc), scrit(nloc) real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd) real wgh real zm(nloc,na) logical lwork(nloc) c===================================================================== c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS c===================================================================== c ori do 360 i=1,ncum*nlp do 361 j=1,nl do 360 i=1,ncum nent(i,j)=0 c in convect3, m is computed in cv3_closure c ori m(i,1)=0.0 360 continue 361 continue c ori do 400 k=1,nlp c ori do 390 j=1,nlp do 400 j=1,nl do 390 k=1,nl do 385 i=1,ncum qent(i,k,j)=rr(i,j) uent(i,k,j)=u(i,j) vent(i,k,j)=v(i,j) elij(i,k,j)=0.0 cym ment(i,k,j)=0.0 cym sij(i,k,j)=0.0 385 continue 390 continue 400 continue cym ment(1:ncum,1:nd,1:nd)=0.0 sij(1:ncum,1:nd,1:nd)=0.0 c do k=1,ntra c do j=1,nd ! instead nlp c do i=1,nd ! instead nlp c do il=1,ncum c traent(il,i,j,k)=tra(il,j,k) c enddo c enddo c enddo c enddo zm(:,:)=0. c===================================================================== c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING c --- FRACTION (sij) c===================================================================== do 750 i=minorig+1, nl do 710 j=minorig,nl do 700 il=1,ncum if( (i.ge.icb(il)).and.(i.le.inb(il)).and. : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then rti=rr(il,1)-ep(il,i)*clw(il,i) bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd) anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j)) denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j) dei=denom if(abs(dei).lt.0.01)dei=0.01 sij(il,i,j)=anum/dei sij(il,i,i)=1.0 altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j) altem=altem/bf2 cwat=clw(il,j)*(1.-ep(il,j)) stemp=sij(il,i,j) if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) : .and.j.gt.i)then anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2) denom=denom+lv(il,j)*(rr(il,i)-rti) if(abs(denom).lt.0.01)denom=0.01 sij(il,i,j)=anum/denom altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j) altem=altem-(bf2-1.)*cwat end if if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il)) vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il)) c!!! do k=1,ntra c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k) c!!! end do elij(il,i,j)=altem elij(il,i,j)=amax1(0.0,elij(il,i,j)) ment(il,i,j)=m(il,i)/(1.-sij(il,i,j)) nent(il,i)=nent(il,i)+1 end if sij(il,i,j)=amax1(0.0,sij(il,i,j)) sij(il,i,j)=amin1(1.0,sij(il,i,j)) endif ! new 700 continue 710 continue c do k=1,ntra c do j=minorig,nl c do il=1,ncum c if( (i.ge.icb(il)).and.(i.le.inb(il)).and. c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) c : +(1.-sij(il,i,j))*tra(il,nk(il),k) c endif c enddo c enddo c enddo c c *** if no air can entrain at level i assume that updraft detrains *** c *** at that level and calculate detrained air flux and properties *** c c@ do 170 i=icb(il),inb(il) do 740 il=1,ncum if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then c@ if(nent(il,i).eq.0)then ment(il,i,i)=m(il,i) qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i) uent(il,i,i)=u(il,nk(il)) vent(il,i,i)=v(il,nk(il)) elij(il,i,i)=clw(il,i) cMAF sij(il,i,i)=1.0 sij(il,i,i)=0.0 end if 740 continue 750 continue c do j=1,ntra c do i=minorig+1,nl c do il=1,ncum c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then c traent(il,i,i,j)=tra(il,nk(il),j) c endif c enddo c enddo c enddo do 100 j=minorig,nl do 101 i=minorig,nl do 102 il=1,ncum if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) : .and.(i.ge.icb(il)).and.(i.le.inb(il)))then sigij(il,i,j)=sij(il,i,j) endif 102 continue 101 continue 100 continue c@ enddo c@170 continue c===================================================================== c --- NORMALIZE ENTRAINED AIR MASS FLUXES c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING c===================================================================== cym call zilch(asum,ncum*nd) cym call zilch(bsum,ncum*nd) cym call zilch(csum,ncum*nd) call zilch(asum,nloc*nd) call zilch(csum,nloc*nd) call zilch(csum,nloc*nd) do il=1,ncum lwork(il) = .FALSE. enddo DO 789 i=minorig+1,nl num1=0 do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1 enddo if (num1.le.0) goto 789 do 781 il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) ) then lwork(il)=(nent(il,i).ne.0) qp=rr(il,1)-ep(il,i)*clw(il,i) anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) : +(cpv-cpd)*t(il,i)*(qp-rr(il,i)) denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) : +(cpd-cpv)*t(il,i)*(rr(il,i)-qp) if(abs(denom).lt.0.01)denom=0.01 scrit(il)=anum/denom alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp) if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0 smax(il)=0.0 asij(il)=0.0 endif 781 continue do 175 j=nl,minorig,-1 num2=0 do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. : j.ge.(icb(il)-1) .and. j.le.inb(il) : .and. lwork(il) ) num2=num2+1 enddo if (num2.le.0) goto 175 do 782 il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. : j.ge.(icb(il)-1) .and. j.le.inb(il) : .and. lwork(il) ) then if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then wgh=1.0 if(j.gt.i)then sjmax=amax1(sij(il,i,j+1),smax(il)) sjmax=amin1(sjmax,scrit(il)) smax(il)=amax1(sij(il,i,j),smax(il)) sjmin=amax1(sij(il,i,j-1),smax(il)) sjmin=amin1(sjmin,scrit(il)) if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0 smid=amin1(sij(il,i,j),scrit(il)) else sjmax=amax1(sij(il,i,j+1),scrit(il)) smid=amax1(sij(il,i,j),scrit(il)) sjmin=0.0 if(j.gt.1)sjmin=sij(il,i,j-1) sjmin=amax1(sjmin,scrit(il)) endif delp=abs(sjmax-smid) delm=abs(sjmin-smid) asij(il)=asij(il)+wgh*(delp+delm) ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh endif endif 782 continue 175 continue do il=1,ncum if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then asij(il)=amax1(1.0e-16,asij(il)) asij(il)=1.0/asij(il) asum(il,i)=0.0 bsum(il,i)=0.0 csum(il,i)=0.0 endif enddo do 180 j=minorig,nl do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then ment(il,i,j)=ment(il,i,j)*asij(il) endif enddo 180 continue do 190 j=minorig,nl do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then asum(il,i)=asum(il,i)+ment(il,i,j) ment(il,i,j)=ment(il,i,j)*sig(il,j) bsum(il,i)=bsum(il,i)+ment(il,i,j) endif enddo 190 continue do il=1,ncum if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then bsum(il,i)=amax1(bsum(il,i),1.0e-16) bsum(il,i)=1.0/bsum(il,i) endif enddo do 195 j=minorig,nl do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i) endif enddo 195 continue do 197 j=minorig,nl do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then csum(il,i)=csum(il,i)+ment(il,i,j) endif enddo 197 continue do il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) : .and. csum(il,i).lt.m(il,i) ) then nent(il,i)=0 ment(il,i,i)=m(il,i) qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i) uent(il,i,i)=u(il,nk(il)) vent(il,i,i)=v(il,nk(il)) elij(il,i,i)=clw(il,i) cMAF sij(il,i,i)=1.0 sij(il,i,i)=0.0 endif enddo ! il c do j=1,ntra c do il=1,ncum c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) c : .and. csum(il,i).lt.m(il,i) ) then c traent(il,i,i,j)=tra(il,nk(il),j) c endif c enddo c enddo 789 continue c c MAF: renormalisation de MENT do jm=1,nd do im=1,nd do il=1,ncum zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm) end do end do end do c do jm=1,nd do im=1,nd do il=1,ncum if(zm(il,im).ne.0.) then ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im) endif end do end do end do c do jm=1,nd do im=1,nd do 999 il=1,ncum qents(il,im,jm)=qent(il,im,jm) ments(il,im,jm)=ment(il,im,jm) 999 continue enddo enddo return end SUBROUTINE cv30_unsat(nloc,ncum,nd,na,ntra,icb,inb : ,t,rr,rs,gz,u,v,tra,p,ph : ,th,tv,lv,cpn,ep,sigp,clw : ,m,ment,elij,delt,plcl : ,mp,rp,up,vp,trap,wt,water,evap,b) implicit none #include "cvthermo.h" #include "cv30param.h" #include "cvflag.h" c inputs: integer ncum, nd, na, ntra, nloc integer icb(nloc), inb(nloc) real delt, plcl(nloc) real t(nloc,nd), rr(nloc,nd), rs(nloc,nd) real u(nloc,nd), v(nloc,nd) real tra(nloc,nd,ntra) real p(nloc,nd), ph(nloc,nd+1) real th(nloc,na), gz(nloc,na) real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na) real cpn(nloc,na), tv(nloc,na) real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na) c outputs: real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na) real water(nloc,na), evap(nloc,na), wt(nloc,na) real trap(nloc,na,ntra) real b(nloc,na) c local variables integer i,j,k,il,num1 real tinv, delti real awat, afac, afac1, afac2, bfac real pr1, pr2, sigt, b6, c6, revap, tevap, delth real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf real ampmax real lvcp(nloc,na) real wdtrain(nloc) logical lwork(nloc) c------------------------------------------------------ delti = 1./delt tinv=1./3. mp(:,:)=0. do i=1,nl do il=1,ncum mp(il,i)=0.0 rp(il,i)=rr(il,i) up(il,i)=u(il,i) vp(il,i)=v(il,i) wt(il,i)=0.001 water(il,i)=0.0 evap(il,i)=0.0 b(il,i)=0.0 lvcp(il,i)=lv(il,i)/cpn(il,i) enddo enddo c do k=1,ntra c do i=1,nd c do il=1,ncum c trap(il,i,k)=tra(il,i,k) c enddo c enddo c enddo c c *** check whether ep(inb)=0, if so, skip precipitating *** c *** downdraft calculation *** c do il=1,ncum lwork(il)=.TRUE. if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE. enddo call zilch(wdtrain,ncum) DO 400 i=nl+1,1,-1 num1=0 do il=1,ncum if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1 enddo if (num1.le.0) goto 400 c c *** integrate liquid water equation to find condensed water *** c *** and condensed water flux *** c c c *** begin downdraft loop *** c c c *** calculate detrained precipitation *** c do il=1,ncum if (i.le.inb(il) .and. lwork(il)) then if (cvflag_grav) then wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i) else wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i) endif endif enddo if(i.gt.1)then do 320 j=1,i-1 do il=1,ncum if (i.le.inb(il) .and. lwork(il)) then awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i) awat=amax1(awat,0.0) if (cvflag_grav) then wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i) else wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i) endif endif enddo 320 continue endif c c *** find rain water and evaporation using provisional *** c *** estimates of rp(i)and rp(i-1) *** c do 999 il=1,ncum if (i.le.inb(il) .and. lwork(il)) then wt(il,i)=45.0 if(i.lt.inb(il))then rp(il,i)=rp(il,i+1) : +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i) rp(il,i)=0.5*(rp(il,i)+rr(il,i)) endif rp(il,i)=amax1(rp(il,i),0.0) rp(il,i)=amin1(rp(il,i),rs(il,i)) rp(il,inb(il))=rr(il,inb(il)) if(i.eq.1)then afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1)) else rp(il,i-1)=rp(il,i) : +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i) rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1)) rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1)) rp(il,i-1)=amax1(rp(il,i-1),0.0) afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i)) afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1)) : /(1.0e4+2000.0*p(il,i-1)*rs(il,i-1)) afac=0.5*(afac1+afac2) endif if(i.eq.inb(il))afac=0.0 afac=amax1(afac,0.0) bfac=1./(sigd*wt(il,i)) c cjyg1 ccc sigt=1.0 ccc if(i.ge.icb)sigt=sigp(i) c prise en compte de la variation progressive de sigt dans c les couches icb et icb-1: c pour plclph(i), pr1=1 & pr2=0 c pour ph(i+1)