/wrfv2_fire/phys/module_cu_gd.F
FORTRAN Legacy | 4555 lines | 3062 code | 317 blank | 1176 comment | 192 complexity | f99ad882e1f460753dbe05b590e319d4 MD5 | raw file
Possible License(s): AGPL-1.0
- !WRF:MODEL_LAYER:PHYSICS
- !
- MODULE module_cu_gd
- CONTAINS
- !-------------------------------------------------------------
- SUBROUTINE GRELLDRV( &
- DT,itimestep,DX &
- ,rho,RAINCV,PRATEC &
- ,U,V,t,W,q,p,pi &
- ,dz8w,p8w,XLV,CP,G,r_v &
- ,STEPCU,htop,hbot &
- ,CU_ACT_FLAG,warm_rain &
- ,APR_GR,APR_W,APR_MC,APR_ST,APR_AS &
- ,APR_CAPMA,APR_CAPME,APR_CAPMI &
- ,MASS_FLUX,XF_ENS,PR_ENS,HT,XLAND,gsw &
- ,GDC,GDC2 &
- ,ensdim,maxiens,maxens,maxens2,maxens3 &
- ,ids,ide, jds,jde, kds,kde &
- ,ims,ime, jms,jme, kms,kme &
- ,its,ite, jts,jte, kts,kte &
- ,periodic_x,periodic_y &
- ,RQVCUTEN,RQCCUTEN,RQICUTEN &
- ,RQVFTEN,RQVBLTEN &
- ,RTHFTEN,RTHCUTEN,RTHRATEN,RTHBLTEN &
- ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
- ,CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,f_flux )
- !-------------------------------------------------------------
- IMPLICIT NONE
- !-------------------------------------------------------------
- INTEGER, INTENT(IN ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- LOGICAL periodic_x,periodic_y
- integer, intent (in ) :: &
- ensdim,maxiens,maxens,maxens2,maxens3
-
- INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP
- LOGICAL, INTENT(IN ) :: warm_rain
- REAL, INTENT(IN ) :: XLV, R_v
- REAL, INTENT(IN ) :: CP,G
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
- INTENT(IN ) :: &
- U, &
- V, &
- W, &
- pi, &
- t, &
- q, &
- p, &
- dz8w, &
- p8w, &
- rho
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
- OPTIONAL , &
- INTENT(INOUT ) :: &
- GDC,GDC2
- REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
- !
- REAL, INTENT(IN ) :: DT, DX
- !
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: RAINCV, PRATEC, MASS_FLUX, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI,htop,hbot
- !+lxz
- ! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
- ! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
- ! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
- ! ! HBOT>HTOP follow physics leveling convention
- LOGICAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: CU_ACT_FLAG
- !
- ! Optionals
- !
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: RTHFTEN, &
- RQVFTEN
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(IN ) :: &
- RTHRATEN, &
- RTHBLTEN, &
- RQVBLTEN
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQCCUTEN, &
- RQICUTEN
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: &
- CFU1, &
- CFD1, &
- DFU1, &
- EFU1, &
- DFD1, &
- EFD1
- !
- ! Flags relating to the optional tendency arrays declared above
- ! Models that carry the optional tendencies will provdide the
- ! optional arguments at compile time; these flags all the model
- ! to determine at run-time whether a particular tracer is in
- ! use or not.
- !
- LOGICAL, OPTIONAL :: &
- F_QV &
- ,F_QC &
- ,F_QR &
- ,F_QI &
- ,F_QS
- LOGICAL, intent(in), OPTIONAL :: f_flux
- ! LOCAL VARS
- real, dimension ( ims:ime , jms:jme , 1:ensdim) :: &
- massfln,xf_ens,pr_ens
- real, dimension (its:ite,kts:kte+1) :: &
- OUTT,OUTQ,OUTQC,phh,cupclw, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
- logical :: l_flux
- real, dimension (its:ite) :: &
- pret, ter11, aa0, fp
- !+lxz
- integer, dimension (its:ite) :: &
- kbcon, ktop
- !.lxz
- integer, dimension (its:ite,jts:jte) :: &
- iact_old_gr
- integer :: ichoice,iens,ibeg,iend,jbeg,jend
- !
- ! basic environmental input includes moisture convergence (mconv)
- ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
- ! convection for this call only and at that particular gridpoint
- !
- real, dimension (its:ite,kts:kte+1) :: &
- T2d,TN,q2d,qo,PO,P2d,US,VS,omeg
- real, dimension (its:ite) :: &
- Z1,PSUR,AAEQ,direction,mconv,cuten,umean,vmean,pmean
- INTEGER :: i,j,k,ICLDCK,ipr,jpr
- REAL :: tcrit,dp,dq
- INTEGER :: itf,jtf,ktf
- REAL :: rkbcon,rktop !-lxz
- l_flux=.FALSE.
- if (present(f_flux)) l_flux=f_flux
- if (l_flux) then
- l_flux = l_flux .and. present(cfu1) .and. present(cfd1) &
- .and. present(dfu1) .and. present(efu1) &
- .and. present(dfd1) .and. present(efd1)
- endif
- ichoice=0
- iens=1
- ipr=0
- jpr=0
- IF ( periodic_x ) THEN
- ibeg=max(its,ids)
- iend=min(ite,ide-1)
- ELSE
- ibeg=max(its,ids+4)
- iend=min(ite,ide-5)
- END IF
- IF ( periodic_y ) THEN
- jbeg=max(jts,jds)
- jend=min(jte,jde-1)
- ELSE
- jbeg=max(jts,jds+4)
- jend=min(jte,jde-5)
- END IF
- tcrit=258.
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- jtf=MIN(jte,jde-1)
- !
- DO 100 J = jts,jtf
- DO I= its,itf
- cuten(i)=0.
- iact_old_gr(i,j)=0
- mass_flux(i,j)=0.
- pratec(i,j) = 0.
- raincv(i,j)=0.
- CU_ACT_FLAG(i,j) = .true.
- ENDDO
- DO k=1,ensdim
- DO I= its,itf
- massfln(i,j,k)=0.
- ENDDO
- ENDDO
- #if ( EM_CORE == 1 )
- DO k= kts,ktf
- DO I= its,itf
- RTHFTEN(i,k,j)=(RTHFTEN(i,k,j)+RTHRATEN(i,k,j)+RTHBLTEN(i,k,j))*pi(i,k,j)
- RQVFTEN(i,k,j)=RQVFTEN(i,k,j)+RQVBLTEN(i,k,j)
- ENDDO
- ENDDO
- #endif
- ! put hydrostatic pressure on half levels
- DO K=kts,ktf
- DO I=ITS,ITF
- phh(i,k) = p(i,k,j)
- ENDDO
- ENDDO
- DO I=ITS,ITF
- PSUR(I)=p8w(I,1,J)*.01
- TER11(I)=HT(i,j)
- mconv(i)=0.
- aaeq(i)=0.
- direction(i)=0.
- pret(i)=0.
- umean(i)=0.
- vmean(i)=0.
- pmean(i)=0.
- ENDDO
- DO K=kts,ktf
- DO I=ITS,ITF
- omeg(i,k)=0.
- ! cupclw(i,k)=0.
- po(i,k)=phh(i,k)*.01
- P2d(I,K)=PO(i,k)
- US(I,K) =u(i,k,j)
- VS(I,K) =v(i,k,j)
- T2d(I,K)=t(i,k,j)
- q2d(I,K)=q(i,k,j)
- omeg(I,K)= -g*rho(i,k,j)*w(i,k,j)
- TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
- IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
- QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
- IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
- IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
- OUTT(I,K)=0.
- OUTQ(I,K)=0.
- OUTQC(I,K)=0.
- ! RTHFTEN(i,k,j)=0.
- ! RQVFTEN(i,k,j)=0.
- ENDDO
- ENDDO
- do k= kts+1,ktf-1
- DO I = its,itf
- if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
- dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
- umean(i)=umean(i)+us(i,k)*dp
- vmean(i)=vmean(i)+vs(i,k)*dp
- pmean(i)=pmean(i)+dp
- endif
- enddo
- enddo
- DO I = its,itf
- if(pmean(i).gt.0)then
- umean(i)=umean(i)/pmean(i)
- vmean(i)=vmean(i)/pmean(i)
- direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
- if(direction(i).gt.360.)direction(i)=direction(i)-360.
- endif
- ENDDO
- DO K=kts,ktf-1
- DO I = its,itf
- dq=(q2d(i,k+1)-q2d(i,k))
- mconv(i)=mconv(i)+omeg(i,k)*dq/g
- ENDDO
- ENDDO
- DO I = its,itf
- if(mconv(i).lt.0.)mconv(i)=0.
- ENDDO
- !
- !---- CALL CUMULUS PARAMETERIZATION
- !
- CALL CUP_enss(outqc,j,AAEQ,T2d,Q2d,TER11,TN,QO,PO,PRET, &
- P2d,OUTT,OUTQ,DT,PSUR,US,VS,tcrit,iens, &
- mconv,massfln,iact_old_gr,omeg,direction,MASS_FLUX, &
- maxiens,maxens,maxens2,maxens3,ensdim, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, &
- xf_ens,pr_ens,XLAND,gsw,cupclw, &
- xlv,r_v,cp,g,ichoice,ipr,jpr, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1,l_flux,&
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- CALL neg_check(dt,q2d,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
- if(j.ge.jbeg.and.j.le.jend)then
- DO I=its,itf
- ! cuten(i)=0.
- if(i.ge.ibeg.and.i.le.iend)then
- if(pret(i).gt.0.)then
- pratec(i,j)=pret(i)
- raincv(i,j)=pret(i)*dt
- cuten(i)=1.
- rkbcon = kte+kts - kbcon(i)
- rktop = kte+kts - ktop(i)
- if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
- if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
- endif
- else
- pret(i)=0.
- endif
- ENDDO
- DO K=kts,ktf
- DO I=its,itf
- RTHCUTEN(I,K,J)=outt(i,k)*cuten(i)/pi(i,k,j)
- RQVCUTEN(I,K,J)=outq(i,k)*cuten(i)
- ENDDO
- ENDDO
- IF(PRESENT(RQCCUTEN)) THEN
- IF ( F_QC ) THEN
- DO K=kts,ktf
- DO I=its,itf
- RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
- IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
- IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
- ENDDO
- ENDDO
- ENDIF
- ENDIF
- !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
- IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
- IF (F_QI) THEN
- DO K=kts,ktf
- DO I=its,itf
- if(t2d(i,k).lt.258.)then
- RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
- RQCCUTEN(I,K,J)=0.
- IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
- else
- RQICUTEN(I,K,J)=0.
- RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
- IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
- endif
- ENDDO
- ENDDO
- ENDIF
- ENDIF
- if (l_flux) then
- DO K=kts,ktf
- DO I=its,itf
- cfu1(i,k,j)=outcfu1(i,k)*cuten(i)
- cfd1(i,k,j)=outcfd1(i,k)*cuten(i)
- dfu1(i,k,j)=outdfu1(i,k)*cuten(i)
- efu1(i,k,j)=outefu1(i,k)*cuten(i)
- dfd1(i,k,j)=outdfd1(i,k)*cuten(i)
- efd1(i,k,j)=outefd1(i,k)*cuten(i)
- enddo
- enddo
- endif
- endif !jbeg,jend
- 100 continue
- END SUBROUTINE GRELLDRV
- SUBROUTINE CUP_enss(OUTQC,J,AAEQ,T,Q,Z1, &
- TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,PSUR,US,VS, &
- TCRIT,iens,mconv,massfln,iact, &
- omeg,direction,massflx,maxiens, &
- maxens,maxens2,maxens3,ensdim, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop, & !-lxz
- xf_ens,pr_ens,xland,gsw,cupclw, &
- xl,rv,cp,g,ichoice,ipr,jpr, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1,l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte,ipr,jpr
- integer, intent (in ) :: &
- j,ensdim,maxiens,maxens,maxens2,maxens3,ichoice,iens
- !
- !
- !
- real, dimension (ims:ime,jms:jme,1:ensdim) &
- ,intent (inout) :: &
- massfln,xf_ens,pr_ens
- real, dimension (ims:ime,jms:jme) &
- ,intent (inout ) :: &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
- APR_CAPME,APR_CAPMI,massflx
- real, dimension (ims:ime,jms:jme) &
- ,intent (in ) :: &
- xland,gsw
- integer, dimension (its:ite,jts:jte) &
- ,intent (in ) :: &
- iact
- ! outtem = output temp tendency (per s)
- ! outq = output q tendency (per s)
- ! outqc = output qc tendency (per s)
- ! pre = output precip
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- OUTT,OUTQ,OUTQC,CUPCLW, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
- logical, intent(in) :: l_flux
- real, dimension (its:ite) &
- ,intent (out ) :: &
- pre
- !+lxz
- integer, dimension (its:ite) &
- ,intent (out ) :: &
- kbcon,ktop
- !.lxz
- !
- ! basic environmental input includes moisture convergence (mconv)
- ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
- ! convection for this call only and at that particular gridpoint
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- T,TN,PO,P,US,VS,omeg
- real, dimension (its:ite,kts:kte) &
- ,intent (inout) :: &
- Q,QO
- real, dimension (its:ite) &
- ,intent (in ) :: &
- Z1,PSUR,AAEQ,direction,mconv
-
- real &
- ,intent (in ) :: &
- dtime,tcrit,xl,cp,rv,g
- !
- ! local ensemble dependent variables in this routine
- !
- real, dimension (its:ite,1:maxens) :: &
- xaa0_ens
- real, dimension (1:maxens) :: &
- mbdt_ens
- real, dimension (1:maxens2) :: &
- edt_ens
- real, dimension (its:ite,1:maxens2) :: &
- edtc
- real, dimension (its:ite,kts:kte,1:maxens2) :: &
- dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens
- real, dimension (its:ite,kts:kte,1:maxens2) :: &
- CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens
- !
- !
- !
- !***************** the following are your basic environmental
- ! variables. They carry a "_cup" if they are
- ! on model cloud levels (staggered). They carry
- ! an "o"-ending (z becomes zo), if they are the forced
- ! variables. They are preceded by x (z becomes xz)
- ! to indicate modification by some typ of cloud
- !
- ! z = heights of model levels
- ! q = environmental mixing ratio
- ! qes = environmental saturation mixing ratio
- ! t = environmental temp
- ! p = environmental pressure
- ! he = environmental moist static energy
- ! hes = environmental saturation moist static energy
- ! z_cup = heights of model cloud levels
- ! q_cup = environmental q on model cloud levels
- ! qes_cup = saturation q on model cloud levels
- ! t_cup = temperature (Kelvin) on model cloud levels
- ! p_cup = environmental pressure
- ! he_cup = moist static energy on model cloud levels
- ! hes_cup = saturation moist static energy on model cloud levels
- ! gamma_cup = gamma on model cloud levels
- !
- !
- ! hcd = moist static energy in downdraft
- ! zd normalized downdraft mass flux
- ! dby = buoancy term
- ! entr = entrainment rate
- ! zd = downdraft normalized mass flux
- ! entr= entrainment rate
- ! hcd = h in model cloud
- ! bu = buoancy term
- ! zd = normalized downdraft mass flux
- ! gamma_cup = gamma on model cloud levels
- ! mentr_rate = entrainment rate
- ! qcd = cloud q (including liquid water) after entrainment
- ! qrch = saturation q in cloud
- ! pwd = evaporate at that level
- ! pwev = total normalized integrated evaoprate (I2)
- ! entr= entrainment rate
- ! z1 = terrain elevation
- ! entr = downdraft entrainment rate
- ! jmin = downdraft originating level
- ! kdet = level above ground where downdraft start detraining
- ! psur = surface pressure
- ! z1 = terrain elevation
- ! pr_ens = precipitation ensemble
- ! xf_ens = mass flux ensembles
- ! massfln = downdraft mass flux ensembles used in next timestep
- ! omeg = omega from large scale model
- ! mconv = moisture convergence from large scale model
- ! zd = downdraft normalized mass flux
- ! zu = updraft normalized mass flux
- ! dir = "storm motion"
- ! mbdt = arbitrary numerical parameter
- ! dtime = dt over which forcing is applied
- ! iact_gr_old = flag to tell where convection was active
- ! kbcon = LFC of parcel from k22
- ! k22 = updraft originating level
- ! icoic = flag if only want one closure (usually set to zero!)
- ! dby = buoancy term
- ! ktop = cloud top (output)
- ! xmb = total base mass flux
- ! hc = cloud moist static energy
- ! hkb = moist static energy at originating level
- ! mentr_rate = entrainment rate
- real, dimension (its:ite,kts:kte) :: &
- he,hes,qes,z, &
- heo,heso,qeso,zo, &
- xhe,xhes,xqes,xz,xt,xq, &
- qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
- qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
- tn_cup, &
- xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
- xt_cup, &
- dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
- dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
- xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
- ! cd = detrainment function for updraft
- ! cdd = detrainment function for downdraft
- ! dellat = change of temperature per unit mass flux of cloud ensemble
- ! dellaq = change of q per unit mass flux of cloud ensemble
- ! dellaqc = change of qc per unit mass flux of cloud ensemble
- cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC, &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
- ! aa0 cloud work function for downdraft
- ! edt = epsilon
- ! aa0 = cloud work function without forcing effects
- ! aa1 = cloud work function with forcing effects
- ! xaa0 = cloud work function with cloud effects (ensemble dependent)
- ! edt = epsilon
- real, dimension (its:ite) :: &
- edt,edto,edtx,AA1,AA0,XAA0,HKB,HKBO,aad,XHKB,QKB,QKBO, &
- XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO,PWEVO,BU,BUO,cap_max,xland1, &
- cap_max_increment,closure_n
- integer, dimension (its:ite) :: &
- kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x, & !-lxz
- KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX
- integer :: &
- nall,iedt,nens,nens3,ki,I,K,KK,iresult
- real :: &
- day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
- zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
- massfld,dh,cap_maxs
- integer :: itf,jtf,ktf
- integer :: jmini
- logical :: keep_going
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- jtf=MIN(jte,jde-1)
- !sms$distribute end
- day=86400.
- do i=its,itf
- closure_n(i)=16.
- xland1(i)=1.
- if(xland(i,j).gt.1.5)xland1(i)=0.
- cap_max_increment(i)=25.
- enddo
- !
- !--- specify entrainmentrate and detrainmentrate
- !
- if(iens.le.4)then
- radius=14000.-float(iens)*2000.
- else
- radius=12000.
- endif
- !
- !--- gross entrainment rate (these may be changed later on in the
- !--- program, depending what your detrainment is!!)
- !
- entr_rate=.2/radius
- !
- !--- entrainment of mass
- !
- mentrd_rate=0.
- mentr_rate=entr_rate
- !
- !--- initial detrainmentrates
- !
- do k=kts,ktf
- do i=its,itf
- cupclw(i,k)=0.
- cd(i,k)=0.1*entr_rate
- cdd(i,k)=0.
- enddo
- enddo
- !
- !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
- ! base mass flux
- !
- edtmax=.8
- edtmin=.2
- !
- !--- minimum depth (m), clouds must have
- !
- depth_min=500.
- !
- !--- maximum depth (mb) of capping
- !--- inversion (larger cap = no convection)
- !
- cap_maxs=75.
- !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
- DO 7 i=its,itf
- kbmax(i)=1
- aa0(i)=0.
- aa1(i)=0.
- aad(i)=0.
- edt(i)=0.
- kstabm(i)=ktf-1
- IERR(i)=0
- IERR2(i)=0
- IERR3(i)=0
- if(aaeq(i).lt.-1.)then
- ierr(i)=20
- endif
- 7 CONTINUE
- !
- !--- first check for upstream convection
- !
- do i=its,itf
- cap_max(i)=cap_maxs
- ! if(tkmax(i,j).lt.2.)cap_max(i)=25.
- if(gsw(i,j).lt.1.)cap_max(i)=25.
- iresult=0
- ! massfld=0.
- ! call cup_direction2(i,j,direction,iact, &
- ! cu_mfx,iresult,0,massfld, &
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte)
- ! cap_max(i)=cap_maxs
- if(iresult.eq.1)then
- cap_max(i)=cap_maxs+20.
- endif
- ! endif
- enddo
- !
- !--- max height(m) above ground where updraft air can originate
- !
- zkbmax=4000.
- !
- !--- height(m) above which no downdrafts are allowed to originate
- !
- zcutdown=3000.
- !
- !--- depth(m) over which downdraft detrains all its mass
- !
- z_detr=1250.
- !
- do nens=1,maxens
- mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
- enddo
- do nens=1,maxens2
- edt_ens(nens)=.95-float(nens)*.01
- enddo
- ! if(j.eq.jpr)then
- ! print *,'radius ensemble ',iens,radius
- ! print *,mbdt_ens
- ! print *,edt_ens
- ! endif
- !
- !--- environmental conditions, FIRST HEIGHTS
- !
- do i=its,itf
- if(ierr(i).ne.20)then
- do k=1,maxens*maxens2*maxens3
- xf_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
- pr_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)=0.
- enddo
- endif
- enddo
- !
- !--- calculate moist static energy, heights, qes
- !
- call cup_env(z,qes,he,hes,t,q,p,z1, &
- psur,ierr,tcrit,0,xl,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
- psur,ierr,tcrit,0,xl,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- environmental values on cloud levels
- !
- call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
- hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
- ierr,z1,xl,rv,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
- heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
- ierr,z1,xl,rv,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- if(ierr(i).eq.0)then
- !
- do k=kts,ktf-2
- if(zo_cup(i,k).gt.zkbmax+z1(i))then
- kbmax(i)=k
- go to 25
- endif
- enddo
- 25 continue
- !
- !
- !--- level where detrainment for downdraft starts
- !
- do k=kts,ktf
- if(zo_cup(i,k).gt.z_detr+z1(i))then
- kdet(i)=k
- go to 26
- endif
- enddo
- 26 continue
- !
- endif
- enddo
- !
- !
- !
- !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
- !
- CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- DO 36 i=its,itf
- IF(ierr(I).eq.0.)THEN
- IF(K22(I).GE.KBMAX(i))ierr(i)=2
- endif
- 36 CONTINUE
- !
- !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
- !
- call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
- ierr,kbmax,po_cup,cap_max, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- ! call cup_kbcon_cin(1,k22,kbcon,heo_cup,heso_cup,z,tn_cup, &
- ! qeso_cup,ierr,kbmax,po_cup,cap_max,xl,cp,&
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte)
- !
- !--- increase detrainment in stable layers
- !
- CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- IF(ierr(I).eq.0.)THEN
- if(kstabm(i)-1.gt.kstabi(i))then
- do k=kstabi(i),kstabm(i)-1
- cd(i,k)=cd(i,k-1)+1.5*entr_rate
- if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
- enddo
- ENDIF
- ENDIF
- ENDDO
- !
- !--- calculate incloud moist static energy
- !
- call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
- kbcon,ierr,dby,he,hes_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
- kbcon,ierr,dbyo,heo,heso_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !--- DETERMINE CLOUD TOP - KTOP
- !
- call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- DO 37 i=its,itf
- kzdown(i)=0
- if(ierr(i).eq.0)then
- zktop=(zo_cup(i,ktop(i))-z1(i))*.6
- zktop=min(zktop+z1(i),zcutdown+z1(i))
- do k=kts,kte
- if(zo_cup(i,k).gt.zktop)then
- kzdown(i)=k
- go to 37
- endif
- enddo
- endif
- 37 CONTINUE
- !
- !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
- !
- call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- DO 100 i=its,ite
- IF(ierr(I).eq.0.)THEN
- !
- !--- check whether it would have buoyancy, if there where
- !--- no entrainment/detrainment
- !
- !jm begin 20061212: the following code replaces code that
- ! was too complex and causing problem for optimization.
- ! Done in consultation with G. Grell.
- jmini = jmin(i)
- keep_going = .TRUE.
- DO WHILE ( keep_going )
- keep_going = .FALSE.
- if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
- if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
- ki = jmini
- hcdo(i,ki)=heso_cup(i,ki)
- DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
- dh=0.
- DO k=ki-1,1,-1
- hcdo(i,k)=heso_cup(i,jmini)
- DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
- dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
- IF(dh.gt.0.)THEN
- jmini=jmini-1
- IF ( jmini .gt. 3 ) THEN
- keep_going = .TRUE.
- ELSE
- ierr(i) = 9
- EXIT
- ENDIF
- ENDIF
- ENDDO
- ENDDO
- jmin(i) = jmini
- IF ( jmini .le. 3 ) THEN
- ierr(i)=4
- ENDIF
- !jm end 20061212
- ENDIF
- 100 CONTINUE
- !
- ! - Must have at least depth_min m between cloud convective base
- ! and cloud top.
- !
- do i=its,itf
- IF(ierr(I).eq.0.)THEN
- IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
- ierr(i)=6
- endif
- endif
- enddo
- !
- !c--- normalized updraft mass flux profile
- !
- call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !c--- normalized downdraft mass flux profile,also work on bottom detrainment
- !--- in this routine
- !
- call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
- 0,kdet,z1, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
- 1,kdet,z1, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- downdraft moist static energy
- !
- call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
- jmin,ierr,he,dbyd,he_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
- jmin,ierr,heo,dbydo,he_cup,&
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- calculate moisture properties of downdraft
- !
- call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
- pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
- pwev,bu,qrcd,q,he,t_cup,2,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_moisture(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
- pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
- pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- calculate moisture properties of updraft
- !
- call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
- kbcon,ktop,cd,dby,mentr_rate,clw_all, &
- q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do k=kts,ktf
- do i=its,itf
- cupclw(i,k)=qrc(i,k)
- enddo
- enddo
- call cup_up_moisture(ierr,zo_cup,qco,qrco,pwo,pwavo, &
- kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
- qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- calculate workfunctions for updrafts
- !
- call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
- kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
- kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- if(ierr(i).eq.0)then
- if(aa1(i).eq.0.)then
- ierr(i)=17
- endif
- endif
- enddo
- !
- !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
- !
- call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
- pwevo,edtmax,edtmin,maxens2,edtc, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do 250 iedt=1,maxens2
- do i=its,itf
- if(ierr(i).eq.0)then
- edt(i)=edtc(i,iedt)
- edto(i)=edtc(i,iedt)
- edtx(i)=edtc(i,iedt)
- endif
- enddo
- do k=kts,ktf
- do i=its,itf
- dellat_ens(i,k,iedt)=0.
- dellaq_ens(i,k,iedt)=0.
- dellaqc_ens(i,k,iedt)=0.
- pwo_ens(i,k,iedt)=0.
- enddo
- enddo
- if (l_flux) then
- do k=kts,ktf
- do i=its,itf
- cfu1_ens(i,k,iedt)=0.
- cfd1_ens(i,k,iedt)=0.
- dfu1_ens(i,k,iedt)=0.
- efu1_ens(i,k,iedt)=0.
- dfd1_ens(i,k,iedt)=0.
- efd1_ens(i,k,iedt)=0.
- enddo
- enddo
- endif
- !
- do i=its,itf
- aad(i)=0.
- enddo
- ! do i=its,itf
- ! if(ierr(i).eq.0)then
- ! eddt(i,j)=edt(i)
- ! EDTX(I)=EDT(I)
- ! BU(I)=0.
- ! BUO(I)=0.
- ! endif
- ! enddo
- !
- !---downdraft workfunctions
- !
- ! call cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
- ! hcd,hes_cup,z,zd, &
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte)
- ! call cup_dd_aa0(edto,ierr,aad,jmin,gammao_cup,tn_cup, &
- ! hcdo,heso_cup,zo,zdo, &
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte)
- !
- !--- change per unit mass that a model cloud would modify the environment
- !
- !--- 1. in bottom layer
- !
- call cup_dellabot(ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto, &
- zuo,zdo,cdd,heo,dellah,j,mentrd_rate,zo,g, &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dellabot(ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto, &
- zuo,zdo,cdd,qo,dellaq,j,mentrd_rate,zo,g,&
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,.FALSE., &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- 2. everywhere else
- !
- call cup_dellas(ierr,zo_cup,po_cup,hcdo,edto,zdo,cdd, &
- heo,dellah,j,mentrd_rate,zuo,g, &
- cd,hco,ktop,k22,kbcon,mentr_rate,jmin,heo_cup,kdet, &
- k22,ipr,jpr,'deep', &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !-- take out cloud liquid water for detrainment
- !
- !?? do k=kts,ktf
- do k=kts,ktf-1
- do i=its,itf
- scr1(i,k)=0.
- dellaqc(i,k)=0.
- if(ierr(i).eq.0)then
- ! print *,'in vupnewg, after della ',ierr(i),aa0(i),i,j
- scr1(i,k)=qco(i,k)-qrco(i,k)
- if(k.eq.ktop(i)-0)dellaqc(i,k)= &
- .01*zuo(i,ktop(i))*qrco(i,ktop(i))* &
- 9.81/(po_cup(i,k)-po_cup(i,k+1))
- if(k.lt.ktop(i).and.k.gt.kbcon(i))then
- dz=zo_cup(i,k+1)-zo_cup(i,k)
- dellaqc(i,k)=.01*9.81*cd(i,k)*dz*zuo(i,k) &
- *.5*(qrco(i,k)+qrco(i,k+1))/ &
- (po_cup(i,k)-po_cup(i,k+1))
- endif
- endif
- enddo
- enddo
- call cup_dellas(ierr,zo_cup,po_cup,qrcdo,edto,zdo,cdd, &
- qo,dellaq,j,mentrd_rate,zuo,g, &
- cd,scr1,ktop,k22,kbcon,mentr_rate,jmin,qo_cup,kdet, &
- k22,ipr,jpr,'deep', &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,.FALSE., &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- !
- !--- using dellas, calculate changed environmental profiles
- !
- ! do 200 nens=1,maxens
- mbdt=mbdt_ens(2)
- do i=its,itf
- xaa0_ens(i,1)=0.
- xaa0_ens(i,2)=0.
- xaa0_ens(i,3)=0.
- enddo
- ! mbdt=mbdt_ens(nens)
- ! do i=its,itf
- ! xaa0_ens(i,nens)=0.
- ! enddo
- do k=kts,ktf
- do i=its,itf
- dellat(i,k)=0.
- if(ierr(i).eq.0)then
- XHE(I,K)=DELLAH(I,K)*MBDT+HEO(I,K)
- XQ(I,K)=DELLAQ(I,K)*MBDT+QO(I,K)
- DELLAT(I,K)=(1./cp)*(DELLAH(I,K)-xl*DELLAQ(I,K))
- XT(I,K)= DELLAT(I,K)*MBDT+TN(I,K)
- IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! print *,k,DELLAH(I,K),DELLAQ(I,K),DELLAT(I,K)
- ! endif
- ENDIF
- enddo
- enddo
- do i=its,itf
- if(ierr(i).eq.0)then
- XHE(I,ktf)=HEO(I,ktf)
- XQ(I,ktf)=QO(I,ktf)
- XT(I,ktf)=TN(I,ktf)
- IF(XQ(I,ktf).LE.0.)XQ(I,ktf)=1.E-08
- endif
- enddo
- !
- !--- calculate moist static energy, heights, qes
- !
- call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
- psur,ierr,tcrit,2,xl,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- environmental values on cloud levels
- !
- call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
- xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
- ierr,z1,xl,rv,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !
- !**************************** static control
- !
- !--- moist static energy inside cloud
- !
- do i=its,itf
- if(ierr(i).eq.0)then
- xhkb(i)=xhe(i,k22(i))
- endif
- enddo
- call cup_up_he(k22,xhkb,xz_cup,cd,mentr_rate,xhe_cup,xhc, &
- kbcon,ierr,xdby,xhe,xhes_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !c--- normalized mass flux profile
- !
- call cup_up_nms(xzu,xz_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- moisture downdraft
- !
- call cup_dd_nms(xzd,xz_cup,cdd,mentrd_rate,jmin,ierr, &
- 1,kdet,z1, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_he(xhes_cup,xzd,xhcd,xz_cup,cdd,mentrd_rate, &
- jmin,ierr,xhe,dbyd,xhe_cup,&
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_moisture(xzd,xhcd,xhes_cup,xqcd,xqes_cup, &
- xpwd,xq_cup,xz_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
- xpwev,bu,xqrcd,xq,xhe,xt_cup,3,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !------- MOISTURE updraft
- !
- call cup_up_moisture(ierr,xz_cup,xqc,xqrc,xpw,xpwav, &
- kbcon,ktop,cd,xdby,mentr_rate,clw_all, &
- xq,GAMMA_cup,xzu,xqes_cup,k22,xq_cup,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- workfunctions for updraft
- !
- call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
- kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- workfunctions for downdraft
- !
- !
- ! call cup_dd_aa0(edtx,ierr,xaa0,jmin,gamma_cup,xt_cup, &
- ! xhcd,xhes_cup,xz,xzd, &
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte)
- do 200 nens=1,maxens
- do i=its,itf
- if(ierr(i).eq.0)then
- xaa0_ens(i,nens)=xaa0(i)
- nall=(iens-1)*maxens3*maxens*maxens2 &
- +(iedt-1)*maxens*maxens3 &
- +(nens-1)*maxens3
- do k=kts,ktf
- if(k.le.ktop(i))then
- do nens3=1,maxens3
- if(nens3.eq.7)then
- !--- b=0
- pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
- pwo(i,k)
- ! +edto(i)*pwdo(i,k)
- !--- b=beta
- else if(nens3.eq.8)then
- pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
- pwo(i,k)
- !--- b=beta/2
- else if(nens3.eq.9)then
- pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
- pwo(i,k)
- ! +.5*edto(i)*pwdo(i,k)
- else
- pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
- pwo(i,k)+edto(i)*pwdo(i,k)
- endif
- enddo
- endif
- enddo
- if(pr_ens(i,j,nall+7).lt.1.e-6)then
- ierr(i)=18
- do nens3=1,maxens3
- pr_ens(i,j,nall+nens3)=0.
- enddo
- endif
- do nens3=1,maxens3
- if(pr_ens(i,j,nall+nens3).lt.1.e-4)then
- pr_ens(i,j,nall+nens3)=0.
- endif
- enddo
- endif
- enddo
- 200 continue
- !
- !--- LARGE SCALE FORCING
- !
- !
- !------- CHECK wether aa0 should have been zero
- !
- !
- CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22x,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- ierr2(i)=ierr(i)
- ierr3(i)=ierr(i)
- enddo
- call cup_kbcon(cap_max_increment,2,k22x,kbconx,heo_cup, &
- heso_cup,ierr2,kbmax,po_cup,cap_max, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- call cup_kbcon(cap_max_increment,3,k22x,kbconx,heo_cup, &
- heso_cup,ierr3,kbmax,po_cup,cap_max, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
- !
- call cup_forcing_ens(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt_ens,dtime, &
- ierr,ierr2,ierr3,xf_ens,j,'deeps', &
- maxens,iens,iedt,maxens2,maxens3,mconv, &
- po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
- massflx,iact,direction,ensdim,massfln,ichoice, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !
- do k=kts,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- dellat_ens(i,k,iedt)=dellat(i,k)
- dellaq_ens(i,k,iedt)=dellaq(i,k)
- dellaqc_ens(i,k,iedt)=dellaqc(i,k)
- pwo_ens(i,k,iedt)=pwo(i,k)+edt(i)*pwdo(i,k)
- else
- dellat_ens(i,k,iedt)=0.
- dellaq_ens(i,k,iedt)=0.
- dellaqc_ens(i,k,iedt)=0.
- pwo_ens(i,k,iedt)=0.
- endif
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! print *,iens,iedt,dellat(i,k),dellat_ens(i,k,iedt), &
- ! dellaq(i,k), dellaqc(i,k)
- ! endif
- enddo
- enddo
- if (l_flux) then
- do k=kts,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- cfu1_ens(i,k,iedt)=cfu1(i,k)
- cfd1_ens(i,k,iedt)=cfd1(i,k)
- dfu1_ens(i,k,iedt)=dfu1(i,k)
- efu1_ens(i,k,iedt)=efu1(i,k)
- dfd1_ens(i,k,iedt)=dfd1(i,k)
- efd1_ens(i,k,iedt)=efd1(i,k)
- else
- cfu1_ens(i,k,iedt)=0.
- cfd1_ens(i,k,iedt)=0.
- dfu1_ens(i,k,iedt)=0.
- efu1_ens(i,k,iedt)=0.
- dfd1_ens(i,k,iedt)=0.
- efd1_ens(i,k,iedt)=0.
- end if
- end do
- end do
- end if
- 250 continue
- !
- !--- FEEDBACK
- !
- call cup_output_ens(xf_ens,ierr,dellat_ens,dellaq_ens, &
- dellaqc_ens,outt,outq,outqc,pre,pwo_ens,xmb,ktop, &
- j,'deep',maxens2,maxens,iens,ierr2,ierr3, &
- pr_ens,maxens3,ensdim,massfln, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1, &
- CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens, &
- l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- PRE(I)=MAX(PRE(I),0.)
- enddo
- !
- !---------------------------done------------------------------
- !
- END SUBROUTINE CUP_enss
- SUBROUTINE cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
- hcd,hes_cup,z,zd, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! aa0 cloud work function for downdraft
- ! gamma_cup = gamma on model cloud levels
- ! t_cup = temperature (Kelvin) on model cloud levels
- ! hes_cup = saturation moist static energy on model cloud levels
- ! hcd = moist static energy in downdraft
- ! edt = epsilon
- ! zd normalized downdraft mass flux
- ! z = heights of model levels
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z,zd,gamma_cup,t_cup,hes_cup,hcd
- real, dimension (its:ite) &
- ,intent (in ) :: &
- edt
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- jmin
- !
- ! input and output
- !
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite) &
- ,intent (out ) :: &
- aa0
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,kk
- real :: &
- dz
- !
- integer :: itf, ktf
- !
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !?? DO k=kts,kte-1
- DO k=kts,ktf-1
- do i=its,itf
- IF(ierr(I).eq.0.and.k.lt.jmin(i))then
- KK=JMIN(I)-K
- !
- !--- ORIGINAL
- !
- DZ=(Z(I,KK)-Z(I,KK+1))
- AA0(I)=AA0(I)+zd(i,kk)*EDT(I)*DZ*(9.81/(1004.*T_cup(I,KK))) &
- *((hcd(i,kk)-hes_cup(i,kk))/(1.+GAMMA_cup(i,kk)))
- endif
- enddo
- enddo
- END SUBROUTINE CUP_dd_aa0
- SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
- pwev,edtmax,edtmin,maxens2,edtc, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- maxens2
- !
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- us,vs,z,p
- real, dimension (its:ite,1:maxens2) &
- ,intent (out ) :: &
- edtc
- real, dimension (its:ite) &
- ,intent (out ) :: &
- edt
- real, dimension (its:ite) &
- ,intent (in ) :: &
- pwav,pwev
- real &
- ,intent (in ) :: &
- edtmax,edtmin
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- ktop,kbcon
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- !
- ! local variables in this routine
- !
- integer i,k,kk
- real einc,pef,pefb,prezk,zkbc
- real, dimension (its:ite) :: &
- vshear,sdp,vws
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
- !
- ! */ calculate an average wind shear over the depth of the cloud
- !
- do i=its,itf
- edt(i)=0.
- vws(i)=0.
- sdp(i)=0.
- vshear(i)=0.
- enddo
- do kk = kts,ktf-1
- do 62 i=its,itf
- IF(ierr(i).ne.0)GO TO 62
- if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
- vws(i) = vws(i)+ &
- (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
- + abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
- (p(i,kk) - p(i,kk+1))
- sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
- endif
- if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
- 62 continue
- end do
- do i=its,itf
- IF(ierr(i).eq.0)then
- pef=(1.591-.639*VSHEAR(I)+.0953*(VSHEAR(I)**2) &
- -.00496*(VSHEAR(I)**3))
- if(pef.gt.edtmax)pef=edtmax
- if(pef.lt.edtmin)pef=edtmin
- !
- !--- cloud base precip efficiency
- !
- zkbc=z(i,kbcon(i))*3.281e-3
- prezk=.02
- if(zkbc.gt.3.)then
- prezk=.96729352+zkbc*(-.70034167+zkbc*(.162179896+zkbc &
- *(- 1.2569798E-2+zkbc*(4.2772E-4-zkbc*5.44E-6))))
- endif
- if(zkbc.gt.25)then
- prezk=2.4
- endif
- pefb=1./(1.+prezk)
- if(pefb.gt.edtmax)pefb=edtmax
- if(pefb.lt.edtmin)pefb=edtmin
- EDT(I)=1.-.5*(pefb+pef)
- !--- edt here is 1-precipeff!
- ! einc=(1.-edt(i))/float(maxens2)
- ! einc=edt(i)/float(maxens2+1)
- !--- 20 percent
- einc=.2*edt(i)
- do k=1,maxens2
- edtc(i,k)=edt(i)+float(k-2)*einc
- enddo
- endif
- enddo
- do i=its,itf
- IF(ierr(i).eq.0)then
- do k=1,maxens2
- EDTC(I,K)=-EDTC(I,K)*PWAV(I)/PWEV(I)
- IF(EDTC(I,K).GT.edtmax)EDTC(I,K)=edtmax
- IF(EDTC(I,K).LT.edtmin)EDTC(I,K)=edtmin
- enddo
- endif
- enddo
- END SUBROUTINE cup_dd_edt
- SUBROUTINE cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,entr, &
- jmin,ierr,he,dby,he_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! hcd = downdraft moist static energy
- ! he = moist static energy on model levels
- ! he_cup = moist static energy on model cloud levels
- ! hes_cup = saturation moist static energy on model cloud levels
- ! dby = buoancy term
- ! cdd= detrainment function
- ! z_cup = heights of model cloud levels
- ! entr = entrainment rate
- ! zd = downdraft normalized mass flux
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- he,he_cup,hes_cup,z_cup,cdd,zd
- ! entr= entrainment rate
- real &
- ,intent (in ) :: &
- entr
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- jmin
- !
- ! input and output
- !
- ! ierr error value, maybe modified in this routine
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- hcd,dby
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,ki
- real :: &
- dz
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- do k=kts+1,ktf
- do i=its,itf
- dby(i,k)=0.
- IF(ierr(I).eq.0)then
- hcd(i,k)=hes_cup(i,k)
- endif
- enddo
- enddo
- !
- do 100 i=its,itf
- IF(ierr(I).eq.0)then
- k=jmin(i)
- hcd(i,k)=hes_cup(i,k)
- dby(i,k)=hcd(i,jmin(i))-hes_cup(i,k)
- !
- do ki=jmin(i)-1,1,-1
- DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
- HCD(i,Ki)=(HCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
- +entr*DZ*HE(i,Ki) &
- )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
- dby(i,ki)=HCD(i,Ki)-hes_cup(i,ki)
- enddo
- !
- endif
- !--- end loop over i
- 100 continue
- END SUBROUTINE cup_dd_he
- SUBROUTINE cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
- pwd,q_cup,z_cup,cdd,entr,jmin,ierr, &
- gamma_cup,pwev,bu,qrcd, &
- q,he,t_cup,iloop,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! cdd= detrainment function
- ! q = environmental q on model levels
- ! q_cup = environmental q on model cloud levels
- ! qes_cup = saturation q on model cloud levels
- ! hes_cup = saturation h on model cloud levels
- ! hcd = h in model cloud
- ! bu = buoancy term
- ! zd = normalized downdraft mass flux
- ! gamma_cup = gamma on model cloud levels
- ! mentr_rate = entrainment rate
- ! qcd = cloud q (including liquid water) after entrainment
- ! qrch = saturation q in cloud
- ! pwd = evaporate at that level
- ! pwev = total normalized integrated evaoprate (I2)
- ! entr= entrainment rate
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- zd,t_cup,hes_cup,hcd,qes_cup,q_cup,z_cup,cdd,gamma_cup,q,he
- real &
- ,intent (in ) :: &
- entr,xl
- integer &
- ,intent (in ) :: &
- iloop
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- jmin
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- qcd,qrcd,pwd
- real, dimension (its:ite) &
- ,intent (out ) :: &
- pwev,bu
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,ki
- real :: &
- dh,dz,dqeva
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- do i=its,itf
- bu(i)=0.
- pwev(i)=0.
- enddo
- do k=kts,ktf
- do i=its,itf
- qcd(i,k)=0.
- qrcd(i,k)=0.
- pwd(i,k)=0.
- enddo
- enddo
- !
- !
- !
- do 100 i=its,itf
- IF(ierr(I).eq.0)then
- k=jmin(i)
- DZ=Z_cup(i,K+1)-Z_cup(i,K)
- qcd(i,k)=q_cup(i,k)
- ! qcd(i,k)=.5*(qes_cup(i,k)+q_cup(i,k))
- qrcd(i,k)=qes_cup(i,k)
- pwd(i,jmin(i))=min(0.,qcd(i,k)-qrcd(i,k))
- pwev(i)=pwev(i)+pwd(i,jmin(i))
- qcd(i,k)=qes_cup(i,k)
- !
- DH=HCD(I,k)-HES_cup(I,K)
- bu(i)=dz*dh
- do ki=jmin(i)-1,1,-1
- DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
- QCD(i,Ki)=(qCD(i,Ki+1)*(1.-.5*CDD(i,Ki)*DZ) &
- +entr*DZ*q(i,Ki) &
- )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
- !
- !--- to be negatively buoyant, hcd should be smaller than hes!
- !
- DH=HCD(I,ki)-HES_cup(I,Ki)
- bu(i)=bu(i)+dz*dh
- QRCD(I,Ki)=qes_cup(i,ki)+(1./XL)*(GAMMA_cup(i,ki) &
- /(1.+GAMMA_cup(i,ki)))*DH
- dqeva=qcd(i,ki)-qrcd(i,ki)
- if(dqeva.gt.0.)dqeva=0.
- pwd(i,ki)=zd(i,ki)*dqeva
- qcd(i,ki)=qrcd(i,ki)
- pwev(i)=pwev(i)+pwd(i,ki)
- ! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
- ! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
- ! endif
- enddo
- !
- !--- end loop over i
- if(pwev(I).eq.0.and.iloop.eq.1)then
- ! print *,'problem with buoy in cup_dd_moisture',i
- ierr(i)=7
- endif
- if(BU(I).GE.0.and.iloop.eq.1)then
- ! print *,'problem with buoy in cup_dd_moisture',i
- ierr(i)=7
- endif
- endif
- 100 continue
- END SUBROUTINE cup_dd_moisture
- SUBROUTINE cup_dd_nms(zd,z_cup,cdd,entr,jmin,ierr, &
- itest,kdet,z1, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! z_cup = height of cloud model level
- ! z1 = terrain elevation
- ! entr = downdraft entrainment rate
- ! jmin = downdraft originating level
- ! kdet = level above ground where downdraft start detraining
- ! itest = flag to whether to calculate cdd
-
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z_cup
- real, dimension (its:ite) &
- ,intent (in ) :: &
- z1
- real &
- ,intent (in ) :: &
- entr
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- jmin,kdet
- integer &
- ,intent (in ) :: &
- itest
- !
- ! input and output
- !
- ! ierr error value, maybe modified in this routine
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- ! zd is the normalized downdraft mass flux
- ! cdd is the downdraft detrainmen function
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- zd
- real, dimension (its:ite,kts:kte) &
- ,intent (inout) :: &
- cdd
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,ki
- real :: &
- a,perc,dz
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !--- perc is the percentage of mass left when hitting the ground
- !
- perc=.03
- do k=kts,ktf
- do i=its,itf
- zd(i,k)=0.
- if(itest.eq.0)cdd(i,k)=0.
- enddo
- enddo
- a=1.-perc
- !
- !
- !
- do 100 i=its,itf
- IF(ierr(I).eq.0)then
- zd(i,jmin(i))=1.
- !
- !--- integrate downward, specify detrainment(cdd)!
- !
- do ki=jmin(i)-1,1,-1
- DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
- if(ki.le.kdet(i).and.itest.eq.0)then
- cdd(i,ki)=entr+(1.- (a*(z_cup(i,ki)-z1(i)) &
- +perc*(z_cup(i,kdet(i))-z1(i)) ) &
- /(a*(z_cup(i,ki+1)-z1(i)) &
- +perc*(z_cup(i,kdet(i))-z1(i))))/dz
- endif
- zd(i,ki)=zd(i,ki+1)*(1.+(entr-cdd(i,ki))*dz)
- enddo
- !
- endif
- !--- end loop over i
- 100 continue
- END SUBROUTINE cup_dd_nms
- SUBROUTINE cup_dellabot(ipr,jpr,he_cup,ierr,z_cup,p_cup, &
- hcd,edt,zu,zd,cdd,he,della,j,mentrd_rate,z,g, &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- j,ipr,jpr
- !
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- della
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z_cup,p_cup,hcd,zu,zd,cdd,he,z,he_cup
- real, dimension (its:ite) &
- ,intent (in ) :: &
- edt
- real &
- ,intent (in ) :: &
- g,mentrd_rate
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite,kts:kte) &
- ,intent (inout ) :: &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
- logical, intent(in) :: l_flux
- !
- ! local variables in this routine
- !
- integer i
- real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
- totmas
- !
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !
- ! if(j.eq.jpr)print *,'in cup dellabot '
- do 100 i=its,itf
- if (l_flux) then
- cfu1(i,1)=0.
- cfd1(i,1)=0.
- cfu1(i,2)=0.
- cfd1(i,2)=0.
- dfu1(i,1)=0.
- efu1(i,1)=0.
- dfd1(i,1)=0.
- efd1(i,1)=0.
- endif
- della(i,1)=0.
- if(ierr(i).ne.0)go to 100
- dz=z_cup(i,2)-z_cup(i,1)
- DP=100.*(p_cup(i,1)-P_cup(i,2))
- detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
- detdo2=edt(i)*zd(i,1)
- entdo=edt(i)*zd(i,2)*mentrd_rate*dz
- subin=-EDT(I)*zd(i,2)
- detdo=detdo1+detdo2-entdo+subin
- DELLA(I,1)=(detdo1*.5*(HCD(i,1)+HCD(i,2)) &
- +detdo2*hcd(i,1) &
- +subin*he_cup(i,2) &
- -entdo*he(i,1))*g/dp
- if (l_flux) then
- cfd1(i,2) = -edt(i)*zd(i,2) !only contribution to subin, subdown=0
- dfd1(i,1) = detdo1+detdo2
- efd1(i,1) = -entdo
- endif
- 100 CONTINUE
- END SUBROUTINE cup_dellabot
- SUBROUTINE cup_dellas(ierr,z_cup,p_cup,hcd,edt,zd,cdd, &
- he,della,j,mentrd_rate,zu,g, &
- cd,hc,ktop,k22,kbcon,mentr_rate,jmin,he_cup,kdet,kpbl, &
- ipr,jpr,name, &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1,l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- j,ipr,jpr
- !
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- della
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z_cup,p_cup,hcd,zd,cdd,he,hc,cd,zu,he_cup
- real, dimension (its:ite) &
- ,intent (in ) :: &
- edt
- real &
- ,intent (in ) :: &
- g,mentrd_rate,mentr_rate
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon,ktop,k22,jmin,kdet,kpbl
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- character *(*), intent (in) :: &
- name
- real, dimension (its:ite,kts:kte) &
- ,intent (inout ) :: &
- CFU1,CFD1,DFU1,EFU1,DFD1,EFD1
- logical, intent(in) :: l_flux
- !
- ! local variables in this routine
- !
- integer i,k
- real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
- detup,subdown,entdoj,entupk,detupk,totmas
- !
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !
- i=ipr
- ! if(j.eq.jpr)then
- ! print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
- ! print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
- ! endif
- DO K=kts+1,ktf
- do i=its,itf
- della(i,k)=0.
- enddo
- enddo
- if (l_flux) then
- DO K=kts+1,ktf-1
- do i=its,itf
- cfu1(i,k+1)=0.
- cfd1(i,k+1)=0.
- enddo
- enddo
- DO K=kts+1,ktf
- do i=its,itf
- dfu1(i,k)=0.
- efu1(i,k)=0.
- dfd1(i,k)=0.
- efd1(i,k)=0.
- enddo
- enddo
- endif
- !
- DO 100 k=kts+1,ktf-1
- DO 100 i=its,ite
- IF(ierr(i).ne.0)GO TO 100
- IF(K.Gt.KTOP(I))GO TO 100
- !
- !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
- !--- WITH ZD CALCULATIONS IN SOUNDD.
- !
- DZ=Z_cup(I,K+1)-Z_cup(I,K)
- detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
- entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
- subin=zu(i,k+1)-zd(i,k+1)*edt(i)
- entup=0.
- detup=0.
- if(k.ge.kbcon(i).and.k.lt.ktop(i))then
- entup=mentr_rate*dz*zu(i,k)
- detup=CD(i,K+1)*DZ*ZU(i,k)
- endif
- subdown=(zu(i,k)-zd(i,k)*edt(i))
- entdoj=0.
- entupk=0.
- detupk=0.
- !
- if(k.eq.jmin(i))then
- entdoj=edt(i)*zd(i,k)
- endif
- if(k.eq.k22(i)-1)then
- entupk=zu(i,kpbl(i))
- endif
- if(k.gt.kdet(i))then
- detdo=0.
- endif
- if(k.eq.ktop(i)-0)then
- detupk=zu(i,ktop(i))
- subin=0.
- endif
- if(k.lt.kbcon(i))then
- detup=0.
- endif
- if (l_flux) then
- ! z_cup(k+1): zu(k+1), -zd(k+1) ==> subin ==> cf[du]1 (k+1) (full-eta level k+1)
- !
- ! z(k) : detup, detdo, entup, entdo ==> [de]f[du]1 (k) (half-eta level k )
- !
- ! z_cup(k) : zu(k), -zd(k) ==> subdown ==> cf[du]1 (k) (full-eta level k )
- ! Store downdraft/updraft mass fluxes at full eta level k (z_cup(k)) in cf[ud]1(k):
- cfu1(i,k+1) = zu(i,k+1)
- cfd1(i,k+1) = -edt(i)*zd(i,k+1)
- ! Store detrainment/entrainment mass fluxes at half eta level k (z(k)) in [de]f[du]1(k):
- dfu1(i,k) = detup+detupk
- efu1(i,k) = -(entup+entupk)
- dfd1(i,k) = detdo
- efd1(i,k) = -(entdo+entdoj)
- endif
- !C
- !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
- !C
- totmas=subin-subdown+detup-entup-entdo+ &
- detdo-entupk-entdoj+detupk
- ! if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k,
- ! 1 totmas,subin,subdown
- ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
- ! 1 entup,entupk,detupk
- ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
- ! 1 detdo,entdoj
- if(abs(totmas).gt.1.e-6)then
- ! print *,'*********************',i,j,k,totmas,name
- ! print *,kpbl(i),k22(i),kbcon(i),ktop(i)
- !c print *,'updr stuff = ',subin,
- !c 1 subdown,detup,entup,entupk,detupk
- !c print *,'dddr stuff = ',entdo,
- !c 1 detdo,entdoj
- ! call wrf_error_fatal ( 'totmas .gt.1.e-6' )
- endif
- dp=100.*(p_cup(i,k-1)-p_cup(i,k))
- della(i,k)=(subin*he_cup(i,k+1) &
- -subdown*he_cup(i,k) &
- +detup*.5*(HC(i,K+1)+HC(i,K)) &
- +detdo*.5*(HCD(i,K+1)+HCD(i,K)) &
- -entup*he(i,k) &
- -entdo*he(i,k) &
- -entupk*he_cup(i,k22(i)) &
- -entdoj*he_cup(i,jmin(i)) &
- +detupk*hc(i,ktop(i)) &
- )*g/dp
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! print *,k,della(i,k),subin*he_cup(i,k+1),subdown*he_cup(i,k),
- ! 1 detdo*.5*(HCD(i,K+1)+HCD(i,K))
- ! print *,k,detup*.5*(HC(i,K+1)+HC(i,K)),detupk*hc(i,ktop(i)),
- ! 1 entup*he(i,k),entdo*he(i,k)
- ! print *,k,he_cup(i,k+1),he_cup(i,k),entupk*he_cup(i,k)
- ! endif
- 100 CONTINUE
- END SUBROUTINE cup_dellas
- SUBROUTINE cup_direction2(i,j,dir,id,massflx, &
- iresult,imass,massfld, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- i,j,imass
- integer, intent (out ) :: &
- iresult
- !
- ! ierr error value, maybe modified in this routine
- !
- integer, dimension (ims:ime,jms:jme) &
- ,intent (in ) :: &
- id
- real, dimension (ims:ime,jms:jme) &
- ,intent (in ) :: &
- massflx
- real, dimension (its:ite) &
- ,intent (inout) :: &
- dir
- real &
- ,intent (out ) :: &
- massfld
- !
- ! local variables in this routine
- !
- integer k,ia,ja,ib,jb
- real diff
- !
- !
- !
- if(imass.eq.1)then
- massfld=massflx(i,j)
- endif
- iresult=0
- ! return
- diff=22.5
- if(dir(i).lt.22.5)dir(i)=360.+dir(i)
- if(id(i,j).eq.1)iresult=1
- ! ja=max(2,j-1)
- ! ia=max(2,i-1)
- ! jb=min(mjx-1,j+1)
- ! ib=min(mix-1,i+1)
- ja=j-1
- ia=i-1
- jb=j+1
- ib=i+1
- if(dir(i).gt.90.-diff.and.dir(i).le.90.+diff)then
- !--- steering flow from the east
- if(id(ib,j).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ib,j),massflx(i,j))
- endif
- return
- endif
- else if(dir(i).gt.135.-diff.and.dir(i).le.135.+diff)then
- !--- steering flow from the south-east
- if(id(ib,ja).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ib,ja),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the south
- else if(dir(i).gt.180.-diff.and.dir(i).le.180.+diff)then
- if(id(i,ja).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(i,ja),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the south west
- else if(dir(i).gt.225.-diff.and.dir(i).le.225.+diff)then
- if(id(ia,ja).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ia,ja),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the west
- else if(dir(i).gt.270.-diff.and.dir(i).le.270.+diff)then
- if(id(ia,j).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ia,j),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the north-west
- else if(dir(i).gt.305.-diff.and.dir(i).le.305.+diff)then
- if(id(ia,jb).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ia,jb),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the north
- else if(dir(i).gt.360.-diff.and.dir(i).le.360.+diff)then
- if(id(i,jb).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(i,jb),massflx(i,j))
- endif
- return
- endif
- !--- steering flow from the north-east
- else if(dir(i).gt.45.-diff.and.dir(i).le.45.+diff)then
- if(id(ib,jb).eq.1)then
- iresult=1
- if(imass.eq.1)then
- massfld=max(massflx(ib,jb),massflx(i,j))
- endif
- return
- endif
- endif
- END SUBROUTINE cup_direction2
- SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
- psur,ierr,tcrit,itest,xl,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- !
- ! ierr error value, maybe modified in this routine
- ! q = environmental mixing ratio
- ! qes = environmental saturation mixing ratio
- ! t = environmental temp
- ! tv = environmental virtual temp
- ! p = environmental pressure
- ! z = environmental heights
- ! he = environmental moist static energy
- ! hes = environmental saturation moist static energy
- ! psur = surface pressure
- ! z1 = terrain elevation
- !
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- p,t
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- he,hes,qes
- real, dimension (its:ite,kts:kte) &
- ,intent (inout) :: &
- z,q
- real, dimension (its:ite) &
- ,intent (in ) :: &
- psur,z1
- real &
- ,intent (in ) :: &
- xl,cp
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- integer &
- ,intent (in ) :: &
- itest
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,iph
- real, dimension (1:2) :: AE,BE,HT
- real, dimension (its:ite,kts:kte) :: tv
- real :: tcrit,e,tvbar
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- HT(1)=XL/CP
- HT(2)=2.834E6/CP
- BE(1)=.622*HT(1)/.286
- AE(1)=BE(1)/273.+ALOG(610.71)
- BE(2)=.622*HT(2)/.286
- AE(2)=BE(2)/273.+ALOG(610.71)
- ! print *, 'TCRIT = ', tcrit,its,ite
- DO k=kts,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- !Csgb - IPH is for phase, dependent on TCRIT (water or ice)
- IPH=1
- IF(T(I,K).LE.TCRIT)IPH=2
- ! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
- E=EXP(AE(IPH)-BE(IPH)/T(I,K))
- ! print *, 'P, E = ', P(I,K), E
- QES(I,K)=.622*E/(100.*P(I,K)-E)
- IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
- IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
- TV(I,K)=T(I,K)+.608*Q(I,K)*T(I,K)
- endif
- enddo
- enddo
- !
- !--- z's are calculated with changed h's and q's and t's
- !--- if itest=2
- !
- if(itest.ne.2)then
- do i=its,itf
- if(ierr(i).eq.0)then
- Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
- ALOG(PSUR(I)))*287.*TV(I,1)/9.81
- endif
- enddo
- ! --- calculate heights
- DO K=kts+1,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- TVBAR=.5*TV(I,K)+.5*TV(I,K-1)
- Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
- ALOG(P(I,K-1)))*287.*TVBAR/9.81
- endif
- enddo
- enddo
- else
- do k=kts,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- z(i,k)=(he(i,k)-1004.*t(i,k)-2.5e6*q(i,k))/9.81
- z(i,k)=max(1.e-3,z(i,k))
- endif
- enddo
- enddo
- endif
- !
- !--- calculate moist static energy - HE
- ! saturated moist static energy - HES
- !
- DO k=kts,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- if(itest.eq.0)HE(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*Q(I,K)
- HES(I,K)=9.81*Z(I,K)+1004.*T(I,K)+2.5E06*QES(I,K)
- IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
- ! if(i.eq.2)then
- ! print *,k,z(i,k),t(i,k),p(i,k),he(i,k),hes(i,k)
- ! endif
- endif
- enddo
- enddo
- END SUBROUTINE cup_env
- SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
- he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
- ierr,z1,xl,rv,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- !
- ! ierr error value, maybe modified in this routine
- ! q = environmental mixing ratio
- ! q_cup = environmental mixing ratio on cloud levels
- ! qes = environmental saturation mixing ratio
- ! qes_cup = environmental saturation mixing ratio on cloud levels
- ! t = environmental temp
- ! t_cup = environmental temp on cloud levels
- ! p = environmental pressure
- ! p_cup = environmental pressure on cloud levels
- ! z = environmental heights
- ! z_cup = environmental heights on cloud levels
- ! he = environmental moist static energy
- ! he_cup = environmental moist static energy on cloud levels
- ! hes = environmental saturation moist static energy
- ! hes_cup = environmental saturation moist static energy on cloud levels
- ! gamma_cup = gamma on cloud levels
- ! psur = surface pressure
- ! z1 = terrain elevation
- !
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- qes,q,he,hes,z,p,t
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
- real, dimension (its:ite) &
- ,intent (in ) :: &
- psur,z1
- real &
- ,intent (in ) :: &
- xl,rv,cp
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- do k=kts+1,ktf
- do i=its,itf
- if(ierr(i).eq.0)then
- qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
- q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
- hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
- he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
- if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
- z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
- p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
- t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
- gamma_cup(i,k)=(xl/cp)*(xl/(rv*t_cup(i,k) &
- *t_cup(i,k)))*qes_cup(i,k)
- endif
- enddo
- enddo
- do i=its,itf
- if(ierr(i).eq.0)then
- qes_cup(i,1)=qes(i,1)
- q_cup(i,1)=q(i,1)
- hes_cup(i,1)=hes(i,1)
- he_cup(i,1)=he(i,1)
- z_cup(i,1)=.5*(z(i,1)+z1(i))
- p_cup(i,1)=.5*(p(i,1)+psur(i))
- t_cup(i,1)=t(i,1)
- gamma_cup(i,1)=xl/cp*(xl/(rv*t_cup(i,1) &
- *t_cup(i,1)))*qes_cup(i,1)
- endif
- enddo
- END SUBROUTINE cup_env_clev
- SUBROUTINE cup_forcing_ens(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
- xf_ens,j,name,maxens,iens,iedt,maxens2,maxens3,mconv, &
- p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon,massflx, &
- iact_old_gr,dir,ensdim,massfln,icoic, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- j,ensdim,maxens,iens,iedt,maxens2,maxens3
- !
- ! ierr error value, maybe modified in this routine
- ! pr_ens = precipitation ensemble
- ! xf_ens = mass flux ensembles
- ! massfln = downdraft mass flux ensembles used in next timestep
- ! omeg = omega from large scale model
- ! mconv = moisture convergence from large scale model
- ! zd = downdraft normalized mass flux
- ! zu = updraft normalized mass flux
- ! aa0 = cloud work function without forcing effects
- ! aa1 = cloud work function with forcing effects
- ! xaa0 = cloud work function with cloud effects (ensemble dependent)
- ! edt = epsilon
- ! dir = "storm motion"
- ! mbdt = arbitrary numerical parameter
- ! dtime = dt over which forcing is applied
- ! iact_gr_old = flag to tell where convection was active
- ! kbcon = LFC of parcel from k22
- ! k22 = updraft originating level
- ! icoic = flag if only want one closure (usually set to zero!)
- ! name = deep or shallow convection flag
- !
- real, dimension (ims:ime,jms:jme,1:ensdim) &
- ,intent (inout) :: &
- pr_ens
- real, dimension (ims:ime,jms:jme,1:ensdim) &
- ,intent (out ) :: &
- xf_ens,massfln
- real, dimension (ims:ime,jms:jme) &
- ,intent (in ) :: &
- massflx
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- omeg,zd,zu,p_cup
- real, dimension (its:ite,1:maxens) &
- ,intent (in ) :: &
- xaa0
- real, dimension (its:ite) &
- ,intent (in ) :: &
- aa1,edt,dir,mconv,xland
- real, dimension (its:ite) &
- ,intent (inout) :: &
- aa0,closure_n
- real, dimension (1:maxens) &
- ,intent (in ) :: &
- mbdt
- real &
- ,intent (in ) :: &
- dtime
- integer, dimension (its:ite,jts:jte) &
- ,intent (in ) :: &
- iact_old_gr
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- k22,kbcon,ktop
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr,ierr2,ierr3
- integer &
- ,intent (in ) :: &
- icoic
- character *(*), intent (in) :: &
- name
- !
- ! local variables in this routine
- !
- real, dimension (1:maxens3) :: &
- xff_ens3
- real, dimension (1:maxens) :: &
- xk
- integer :: &
- i,k,nall,n,ne,nens,nens3,iresult,iresultd,iresulte,mkxcrt,kclim
- parameter (mkxcrt=15)
- real :: &
- a1,massfld,xff0,xomg,aclim1,aclim2,aclim3,aclim4
- real, dimension(1:mkxcrt) :: &
- pcrit,acrit,acritt
- integer :: itf,nall2
- itf=MIN(ite,ide-1)
- DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
- 350.,300.,250.,200.,150./
- DATA ACRIT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
- .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
- ! GDAS DERIVED ACRIT
- DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
- .743,.813,.886,.947,1.138,1.377,1.896/
- !
- nens=0
- !--- LARGE SCALE FORCING
- !
- DO 100 i=its,itf
- ! if(i.eq.ipr.and.j.eq.jpr)print *,'ierr = ',ierr(i)
- if(name.eq.'deeps'.and.ierr(i).gt.995)then
- ! print *,i,j,ierr(i),aa0(i)
- aa0(i)=0.
- ierr(i)=0
- endif
- IF(ierr(i).eq.0)then
- ! kclim=0
- do k=mkxcrt,1,-1
- if(p_cup(i,ktop(i)).lt.pcrit(k))then
- kclim=k
- go to 9
- endif
- enddo
- if(p_cup(i,ktop(i)).ge.pcrit(1))kclim=1
- 9 continue
- kclim=max(kclim,1)
- k=max(kclim-1,1)
- aclim1=acrit(kclim)*1.e3
- aclim2=acrit(k)*1.e3
- aclim3=acritt(kclim)*1.e3
- aclim4=acritt(k)*1.e3
- ! print *,'p_cup(ktop(i)),kclim,pcrit(kclim)'
- ! print *,p_cup(i,ktop(i)),kclim,pcrit(kclim)
- ! print *,'aclim1,aclim2,aclim3,aclim4'
- ! print *,aclim1,aclim2,aclim3,aclim4
- ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
- ! print *,dtime,name,ierr(i),aa1(i),aa0(i)
- !
- !--- treatment different for this closure
- !
- if(name.eq.'deeps')then
- !
- xff0= (AA1(I)-AA0(I))/DTIME
- xff_ens3(1)=(AA1(I)-AA0(I))/dtime
- xff_ens3(2)=.9*xff_ens3(1)
- xff_ens3(3)=1.1*xff_ens3(1)
- !
- !--- more original Arakawa-Schubert (climatologic value of aa0)
- !
- !
- !--- omeg is in bar/s, mconv done with omeg in Pa/s
- ! more like Brown (1979), or Frank-Cohen (199?)
- !
- xff_ens3(4)=-omeg(i,k22(i))/9.81
- xff_ens3(5)=-omeg(i,kbcon(i))/9.81
- xff_ens3(6)=-omeg(i,1)/9.81
- do k=2,kbcon(i)-1
- xomg=-omeg(i,k)/9.81
- if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
- enddo
- !
- !--- more like Krishnamurti et al.
- !
- xff_ens3(7)=mconv(i)
- xff_ens3(8)=mconv(i)
- xff_ens3(9)=mconv(i)
- !
- !--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
- !
- xff_ens3(10)=AA1(I)/(60.*20.)
- xff_ens3(11)=AA1(I)/(60.*30.)
- xff_ens3(12)=AA1(I)/(60.*40.)
- !
- !--- more original Arakawa-Schubert (climatologic value of aa0)
- !
- xff_ens3(13)=max(0.,(AA1(I)-aclim1)/dtime)
- xff_ens3(14)=max(0.,(AA1(I)-aclim2)/dtime)
- xff_ens3(15)=max(0.,(AA1(I)-aclim3)/dtime)
- xff_ens3(16)=max(0.,(AA1(I)-aclim4)/dtime)
- ! if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
- ! xff_ens3(10)=0.
- ! xff_ens3(11)=0.
- ! xff_ens3(12)=0.
- ! xff_ens3(13)=0.
- ! xff_ens3(14)=0.
- ! xff_ens3(15)=0.
- ! xff_ens3(16)=0.
- ! endif
- do nens=1,maxens
- XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(2)
- if(xk(nens).le.0.and.xk(nens).gt.-1.e-6) &
- xk(nens)=-1.e-6
- if(xk(nens).gt.0.and.xk(nens).lt.1.e-6) &
- xk(nens)=1.e-6
- enddo
- !
- !--- add up all ensembles
- !
- do 350 ne=1,maxens
- !
- !--- for every xk, we have maxens3 xffs
- !--- iens is from outermost ensemble (most expensive!
- !
- !--- iedt (maxens2 belongs to it)
- !--- is from second, next outermost, not so expensive
- !
- !--- so, for every outermost loop, we have maxens*maxens2*3
- !--- ensembles!!! nall would be 0, if everything is on first
- !--- loop index, then ne would start counting, then iedt, then iens....
- !
- iresult=0
- iresultd=0
- iresulte=0
- nall=(iens-1)*maxens3*maxens*maxens2 &
- +(iedt-1)*maxens*maxens3 &
- +(ne-1)*maxens3
- !
- ! over water, enfor!e small cap for some of the closures
- !
- if(xland(i).lt.0.1)then
- if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
- ! - ierr2 - 75 mb cap thickness, ierr3 - 125 cap thickness
- ! - for larger cap, set Grell closure to zero
- xff_ens3(1) =0.
- massfln(i,j,nall+1)=0.
- xff_ens3(2) =0.
- massfln(i,j,nall+2)=0.
- xff_ens3(3) =0.
- massfln(i,j,nall+3)=0.
- closure_n(i)=closure_n(i)-1.
- xff_ens3(7) =0.
- massfln(i,j,nall+7)=0.
- xff_ens3(8) =0.
- massfln(i,j,nall+8)=0.
- xff_ens3(9) =0.
- ! massfln(i,j,nall+9)=0.
- closure_n(i)=closure_n(i)-1.
- endif
- !
- ! also take out some closures in general
- !
- xff_ens3(4) =0.
- massfln(i,j,nall+4)=0.
- xff_ens3(5) =0.
- massfln(i,j,nall+5)=0.
- xff_ens3(6) =0.
- massfln(i,j,nall+6)=0.
- closure_n(i)=closure_n(i)-3.
- xff_ens3(10)=0.
- massfln(i,j,nall+10)=0.
- xff_ens3(11)=0.
- massfln(i,j,nall+11)=0.
- xff_ens3(12)=0.
- massfln(i,j,nall+12)=0.
- if(ne.eq.1)closure_n(i)=closure_n(i)-3
- xff_ens3(13)=0.
- massfln(i,j,nall+13)=0.
- xff_ens3(14)=0.
- massfln(i,j,nall+14)=0.
- xff_ens3(15)=0.
- massfln(i,j,nall+15)=0.
- massfln(i,j,nall+16)=0.
- if(ne.eq.1)closure_n(i)=closure_n(i)-4
- endif
- !
- ! end water treatment
- !
- !--- check for upwind convection
- ! iresult=0
- massfld=0.
- ! call cup_direction2(i,j,dir,iact_old_gr, &
- ! massflx,iresult,1, &
- ! massfld, &
- ! ids,ide, jds,jde, kds,kde, &
- ! ims,ime, jms,jme, kms,kme, &
- ! its,ite, jts,jte, kts,kte )
- ! if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
- ! if(iedt.eq.1.and.ne.eq.1)then
- ! print *,massfld,ne,iedt,iens
- ! print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
- ! endif
- ! print *,i,j,massfld,aa0(i),aa1(i)
- IF(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
- iresulte=max(iresult,iresultd)
- iresulte=1
- if(iresulte.eq.1)then
- !
- !--- special treatment for stability closures
- !
- if(xff0.gt.0.)then
- xf_ens(i,j,nall+1)=max(0.,-xff_ens3(1)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+2)=max(0.,-xff_ens3(2)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+3)=max(0.,-xff_ens3(3)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+14)=max(0.,-xff_ens3(14)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+15)=max(0.,-xff_ens3(15)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+16)=max(0.,-xff_ens3(16)/xk(ne)) &
- +massfld
- else
- xf_ens(i,j,nall+1)=massfld
- xf_ens(i,j,nall+2)=massfld
- xf_ens(i,j,nall+3)=massfld
- xf_ens(i,j,nall+13)=massfld
- xf_ens(i,j,nall+14)=massfld
- xf_ens(i,j,nall+15)=massfld
- xf_ens(i,j,nall+16)=massfld
- endif
- !
- !--- if iresult.eq.1, following independent of xff0
- !
- xf_ens(i,j,nall+4)=max(0.,xff_ens3(4) &
- +massfld)
- xf_ens(i,j,nall+5)=max(0.,xff_ens3(5) &
- +massfld)
- xf_ens(i,j,nall+6)=max(0.,xff_ens3(6) &
- +massfld)
- a1=max(1.e-3,pr_ens(i,j,nall+7))
- xf_ens(i,j,nall+7)=max(0.,xff_ens3(7) &
- /a1)
- a1=max(1.e-3,pr_ens(i,j,nall+8))
- xf_ens(i,j,nall+8)=max(0.,xff_ens3(8) &
- /a1)
- a1=max(1.e-3,pr_ens(i,j,nall+9))
- xf_ens(i,j,nall+9)=max(0.,xff_ens3(9) &
- /a1)
- if(XK(ne).lt.0.)then
- xf_ens(i,j,nall+10)=max(0., &
- -xff_ens3(10)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+11)=max(0., &
- -xff_ens3(11)/xk(ne)) &
- +massfld
- xf_ens(i,j,nall+12)=max(0., &
- -xff_ens3(12)/xk(ne)) &
- +massfld
- else
- xf_ens(i,j,nall+10)=massfld
- xf_ens(i,j,nall+11)=massfld
- xf_ens(i,j,nall+12)=massfld
- endif
- if(icoic.ge.1)then
- closure_n(i)=0.
- xf_ens(i,j,nall+1)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+2)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+3)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+4)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+5)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+6)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+7)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+8)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+9)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+10)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+11)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+12)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+13)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+14)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+15)=xf_ens(i,j,nall+icoic)
- xf_ens(i,j,nall+16)=xf_ens(i,j,nall+icoic)
- endif
- !
- ! replace 13-16 for now with other stab closures
- ! (13 gave problems for mass model)
- !
- ! xf_ens(i,j,nall+13)=xf_ens(i,j,nall+1)
- if(icoic.eq.0)xf_ens(i,j,nall+14)=xf_ens(i,j,nall+13)
- ! xf_ens(i,j,nall+15)=xf_ens(i,j,nall+11)
- ! xf_ens(i,j,nall+16)=xf_ens(i,j,nall+11)
- ! xf_ens(i,j,nall+7)=xf_ens(i,j,nall+4)
- ! xf_ens(i,j,nall+8)=xf_ens(i,j,nall+5)
- ! xf_ens(i,j,nall+9)=xf_ens(i,j,nall+6)
- !
- !--- store new for next time step
- !
- do nens3=1,maxens3
- massfln(i,j,nall+nens3)=edt(i) &
- *xf_ens(i,j,nall+nens3)
- massfln(i,j,nall+nens3)=max(0., &
- massfln(i,j,nall+nens3))
- enddo
- !
- !
- !--- do some more on the caps!!! ne=1 for 175, ne=2 for 100,....
- !
- ! do not care for caps here for closure groups 1 and 5,
- ! they are fine, do not turn them off here
- !
- !
- if(ne.eq.2.and.ierr2(i).gt.0)then
- xf_ens(i,j,nall+1) =0.
- xf_ens(i,j,nall+2) =0.
- xf_ens(i,j,nall+3) =0.
- xf_ens(i,j,nall+4) =0.
- xf_ens(i,j,nall+5) =0.
- xf_ens(i,j,nall+6) =0.
- xf_ens(i,j,nall+7) =0.
- xf_ens(i,j,nall+8) =0.
- xf_ens(i,j,nall+9) =0.
- xf_ens(i,j,nall+10)=0.
- xf_ens(i,j,nall+11)=0.
- xf_ens(i,j,nall+12)=0.
- xf_ens(i,j,nall+13)=0.
- xf_ens(i,j,nall+14)=0.
- xf_ens(i,j,nall+15)=0.
- xf_ens(i,j,nall+16)=0.
- massfln(i,j,nall+1)=0.
- massfln(i,j,nall+2)=0.
- massfln(i,j,nall+3)=0.
- massfln(i,j,nall+4)=0.
- massfln(i,j,nall+5)=0.
- massfln(i,j,nall+6)=0.
- massfln(i,j,nall+7)=0.
- massfln(i,j,nall+8)=0.
- massfln(i,j,nall+9)=0.
- massfln(i,j,nall+10)=0.
- massfln(i,j,nall+11)=0.
- massfln(i,j,nall+12)=0.
- massfln(i,j,nall+13)=0.
- massfln(i,j,nall+14)=0.
- massfln(i,j,nall+15)=0.
- massfln(i,j,nall+16)=0.
- endif
- if(ne.eq.3.and.ierr3(i).gt.0)then
- xf_ens(i,j,nall+1) =0.
- xf_ens(i,j,nall+2) =0.
- xf_ens(i,j,nall+3) =0.
- xf_ens(i,j,nall+4) =0.
- xf_ens(i,j,nall+5) =0.
- xf_ens(i,j,nall+6) =0.
- xf_ens(i,j,nall+7) =0.
- xf_ens(i,j,nall+8) =0.
- xf_ens(i,j,nall+9) =0.
- xf_ens(i,j,nall+10)=0.
- xf_ens(i,j,nall+11)=0.
- xf_ens(i,j,nall+12)=0.
- xf_ens(i,j,nall+13)=0.
- xf_ens(i,j,nall+14)=0.
- xf_ens(i,j,nall+15)=0.
- xf_ens(i,j,nall+16)=0.
- massfln(i,j,nall+1)=0.
- massfln(i,j,nall+2)=0.
- massfln(i,j,nall+3)=0.
- massfln(i,j,nall+4)=0.
- massfln(i,j,nall+5)=0.
- massfln(i,j,nall+6)=0.
- massfln(i,j,nall+7)=0.
- massfln(i,j,nall+8)=0.
- massfln(i,j,nall+9)=0.
- massfln(i,j,nall+10)=0.
- massfln(i,j,nall+11)=0.
- massfln(i,j,nall+12)=0.
- massfln(i,j,nall+13)=0.
- massfln(i,j,nall+14)=0.
- massfln(i,j,nall+15)=0.
- massfln(i,j,nall+16)=0.
- endif
- endif
- 350 continue
- ! ne=1, cap=175
- !
- nall=(iens-1)*maxens3*maxens*maxens2 &
- +(iedt-1)*maxens*maxens3
- ! ne=2, cap=100
- !
- nall2=(iens-1)*maxens3*maxens*maxens2 &
- +(iedt-1)*maxens*maxens3 &
- +(2-1)*maxens3
- xf_ens(i,j,nall+4) = xf_ens(i,j,nall2+4)
- xf_ens(i,j,nall+5) =xf_ens(i,j,nall2+5)
- xf_ens(i,j,nall+6) =xf_ens(i,j,nall2+6)
- xf_ens(i,j,nall+7) =xf_ens(i,j,nall2+7)
- xf_ens(i,j,nall+8) =xf_ens(i,j,nall2+8)
- xf_ens(i,j,nall+9) =xf_ens(i,j,nall2+9)
- xf_ens(i,j,nall+10)=xf_ens(i,j,nall2+10)
- xf_ens(i,j,nall+11)=xf_ens(i,j,nall2+11)
- xf_ens(i,j,nall+12)=xf_ens(i,j,nall2+12)
- go to 100
- endif
- elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
- do n=1,ensdim
- xf_ens(i,j,n)=0.
- massfln(i,j,n)=0.
- enddo
- endif
- 100 continue
- END SUBROUTINE cup_forcing_ens
- SUBROUTINE cup_kbcon(cap_inc,iloop,k22,kbcon,he_cup,hes_cup, &
- ierr,kbmax,p_cup,cap_max, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- !
- !
- !
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- he_cup,hes_cup,p_cup
- real, dimension (its:ite) &
- ,intent (in ) :: &
- cap_max,cap_inc
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbmax
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- kbcon,k22,ierr
- integer &
- ,intent (in ) :: &
- iloop
- !
- ! local variables in this routine
- !
- integer :: &
- i
- real :: &
- pbcdif,plus
- integer :: itf
- itf=MIN(ite,ide-1)
- !
- !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
- !
- DO 27 i=its,itf
- kbcon(i)=1
- IF(ierr(I).ne.0)GO TO 27
- KBCON(I)=K22(I)
- GO TO 32
- 31 CONTINUE
- KBCON(I)=KBCON(I)+1
- IF(KBCON(I).GT.KBMAX(i)+2)THEN
- if(iloop.lt.4)ierr(i)=3
- ! if(iloop.lt.4)ierr(i)=997
- GO TO 27
- ENDIF
- 32 CONTINUE
- IF(HE_cup(I,K22(I)).LT.HES_cup(I,KBCON(I)))GO TO 31
- ! cloud base pressure and max moist static energy pressure
- ! i.e., the depth (in mb) of the layer of negative buoyancy
- if(KBCON(I)-K22(I).eq.1)go to 27
- PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
- plus=max(25.,cap_max(i)-float(iloop-1)*cap_inc(i))
- if(iloop.eq.4)plus=cap_max(i)
- IF(PBCDIF.GT.plus)THEN
- K22(I)=K22(I)+1
- KBCON(I)=K22(I)
- GO TO 32
- ENDIF
- 27 CONTINUE
- END SUBROUTINE cup_kbcon
- SUBROUTINE cup_kbcon_cin(iloop,k22,kbcon,he_cup,hes_cup, &
- z,tmean,qes,ierr,kbmax,p_cup,cap_max,xl,cp, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- !
- !
- !
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- he_cup,hes_cup,p_cup,z,tmean,qes
- real, dimension (its:ite) &
- ,intent (in ) :: &
- cap_max
- real &
- ,intent (in ) :: &
- xl,cp
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbmax
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- kbcon,k22,ierr
- integer &
- ,intent (in ) :: &
- iloop
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- real :: &
- cin,cin_max,dh,tprim,gamma
- !
- integer :: itf
- itf=MIN(ite,ide-1)
- !
- !
-
- !
- !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
- !
- DO 27 i=its,itf
- cin_max=-cap_max(i)
- kbcon(i)=1
- cin = 0.
- IF(ierr(I).ne.0)GO TO 27
- KBCON(I)=K22(I)
- GO TO 32
- 31 CONTINUE
- KBCON(I)=KBCON(I)+1
- IF(KBCON(I).GT.KBMAX(i)+2)THEN
- if(iloop.eq.1)ierr(i)=3
- ! if(iloop.eq.2)ierr(i)=997
- GO TO 27
- ENDIF
- 32 CONTINUE
- dh = HE_cup(I,K22(I)) - HES_cup(I,KBCON(I))
- if (dh.lt. 0.) then
- GAMMA=(xl/cp)*(xl/(461.525*(Tmean(I,K22(i))**2)))*QES(I,K22(i))
- tprim = dh/(cp*(1.+gamma))
- cin = cin + 9.8066 * tprim &
- *(z(i,k22(i))-z(i,k22(i)-1)) / tmean(i,k22(i))
- go to 31
- end if
- ! If negative energy in negatively buoyant layer
- ! exceeds convective inhibition (CIN) threshold,
- ! then set K22 level one level up and see if that
- ! will work.
- IF(cin.lT.cin_max)THEN
- ! print *,i,cin,cin_max,k22(i),kbcon(i)
- K22(I)=K22(I)+1
- KBCON(I)=K22(I)
- GO TO 32
- ENDIF
- 27 CONTINUE
- END SUBROUTINE cup_kbcon_cin
- SUBROUTINE cup_ktop(ilo,dby,kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! dby = buoancy term
- ! ktop = cloud top (output)
- ! ilo = flag
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (inout) :: &
- dby
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon
- integer &
- ,intent (in ) :: &
- ilo
- integer, dimension (its:ite) &
- ,intent (out ) :: &
- ktop
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- !
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- !
- !
- DO 42 i=its,itf
- ktop(i)=1
- IF(ierr(I).EQ.0)then
- DO 40 K=KBCON(I)+1,ktf-1
- IF(DBY(I,K).LE.0.)THEN
- KTOP(I)=K-1
- GO TO 41
- ENDIF
- 40 CONTINUE
- if(ilo.eq.1)ierr(i)=5
- ! if(ilo.eq.2)ierr(i)=998
- GO TO 42
- 41 CONTINUE
- do k=ktop(i)+1,ktf
- dby(i,k)=0.
- enddo
- endif
- 42 CONTINUE
- END SUBROUTINE cup_ktop
- SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! array input array
- ! x output array with return values
- ! kt output array of levels
- ! ks,kend check-range
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- array
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- ierr,ke
- integer &
- ,intent (in ) :: &
- ks
- integer, dimension (its:ite) &
- ,intent (out ) :: &
- maxx
- real, dimension (its:ite) :: &
- x
- real :: &
- xar
- integer :: &
- i,k
- integer :: itf
- itf=MIN(ite,ide-1)
- DO 200 i=its,itf
- MAXX(I)=KS
- if(ierr(i).eq.0)then
- X(I)=ARRAY(I,KS)
- !
- DO 100 K=KS,KE(i)
- XAR=ARRAY(I,K)
- IF(XAR.GE.X(I)) THEN
- X(I)=XAR
- MAXX(I)=K
- ENDIF
- 100 CONTINUE
- endif
- 200 CONTINUE
- END SUBROUTINE cup_MAXIMI
- SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! array input array
- ! x output array with return values
- ! kt output array of levels
- ! ks,kend check-range
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- array
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- ierr,ks,kend
- integer, dimension (its:ite) &
- ,intent (out ) :: &
- kt
- real, dimension (its:ite) :: &
- x
- integer :: &
- i,k,kstop
- integer :: itf
- itf=MIN(ite,ide-1)
- DO 200 i=its,itf
- KT(I)=KS(I)
- if(ierr(i).eq.0)then
- X(I)=ARRAY(I,KS(I))
- KSTOP=MAX(KS(I)+1,KEND(I))
- !
- DO 100 K=KS(I)+1,KSTOP
- IF(ARRAY(I,K).LT.X(I)) THEN
- X(I)=ARRAY(I,K)
- KT(I)=K
- ENDIF
- 100 CONTINUE
- endif
- 200 CONTINUE
- END SUBROUTINE cup_MINIMI
- SUBROUTINE cup_output_ens(xf_ens,ierr,dellat,dellaq,dellaqc, &
- outtem,outq,outqc,pre,pw,xmb,ktop, &
- j,name,nx,nx2,iens,ierr2,ierr3,pr_ens, &
- maxens3,ensdim,massfln, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1, &
- CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens, &
- l_flux, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- integer, intent (in ) :: &
- j,ensdim,nx,nx2,iens,maxens3
- ! xf_ens = ensemble mass fluxes
- ! pr_ens = precipitation ensembles
- ! dellat = change of temperature per unit mass flux of cloud ensemble
- ! dellaq = change of q per unit mass flux of cloud ensemble
- ! dellaqc = change of qc per unit mass flux of cloud ensemble
- ! outtem = output temp tendency (per s)
- ! outq = output q tendency (per s)
- ! outqc = output qc tendency (per s)
- ! pre = output precip
- ! xmb = total base mass flux
- ! xfac1 = correction factor
- ! pw = pw -epsilon*pd (ensemble dependent)
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (ims:ime,jms:jme,1:ensdim) &
- ,intent (inout) :: &
- xf_ens,pr_ens,massfln
- real, dimension (ims:ime,jms:jme) &
- ,intent (inout) :: &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
- APR_CAPME,APR_CAPMI
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- outtem,outq,outqc
- real, dimension (its:ite) &
- ,intent (out ) :: &
- pre,xmb
- real, dimension (its:ite) &
- ,intent (inout ) :: &
- closure_n,xland1
- real, dimension (its:ite,kts:kte,1:nx) &
- ,intent (in ) :: &
- dellat,dellaqc,dellaq,pw
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- ktop
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr,ierr2,ierr3
- real, dimension (its:ite,kts:kte,1:ensdim) &
- ,intent (in ) :: &
- CFU1_ens,CFD1_ens,DFU1_ens,EFU1_ens,DFD1_ens,EFD1_ens
- real, dimension (its:ite,kts:kte) &
- ,intent (out) :: &
- outCFU1,outCFD1,outDFU1,outEFU1,outDFD1,outEFD1
- logical, intent(in) :: l_flux
- !
- ! local variables in this routine
- !
- integer :: &
- i,k,n,ncount
- real :: &
- outtes,ddtes,dtt,dtq,dtqc,dtpw,tuning,prerate,clos_wei
- real, dimension (its:ite) :: &
- xfac1
- real, dimension (its:ite):: &
- xmb_ske,xmb_ave,xmb_std,xmb_cur,xmbweight
- real, dimension (its:ite):: &
- pr_ske,pr_ave,pr_std,pr_cur
- real, dimension (its:ite,jts:jte):: &
- pr_gr,pr_w,pr_mc,pr_st,pr_as,pr_capma, &
- pr_capme,pr_capmi
- integer :: iedt, kk
- !
- character *(*), intent (in) :: &
- name
- !
- integer :: itf, ktf
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- tuning=0.
- !
- !
- DO k=kts,ktf
- do i=its,itf
- outtem(i,k)=0.
- outq(i,k)=0.
- outqc(i,k)=0.
- enddo
- enddo
- do i=its,itf
- pre(i)=0.
- xmb(i)=0.
- xfac1(i)=1.
- xmbweight(i)=1.
- enddo
- do i=its,itf
- IF(ierr(i).eq.0)then
- do n=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
- if(pr_ens(i,j,n).le.0.)then
- xf_ens(i,j,n)=0.
- endif
- enddo
- endif
- enddo
- !
- !--- calculate ensemble average mass fluxes
- !
- call massflx_stats(xf_ens,ensdim,nx2,nx,maxens3, &
- xmb_ave,xmb_std,xmb_cur,xmb_ske,j,ierr,1, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI, &
- pr_gr,pr_w,pr_mc,pr_st,pr_as, &
- pr_capma,pr_capme,pr_capmi, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- call massflx_stats(pr_ens,ensdim,nx2,nx,maxens3, &
- pr_ave,pr_std,pr_cur,pr_ske,j,ierr,2, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI, &
- pr_gr,pr_w,pr_mc,pr_st,pr_as, &
- pr_capma,pr_capme,pr_capmi, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- !
- !-- now do feedback
- !
- ddtes=200.
- ! if(name.eq.'shal')ddtes=200.
- do i=its,itf
- if(ierr(i).eq.0)then
- if(xmb_ave(i).le.0.)then
- ierr(i)=13
- xmb_ave(i)=0.
- endif
- ! xmb(i)=max(0.,xmb_ave(i))
- xmb(i)=max(.1*xmb_ave(i),xmb_ave(i)-tuning*xmb_std(i))
- ! --- Now use proper count of how many closures were actually
- ! used in cup_forcing_ens (including screening of some
- ! closures over water) to properly normalize xmb
- clos_wei=16./max(1.,closure_n(i))
- if (xland1(i).lt.0.5)xmb(i)=xmb(i)*clos_wei
- if(xmb(i).eq.0.)then
- ierr(i)=19
- endif
- if(xmb(i).gt.100.)then
- ierr(i)=19
- endif
- xfac1(i)=xmb(i)
- endif
- xfac1(i)=xmb_ave(i)
- ENDDO
- DO k=kts,ktf
- do i=its,itf
- dtt=0.
- dtq=0.
- dtqc=0.
- dtpw=0.
- IF(ierr(i).eq.0.and.k.le.ktop(i))then
- do n=1,nx
- dtt=dtt+dellat(i,k,n)
- dtq=dtq+dellaq(i,k,n)
- dtqc=dtqc+dellaqc(i,k,n)
- dtpw=dtpw+pw(i,k,n)
- enddo
- outtes=dtt*XMB(I)*86400./float(nx)
- IF((OUTTES.GT.2.*ddtes.and.k.gt.2))THEN
- XMB(I)= 2.*ddtes/outtes * xmb(i)
- outtes=1.*ddtes
- endif
- if (outtes .lt. -ddtes) then
- XMB(I)= -ddtes/outtes * xmb(i)
- outtes=-ddtes
- endif
- if (outtes .gt. .5*ddtes.and.k.le.2) then
- XMB(I)= ddtes/outtes * xmb(i)
- outtes=.5*ddtes
- endif
- OUTTEM(I,K)=XMB(I)*dtt/float(nx)
- OUTQ(I,K)=XMB(I)*dtq/float(nx)
- OUTQC(I,K)=XMB(I)*dtqc/float(nx)
- PRE(I)=PRE(I)+XMB(I)*dtpw/float(nx)
- endif
- enddo
- enddo
- do i=its,itf
- if(ierr(i).eq.0)then
- prerate=pre(i)*3600.
- if(prerate.lt.0.1)then
- if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
- pre(i)=0.
- ierr(i)=221
- do k=kts,ktf
- outtem(i,k)=0.
- outq(i,k)=0.
- outqc(i,k)=0.
- enddo
- do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
- massfln(i,j,k)=0.
- xf_ens(i,j,k)=0.
- enddo
- endif
- endif
- endif
- ENDDO
- do i=its,itf
- if(ierr(i).eq.0)then
- xfac1(i)=xmb(i)/xfac1(i)
- do k=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
- massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
- xf_ens(i,j,k)=xf_ens(i,j,k)*xfac1(i)
- enddo
- endif
- ENDDO
- if (l_flux) then
- if (iens .eq. 1) then ! Only do deep convection mass fluxes
- do k=kts,ktf
- do i=its,itf
- outcfu1(i,k)=0.
- outcfd1(i,k)=0.
- outdfu1(i,k)=0.
- outefu1(i,k)=0.
- outdfd1(i,k)=0.
- outefd1(i,k)=0.
- if (ierr(i) .eq. 0) then
- do iedt=1,nx
- do kk=1,nx2*maxens3
- n=(iens-1)*nx*nx2*maxens3 + &
- (iedt-1)*nx2*maxens3 + kk
- outcfu1(i,k)=outcfu1(i,k)+cfu1_ens(i,k,iedt)*xf_ens(i,j,n)
- outcfd1(i,k)=outcfd1(i,k)+cfd1_ens(i,k,iedt)*xf_ens(i,j,n)
- outdfu1(i,k)=outdfu1(i,k)+dfu1_ens(i,k,iedt)*xf_ens(i,j,n)
- outefu1(i,k)=outefu1(i,k)+efu1_ens(i,k,iedt)*xf_ens(i,j,n)
- outdfd1(i,k)=outdfd1(i,k)+dfd1_ens(i,k,iedt)*xf_ens(i,j,n)
- outefd1(i,k)=outefd1(i,k)+efd1_ens(i,k,iedt)*xf_ens(i,j,n)
- end do
- end do
- outcfu1(i,k)=outcfu1(i,k)/ensdim
- outcfd1(i,k)=outcfd1(i,k)/ensdim
- outdfu1(i,k)=outdfu1(i,k)/ensdim
- outefu1(i,k)=outefu1(i,k)/ensdim
- outdfd1(i,k)=outdfd1(i,k)/ensdim
- outefd1(i,k)=outefd1(i,k)/ensdim
- end if !ierr
- end do !i
- end do !k
- end if !iens .eq. 1
- end if !l_flux
- END SUBROUTINE cup_output_ens
- SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
- kbcon,ktop,ierr, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! aa0 cloud work function
- ! gamma_cup = gamma on model cloud levels
- ! t_cup = temperature (Kelvin) on model cloud levels
- ! dby = buoancy term
- ! zu= normalized updraft mass flux
- ! z = heights of model levels
- ! ierr error value, maybe modified in this routine
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z,zu,gamma_cup,t_cup,dby
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon,ktop
- !
- ! input and output
- !
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite) &
- ,intent (out ) :: &
- aa0
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- real :: &
- dz,da
- !
- integer :: itf, ktf
- itf = MIN(ite,ide-1)
- ktf = MIN(kte,kde-1)
- do i=its,itf
- aa0(i)=0.
- enddo
- DO 100 k=kts+1,ktf
- DO 100 i=its,itf
- IF(ierr(i).ne.0)GO TO 100
- IF(K.LE.KBCON(I))GO TO 100
- IF(K.Gt.KTOP(I))GO TO 100
- DZ=Z(I,K)-Z(I,K-1)
- da=zu(i,k)*DZ*(9.81/(1004.*( &
- (T_cup(I,K)))))*DBY(I,K-1)/ &
- (1.+GAMMA_CUP(I,K))
- IF(K.eq.KTOP(I).and.da.le.0.)go to 100
- AA0(I)=AA0(I)+da
- if(aa0(i).lt.0.)aa0(i)=0.
- 100 continue
- END SUBROUTINE cup_up_aa0
- SUBROUTINE cup_up_he(k22,hkb,z_cup,cd,entr,he_cup,hc, &
- kbcon,ierr,dby,he,hes_cup, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! hc = cloud moist static energy
- ! hkb = moist static energy at originating level
- ! he = moist static energy on model levels
- ! he_cup = moist static energy on model cloud levels
- ! hes_cup = saturation moist static energy on model cloud levels
- ! dby = buoancy term
- ! cd= detrainment function
- ! z_cup = heights of model cloud levels
- ! entr = entrainment rate
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- he,he_cup,hes_cup,z_cup,cd
- ! entr= entrainment rate
- real &
- ,intent (in ) :: &
- entr
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon,k22
- !
- ! input and output
- !
- ! ierr error value, maybe modified in this routine
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- hc,dby
- real, dimension (its:ite) &
- ,intent (out ) :: &
- hkb
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- real :: &
- dz
- !
- integer :: itf, ktf
- itf = MIN(ite,ide-1)
- ktf = MIN(kte,kde-1)
- !
- !--- moist static energy inside cloud
- !
- do i=its,itf
- if(ierr(i).eq.0.)then
- hkb(i)=he_cup(i,k22(i))
- do k=1,k22(i)
- hc(i,k)=he_cup(i,k)
- DBY(I,K)=0.
- enddo
- do k=k22(i),kbcon(i)-1
- hc(i,k)=hkb(i)
- DBY(I,K)=0.
- enddo
- k=kbcon(i)
- hc(i,k)=hkb(i)
- DBY(I,Kbcon(i))=Hkb(I)-HES_cup(I,K)
- endif
- enddo
- do k=kts+1,ktf
- do i=its,itf
- if(k.gt.kbcon(i).and.ierr(i).eq.0.)then
- DZ=Z_cup(i,K)-Z_cup(i,K-1)
- HC(i,K)=(HC(i,K-1)*(1.-.5*CD(i,K)*DZ)+entr* &
- DZ*HE(i,K-1))/(1.+entr*DZ-.5*cd(i,k)*dz)
- DBY(I,K)=HC(I,K)-HES_cup(I,K)
- endif
- enddo
- enddo
- END SUBROUTINE cup_up_he
- SUBROUTINE cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
- kbcon,ktop,cd,dby,mentr_rate,clw_all, &
- q,GAMMA_cup,zu,qes_cup,k22,qe_cup,xl, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! cd= detrainment function
- ! q = environmental q on model levels
- ! qe_cup = environmental q on model cloud levels
- ! qes_cup = saturation q on model cloud levels
- ! dby = buoancy term
- ! cd= detrainment function
- ! zu = normalized updraft mass flux
- ! gamma_cup = gamma on model cloud levels
- ! mentr_rate = entrainment rate
- !
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- q,zu,gamma_cup,qe_cup,dby,qes_cup,z_cup,cd
- ! entr= entrainment rate
- real &
- ,intent (in ) :: &
- mentr_rate,xl
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon,ktop,k22
- !
- ! input and output
- !
- ! ierr error value, maybe modified in this routine
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- ! qc = cloud q (including liquid water) after entrainment
- ! qrch = saturation q in cloud
- ! qrc = liquid water content in cloud after rainout
- ! pw = condensate that will fall out at that level
- ! pwav = totan normalized integrated condensate (I1)
- ! c0 = conversion rate (cloud to rain)
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- qc,qrc,pw,clw_all
- real, dimension (its:ite) &
- ,intent (out ) :: &
- pwav
- !
- ! local variables in this routine
- !
- integer :: &
- iall,i,k
- real :: &
- dh,qrch,c0,dz,radius
- !
- integer :: itf, ktf
- itf = MIN(ite,ide-1)
- ktf = MIN(kte,kde-1)
- iall=0
- c0=.002
- !
- !--- no precip for small clouds
- !
- if(mentr_rate.gt.0.)then
- radius=.2/mentr_rate
- if(radius.lt.900.)c0=0.
- ! if(radius.lt.900.)iall=0
- endif
- do i=its,itf
- pwav(i)=0.
- enddo
- do k=kts,ktf
- do i=its,itf
- pw(i,k)=0.
- if(ierr(i).eq.0)qc(i,k)=qes_cup(i,k)
- clw_all(i,k)=0.
- qrc(i,k)=0.
- enddo
- enddo
- do i=its,itf
- if(ierr(i).eq.0.)then
- do k=k22(i),kbcon(i)-1
- qc(i,k)=qe_cup(i,k22(i))
- enddo
- endif
- enddo
- DO 100 k=kts+1,ktf
- DO 100 i=its,itf
- IF(ierr(i).ne.0)GO TO 100
- IF(K.Lt.KBCON(I))GO TO 100
- IF(K.Gt.KTOP(I))GO TO 100
- DZ=Z_cup(i,K)-Z_cup(i,K-1)
- !
- !------ 1. steady state plume equation, for what could
- !------ be in cloud without condensation
- !
- !
- QC(i,K)=(QC(i,K-1)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
- DZ*Q(i,K-1))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
- !
- !--- saturation in cloud, this is what is allowed to be in it
- !
- QRCH=QES_cup(I,K)+(1./XL)*(GAMMA_cup(i,k) &
- /(1.+GAMMA_cup(i,k)))*DBY(I,K)
- !
- !------- LIQUID WATER CONTENT IN cloud after rainout
- !
- clw_all(i,k)=QC(I,K)-QRCH
- QRC(I,K)=(QC(I,K)-QRCH)/(1.+C0*DZ*zu(i,k))
- if(qrc(i,k).lt.0.)then
- qrc(i,k)=0.
- endif
- !
- !------- 3.Condensation
- !
- PW(i,k)=c0*dz*QRC(I,K)*zu(i,k)
- if(iall.eq.1)then
- qrc(i,k)=0.
- pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
- if(pw(i,k).lt.0.)pw(i,k)=0.
- endif
- !
- !----- set next level
- !
- QC(I,K)=QRC(I,K)+qrch
- !
- !--- integrated normalized ondensate
- !
- PWAV(I)=PWAV(I)+PW(I,K)
- 100 CONTINUE
- END SUBROUTINE cup_up_moisture
- SUBROUTINE cup_up_nms(zu,z_cup,entr,cd,kbcon,ktop,ierr,k22, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- IMPLICIT NONE
- !
- ! on input
- !
- ! only local wrf dimensions are need as of now in this routine
- integer &
- ,intent (in ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- ! cd= detrainment function
- real, dimension (its:ite,kts:kte) &
- ,intent (in ) :: &
- z_cup,cd
- ! entr= entrainment rate
- real &
- ,intent (in ) :: &
- entr
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kbcon,ktop,k22
- !
- ! input and output
- !
- ! ierr error value, maybe modified in this routine
- integer, dimension (its:ite) &
- ,intent (inout) :: &
- ierr
- ! zu is the normalized mass flux
- real, dimension (its:ite,kts:kte) &
- ,intent (out ) :: &
- zu
- !
- ! local variables in this routine
- !
- integer :: &
- i,k
- real :: &
- dz
- integer :: itf, ktf
- itf = MIN(ite,ide-1)
- ktf = MIN(kte,kde-1)
- !
- ! initialize for this go around
- !
- do k=kts,ktf
- do i=its,itf
- zu(i,k)=0.
- enddo
- enddo
- !
- ! do normalized mass budget
- !
- do i=its,itf
- IF(ierr(I).eq.0)then
- do k=k22(i),kbcon(i)
- zu(i,k)=1.
- enddo
- DO K=KBcon(i)+1,KTOP(i)
- DZ=Z_cup(i,K)-Z_cup(i,K-1)
- ZU(i,K)=ZU(i,K-1)*(1.+(entr-cd(i,k))*DZ)
- enddo
- endif
- enddo
- END SUBROUTINE cup_up_nms
- !====================================================================
- SUBROUTINE gdinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
- MASS_FLUX,cp,restart, &
- P_QC,P_QI,P_FIRST_SCALAR, &
- RTHFTEN, RQVFTEN, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI, &
- allowed_to_read, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- !--------------------------------------------------------------------
- IMPLICIT NONE
- !--------------------------------------------------------------------
- LOGICAL , INTENT(IN) :: restart,allowed_to_read
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
- REAL, INTENT(IN) :: cp
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQCCUTEN, &
- RQICUTEN
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
- RTHFTEN, &
- RQVFTEN
- REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI, &
- MASS_FLUX
- INTEGER :: i, j, k, itf, jtf, ktf
-
- jtf=min0(jte,jde-1)
- ktf=min0(kte,kde-1)
- itf=min0(ite,ide-1)
-
- IF(.not.restart)THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RTHCUTEN(i,k,j)=0.
- RQVCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RTHFTEN(i,k,j)=0.
- RQVFTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- IF (P_QC .ge. P_FIRST_SCALAR) THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RQCCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- IF (P_QI .ge. P_FIRST_SCALAR) THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RQICUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- DO j=jts,jtf
- DO i=its,itf
- mass_flux(i,j)=0.
- ENDDO
- ENDDO
- ENDIF
- DO j=jts,jtf
- DO i=its,itf
- APR_GR(i,j)=0.
- APR_ST(i,j)=0.
- APR_W(i,j)=0.
- APR_MC(i,j)=0.
- APR_AS(i,j)=0.
- APR_CAPMA(i,j)=0.
- APR_CAPME(i,j)=0.
- APR_CAPMI(i,j)=0.
- ENDDO
- ENDDO
- END SUBROUTINE gdinit
- SUBROUTINE massflx_stats(xf_ens,ensdim,maxens,maxens2,maxens3, &
- xt_ave,xt_std,xt_cur,xt_ske,j,ierr,itest, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI, &
- pr_gr,pr_w,pr_mc,pr_st,pr_as, &
- pr_capma,pr_capme,pr_capmi, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- IMPLICIT NONE
- integer, intent (in ) :: &
- j,ensdim,maxens3,maxens,maxens2,itest
- INTEGER, INTENT(IN ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- real, dimension (its:ite) &
- , intent(inout) :: &
- xt_ave,xt_cur,xt_std,xt_ske
- integer, dimension (its:ite), intent (in) :: &
- ierr
- real, dimension (ims:ime,jms:jme,1:ensdim) &
- , intent(in ) :: &
- xf_ens
- real, dimension (ims:ime,jms:jme) &
- , intent(inout) :: &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- APR_CAPMA,APR_CAPME,APR_CAPMI
- real, dimension (its:ite,jts:jte) &
- , intent(inout) :: &
- pr_gr,pr_w,pr_mc,pr_st,pr_as, &
- pr_capma,pr_capme,pr_capmi
- !
- ! local stuff
- !
- real, dimension (its:ite , 1:maxens3 ) :: &
- x_ave,x_cur,x_std,x_ske
- real, dimension (its:ite , 1:maxens ) :: &
- x_ave_cap
- integer, dimension (1:maxens3) :: nc1
- integer :: i,k
- integer :: num,kk,num2,iedt
- real :: a3,a4
- num=ensdim/maxens3
- num2=ensdim/maxens
- if(itest.eq.1)then
- do i=its,ite
- pr_gr(i,j) = 0.
- pr_w(i,j) = 0.
- pr_mc(i,j) = 0.
- pr_st(i,j) = 0.
- pr_as(i,j) = 0.
- pr_capma(i,j) = 0.
- pr_capme(i,j) = 0.
- pr_capmi(i,j) = 0.
- enddo
- endif
- do k=1,maxens
- do i=its,ite
- x_ave_cap(i,k)=0.
- enddo
- enddo
- do k=1,maxens3
- do i=its,ite
- x_ave(i,k)=0.
- x_std(i,k)=0.
- x_ske(i,k)=0.
- x_cur(i,k)=0.
- enddo
- enddo
- do i=its,ite
- xt_ave(i)=0.
- xt_std(i)=0.
- xt_ske(i)=0.
- xt_cur(i)=0.
- enddo
- do kk=1,num
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0)then
- x_ave(i,k)=x_ave(i,k)+xf_ens(i,j,maxens3*(kk-1)+k)
- endif
- enddo
- enddo
- enddo
- do iedt=1,maxens2
- do k=1,maxens
- do kk=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0)then
- x_ave_cap(i,k)=x_ave_cap(i,k) &
- +xf_ens(i,j,maxens3*(k-1)+(iedt-1)*maxens*maxens3+kk)
- endif
- enddo
- enddo
- enddo
- enddo
- do k=1,maxens
- do i=its,ite
- if(ierr(i).eq.0)then
- x_ave_cap(i,k)=x_ave_cap(i,k)/float(num2)
- endif
- enddo
- enddo
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0)then
- x_ave(i,k)=x_ave(i,k)/float(num)
- endif
- enddo
- enddo
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0)then
- xt_ave(i)=xt_ave(i)+x_ave(i,k)
- endif
- enddo
- enddo
- do i=its,ite
- if(ierr(i).eq.0)then
- xt_ave(i)=xt_ave(i)/float(maxens3)
- endif
- enddo
- !
- !--- now do std, skewness,curtosis
- !
- do kk=1,num
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0.and.x_ave(i,k).gt.0.)then
- ! print *,i,j,k,kk,x_std(i,k),xf_ens(i,j,maxens3*(kk-1)+k),x_ave(i,k)
- x_std(i,k)=x_std(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**2
- x_ske(i,k)=x_ske(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**3
- x_cur(i,k)=x_cur(i,k)+(xf_ens(i,j,maxens3*(kk-1)+k)-x_ave(i,k))**4
- endif
- enddo
- enddo
- enddo
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0.and.xt_ave(i).gt.0.)then
- xt_std(i)=xt_std(i)+(x_ave(i,k)-xt_ave(i))**2
- xt_ske(i)=xt_ske(i)+(x_ave(i,k)-xt_ave(i))**3
- xt_cur(i)=xt_cur(i)+(x_ave(i,k)-xt_ave(i))**4
- endif
- enddo
- enddo
- do k=1,maxens3
- do i=its,ite
- if(ierr(i).eq.0.and.x_std(i,k).gt.0.)then
- x_std(i,k)=x_std(i,k)/float(num)
- a3=max(1.e-6,x_std(i,k))
- x_std(i,k)=sqrt(a3)
- a3=max(1.e-6,x_std(i,k)**3)
- a4=max(1.e-6,x_std(i,k)**4)
- x_ske(i,k)=x_ske(i,k)/float(num)/a3
- x_cur(i,k)=x_cur(i,k)/float(num)/a4
- endif
- ! print*,' '
- ! print*,'Some statistics at gridpoint i,j, ierr',i,j,ierr(i)
- ! print*,'statistics for closure number ',k
- ! print*,'Average= ',x_ave(i,k),' Std= ',x_std(i,k)
- ! print*,'Skewness= ',x_ske(i,k),' Curtosis= ',x_cur(i,k)
- ! print*,' '
- enddo
- enddo
- do i=its,ite
- if(ierr(i).eq.0.and.xt_std(i).gt.0.)then
- xt_std(i)=xt_std(i)/float(maxens3)
- a3=max(1.e-6,xt_std(i))
- xt_std(i)=sqrt(a3)
- a3=max(1.e-6,xt_std(i)**3)
- a4=max(1.e-6,xt_std(i)**4)
- xt_ske(i)=xt_ske(i)/float(maxens3)/a3
- xt_cur(i)=xt_cur(i)/float(maxens3)/a4
- ! print*,' '
- ! print*,'Total ensemble independent statistics at i =',i
- ! print*,'Average= ',xt_ave(i),' Std= ',xt_std(i)
- ! print*,'Skewness= ',xt_ske(i),' Curtosis= ',xt_cur(i)
- ! print*,' '
- !
- ! first go around: store massflx for different closures/caps
- !
- if(itest.eq.1)then
- pr_gr(i,j) = .333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))
- pr_w(i,j) = .333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))
- pr_mc(i,j) = .333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))
- pr_st(i,j) = .333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))
- pr_as(i,j) = .25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
- + x_ave(i,16))
- pr_capma(i,j) = x_ave_cap(i,1)
- pr_capme(i,j) = x_ave_cap(i,2)
- pr_capmi(i,j) = x_ave_cap(i,3)
- !
- ! second go around: store preciprates (mm/hour) for different closures/caps
- !
- else if (itest.eq.2)then
- APR_GR(i,j)=.333*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3))* &
- 3600.*pr_gr(i,j) +APR_GR(i,j)
- APR_W(i,j)=.333*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6))* &
- 3600.*pr_w(i,j) +APR_W(i,j)
- APR_MC(i,j)=.333*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9))* &
- 3600.*pr_mc(i,j) +APR_MC(i,j)
- APR_ST(i,j)=.333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))* &
- 3600.*pr_st(i,j) +APR_ST(i,j)
- APR_AS(i,j)=.25*(x_ave(i,13)+x_ave(i,14)+x_ave(i,15) &
- + x_ave(i,16))* &
- 3600.*pr_as(i,j) +APR_AS(i,j)
- APR_CAPMA(i,j) = x_ave_cap(i,1)* &
- 3600.*pr_capma(i,j) +APR_CAPMA(i,j)
- APR_CAPME(i,j) = x_ave_cap(i,2)* &
- 3600.*pr_capme(i,j) +APR_CAPME(i,j)
- APR_CAPMI(i,j) = x_ave_cap(i,3)* &
- 3600.*pr_capmi(i,j) +APR_CAPMI(i,j)
- endif
- endif
- enddo
- END SUBROUTINE massflx_stats
- SUBROUTINE neg_check(dt,q,outq,outt,outqc,pret,its,ite,kts,kte,itf,ktf)
- INTEGER, INTENT(IN ) :: its,ite,kts,kte,itf,ktf
- real, dimension (its:ite,kts:kte ) , &
- intent(inout ) :: &
- q,outq,outt,outqc
- real, dimension (its:ite ) , &
- intent(inout ) :: &
- pret
- real &
- ,intent (in ) :: &
- dt
- real :: thresh,qmem,qmemf,qmem2,qtest,qmem1
- !
- ! first do check on vertical heating rate
- !
- thresh=200.01
- do i=its,itf
- qmemf=1.
- qmem=0.
- do k=kts,ktf
- qmem=outt(i,k)*86400.
- if(qmem.gt.2.*thresh)then
- qmem2=2.*thresh/qmem
- qmemf=min(qmemf,qmem2)
- !
- !
- ! print *,'1',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
- endif
- if(qmem.lt.-thresh)then
- qmem2=-thresh/qmem
- qmemf=min(qmemf,qmem2)
- !
- !
- ! print *,'2',' adjusted massflux by factor ',i,k,qmem,qmem2,qmemf
- endif
- enddo
- do k=kts,ktf
- outq(i,k)=outq(i,k)*qmemf
- outt(i,k)=outt(i,k)*qmemf
- outqc(i,k)=outqc(i,k)*qmemf
- enddo
- pret(i)=pret(i)*qmemf
- enddo
- !
- ! check whether routine produces negative q's. This can happen, since
- ! tendencies are calculated based on forced q's. This should have no
- ! influence on conservation properties, it scales linear through all
- ! tendencies
- !
- thresh=1.e-10
- do i=its,itf
- qmemf=1.
- do k=kts,ktf
- qmem=outq(i,k)
- if(abs(qmem).gt.0.)then
- qtest=q(i,k)+outq(i,k)*dt
- if(qtest.lt.thresh)then
- !
- ! qmem2 would be the maximum allowable tendency
- !
- qmem1=outq(i,k)
- qmem2=(thresh-q(i,k))/dt
- qmemf=min(qmemf,qmem2/qmem1)
- ! qmemf=max(0.,qmemf)
- ! print *,'4 adjusted tendencies ',i,k,qmem,qmem2,qmemf
- endif
- endif
- enddo
- do k=kts,ktf
- outq(i,k)=outq(i,k)*qmemf
- outt(i,k)=outt(i,k)*qmemf
- outqc(i,k)=outqc(i,k)*qmemf
- enddo
- pret(i)=pret(i)*qmemf
- enddo
- END SUBROUTINE neg_check
- !-------------------------------------------------------
- END MODULE module_cu_gd