/wrfv2_fire/phys/module_cu_g3.F
FORTRAN Legacy | 5313 lines | 3667 code | 294 blank | 1352 comment | 275 complexity | 39bd7e8616606edac6d7d990e6a7ef28 MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !WRF:MODEL_LAYER:PHYSICS
- !
- MODULE module_cu_g3
- CONTAINS
- !-------------------------------------------------------------
- SUBROUTINE G3DRV( &
- 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,edt_out &
- ,GDC,GDC2 ,kpbl,k22_shallow,kbcon_shallow &
- ,ktop_shallow,xmb_shallow &
- ,cugd_tten,cugd_qvten ,cugd_qcten &
- ,cugd_ttens,cugd_qvtens,cugd_avedx,imomentum &
- ,ensdim,maxiens,maxens,maxens2,maxens3,ichoice &
- ,ishallow_g3,ids,ide, jds,jde, kds,kde &
- ,ims,ime, jms,jme, kms,kme &
- ,ips,ipe, jps,jpe, kps,kpe &
- ,its,ite, jts,jte, kts,kte &
- ,periodic_x,periodic_y &
- ,RQVCUTEN,RQCCUTEN,RQICUTEN &
- ,RQVFTEN,RTHFTEN,RTHCUTEN &
- ,rqvblten,rthblten &
- ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
- #if ( WRF_DFI_RADAR == 1 )
- ! Optional CAP suppress option
- ,do_capsuppress,cap_suppress_loc &
- #endif
- )
- !-------------------------------------------------------------
- IMPLICIT NONE
- !-------------------------------------------------------------
- INTEGER, INTENT(IN ) :: &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- ips,ipe, jps,jpe, kps,kpe, &
- its,ite, jts,jte, kts,kte
- LOGICAL periodic_x,periodic_y
- integer, parameter :: ens4_spread = 3 ! max(3,cugd_avedx)
- integer, parameter :: ens4=ens4_spread*ens4_spread
- integer, intent (in ) :: &
- ensdim,maxiens,maxens,maxens2,maxens3,ichoice
-
- INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP,cugd_avedx, &
- ishallow_g3,imomentum
- 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
- INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: KPBL
- INTEGER, DIMENSION( ims:ime , jms:jme ),INTENT(INOUT) :: k22_shallow, &
- kbcon_shallow,ktop_shallow
- !
- REAL, INTENT(IN ) :: DT, DX
- !
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: pratec,RAINCV, MASS_FLUX, &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
- edt_out,APR_CAPMA,APR_CAPME,APR_CAPMI, &
- htop,hbot,xmb_shallow
- !+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, &
- cugd_tten,cugd_qvten,cugd_qcten, &
- cugd_ttens,cugd_qvtens, &
- RQVFTEN
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQVBLTEN, &
- RTHBLTEN, &
- RQCCUTEN, &
- RQICUTEN
- !
- ! 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
- #if ( WRF_DFI_RADAR == 1 )
- !
- ! option of cap suppress:
- ! do_capsuppress = 1 do
- ! do_capsuppress = other don't
- !
- !
- INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
- REAL, DIMENSION( ims:ime, jms:jme ),INTENT(IN ),OPTIONAL :: cap_suppress_loc
- REAL, DIMENSION( its:ite ) :: cap_suppress_j
- #endif
- ! LOCAL VARS
- real, dimension(ims:ime,jms:jme,1:ensdim),intent(inout) :: &
- xf_ens,pr_ens
- real, dimension ( its:ite , jts:jte , 1:ensdim) :: &
- massflni,xfi_ens,pri_ens
- REAL, DIMENSION( its:ite , jts:jte ) :: MASSI_FLX, &
- APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
- edti_out,APRi_CAPMA,APRi_CAPME,APRi_CAPMI,gswi
- real, dimension (its:ite,kts:kte) :: &
- SUBT,SUBQ,OUTT,OUTQ,OUTQC,phh,subm,cupclw,dhdt, &
- outts,outqs
- real, dimension (its:ite) :: &
- pret, ter11, aa0, fp,xlandi
- !+lxz
- integer, dimension (its:ite) :: &
- kbcon, ktop,kpbli,k22s,kbcons,ktops
- !.lxz
- integer, dimension (its:ite,jts:jte) :: &
- iact_old_gr
- integer :: iens,ibeg,iend,jbeg,jend,n,nn,ens4n
- integer :: ibegh,iendh,jbegh,jendh
- integer :: ibegc,iendc,jbegc,jendc
- !
- ! 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) :: &
- T2d,q2d,PO,P2d,US,VS,tn,qo,tshall,qshall
- real, dimension (ips-2:ipe+2,kps:kpe,jps-2:jpe+2) :: &
- ave_f_t,ave_f_q
- real, dimension (its:ite,kts:kte,1:ens4) :: &
- omeg,tx,qx
- real, dimension (its:ite) :: &
- Z1,PSUR,AAEQ,direction,cuten,umean,vmean,pmean,xmbs
- real, dimension (its:ite,1:ens4) :: &
- mconv
- INTEGER :: i,j,k,ICLDCK,ipr,jpr
- REAL :: tcrit,tscl_KF,dp,dq,sub_spread,subcenter
- INTEGER :: itf,jtf,ktf,iss,jss,nbegin,nend
- INTEGER :: high_resolution
- REAL :: rkbcon,rktop !-lxz
- ! ruc variable
- real, dimension (its:ite) :: tkm
- ! A. Betts for shallow convection: suggestion for the KF timescale < DELTAX / 25 m/s
- tscl_kf=dx/25.
- !
- ! write(0,*)'ishallow = ',ishallow_g3
- high_resolution=0
- if(cugd_avedx.gt.1) high_resolution=1
- subcenter=0.
- ! subcenter=1./float(cugd_avedx)
- sub_spread=max(1.,float(cugd_avedx*cugd_avedx-1))
- sub_spread=(1.-subcenter)/sub_spread
- iens=1
- ipr=43
- jpr=1
- ipr=0
- jpr=0
- ! if(itimestep.eq.8)then
- ! ipr=37
- ! jpr=16
- ! endif
- IF ( periodic_x ) THEN
- ibeg=max(its,ids)
- iend=min(ite,ide-1)
- ibegc=max(its,ids)
- iendc=min(ite,ide-1)
- ELSE
- ibeg=max(its,ids)
- iend=min(ite,ide-1)
- ibegc=max(its,ids+4)
- iendc=min(ite,ide-5)
- END IF
- IF ( periodic_y ) THEN
- jbeg=max(jts,jds)
- jend=min(jte,jde-1)
- jbegc=max(jts,jds)
- jendc=min(jte,jde-1)
- ELSE
- jbeg=max(jts,jds)
- jend=min(jte,jde-1)
- jbegc=max(jts,jds+4)
- jendc=min(jte,jde-5)
- END IF
- do j=jts,jte
- do i=its,ite
- k22_shallow(i,j)=0
- kbcon_shallow(i,j)=0
- ktop_shallow(i,j)=0
- xmb_shallow(i,j)=0
- enddo
- enddo
- tcrit=258.
- ave_f_t=0.
- ave_f_q=0.
- itf=MIN(ite,ide-1)
- ktf=MIN(kte,kde-1)
- jtf=MIN(jte,jde-1)
- !
- #if ( EM_CORE == 1 )
- if(high_resolution.eq.1)then
- !
- ! calculate these on the halo...the incominh tendencies have been exchanged on a 24pt halo
- ! only neede for high resolution run
- !
- ibegh=its
- jbegh=jts
- iendh=ite
- jendh=jte
- if(its.eq.ips)ibegh=max(its-1,ids)
- if(jts.eq.jps)jbegh=max(jts-1,jds)
- if(jte.eq.jpe)jendh=min(jte+1,jde-1)
- if(ite.eq.ipe)iendh=min(ite+1,ide-1)
- DO J = jbegh,jendh
- DO k= kts,ktf
- DO I= ibegh,iendh
- ave_f_t(i,k,j)=(rthften(i-1,k,j-1)+rthften(i-1,k,j) + rthften(i-1,k,j+1)+ &
- rthften(i,k,j-1) +rthften(i,k,j) +rthften(i,k,j+1)+ &
- rthften(i+1,k,j-1) +rthften(i+1,k,j) +rthften(i+1,k,j+1))/9.
- ave_f_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) + rqvften(i-1,k,j+1)+ &
- rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
- rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
- ! ave_f_t(i,k,j)=rthften(i,k,j)
- ! ave_f_q(i,k,j)=rqvften(i,k,j)
- ENDDO
- ENDDO
- ENDDO
- endif
- #endif
- DO 100 J = jts,jtf
- DO n= 1,ensdim
- DO I= its,itf
- xfi_ens(i,j,n)=0.
- pri_ens(i,j,n)=0.
- ! xfi_ens(i,j,n)=xf_ens(i,j,n)
- ! pri_ens(i,j,n)=pr_ens(i,j,n)
- ENDDO
- ENDDO
- DO I= its,itf
- kbcon(i)=0
- ktop(i)=0
- tkm(i)=0.
- HBOT(I,J) =REAL(KTE)
- HTOP(I,J) =REAL(KTS)
- iact_old_gr(i,j)=0
- mass_flux(i,j)=0.
- massi_flx(i,j)=0.
- raincv(i,j)=0.
- pratec (i,j)=0.
- edt_out(i,j)=0.
- edti_out(i,j)=0.
- gswi(i,j)=gsw(i,j)
- xlandi(i)=xland(i,j)
- APRi_GR(i,j)=apr_gr(i,j)
- APRi_w(i,j)=apr_w(i,j)
- APRi_mc(i,j)=apr_mc(i,j)
- APRi_st(i,j)=apr_st(i,j)
- APRi_as(i,j)=apr_as(i,j)
- APRi_capma(i,j)=apr_capma(i,j)
- APRi_capme(i,j)=apr_capme(i,j)
- APRi_capmi(i,j)=apr_capmi(i,j)
- CU_ACT_FLAG(i,j) = .true.
- ENDDO
- do k=kts,kte
- DO I= its,itf
- cugd_tten(i,k,j)=0.
- cugd_ttens(i,k,j)=0.
- cugd_qvten(i,k,j)=0.
- cugd_qvtens(i,k,j)=0.
- cugd_qcten(i,k,j)=0.
- ENDDO
- ENDDO
- DO n=1,ens4
- DO I= its,itf
- mconv(i,n)=0.
- ENDDO
- do k=kts,kte
- DO I= its,itf
- omeg(i,k,n)=0.
- tx(i,k,n)=0.
- qx(i,k,n)=0.
- ENDDO
- ENDDO
- ENDDO
- DO k=1,ensdim
- DO I= its,itf
- massflni(i,j,k)=0.
- ENDDO
- ENDDO
- ! 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
- ! PSUR(I)=p(I,1,J)*.01
- TER11(I)=HT(i,j)
- aaeq(i)=0.
- direction(i)=0.
- pret(i)=0.
- umean(i)=0.
- vmean(i)=0.
- pmean(i)=0.
- kpbli(i)=kpbl(i,j)
- ENDDO
- ! if(j.eq.jpr)write(0,*)'psur(ipr),ter11(ipr),kpbli(ipr)'
- ! if(j.eq.jpr)write(0,*)psur(ipr),ter11(ipr),kpbli(ipr),r_v
- DO K=kts,ktf
- DO I=ITS,ITF
- po(i,k)=phh(i,k)*.01
- subm(i,k)=0.
- 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)
- IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
- SUBT(I,K)=0.
- SUBQ(I,K)=0.
- OUTT(I,K)=0.
- OUTQ(I,K)=0.
- OUTQC(I,K)=0.
- OUTTS(I,K)=0.
- OUTQS(I,K)=0.
- TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
- QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
- TSHALL(I,K)=t2d(i,k)+RTHBLTEN(i,k,j)*pi(i,k,j)*dt
- DHDT(I,K)=cp*RTHBLTEN(i,k,j)*pi(i,k,j)+ XLV*RQVBLTEN(i,k,j)
- QSHALL(I,K)=q2d(i,k)+RQVBLTEN(i,k,j)*dt
- if(high_resolution.eq.1)then
- TN(I,K)=t2d(i,k)+ave_f_t(i,k,j)*dt
- QO(I,K)=q2d(i,k)+ave_f_q(i,k,j)*dt
- endif
- IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
- IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! write(0,123)k,p2d(i,k),t2d(i,k),tn(i,k),q2d(i,k),QO(i,k),RTHBLTEN(i,k,j),RQVBLTEN(i,k,j)
- ! endif
- ENDDO
- ENDDO
- 123 format(1x,i2,f8.0,1x,2(1x,f8.3),4(1x,e12.4))
- ens4n=0
- nbegin=0
- nend=0
- if(ens4_spread.gt.1)then
- nbegin=-ens4_spread/2
- nend=ens4_spread/2
- endif
- do nn=nbegin,nend,1
- jss=max(j+nn,jds+0)
- jss=min(jss,jde-1)
- do n=nbegin,nend,1
- ens4n=ens4n+1
- DO K=kts,ktf
- DO I=ITS,ITF
- iss=max(i+n,ids+0)
- iss=min(iss,ide-1)
- omeg(I,K,ens4n)= -g*rho(i,k,j)*w(iss,k,jss)
- ! omeg(I,K,ens4n)= -g*rho(i,k,j)*w(i,k,j)
- Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(iss,k,jss)*dt
- ! Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(i,k,j)*dt
- if(high_resolution.eq.1)Tx(I,K,ens4n)=t2d(i,k)+ave_f_t(iss,k,jss)*dt
- IF(Tx(I,K,ens4n).LT.200.)Tx(I,K,ens4n)=T2d(I,K)
- Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(iss,k,jss)*dt
- Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(i,k,j)*dt
- if(high_resolution.eq.1)qx(I,K,ens4n)=q2d(i,k)+ave_f_q(iss,k,jss)*dt
- IF(Qx(I,K,ens4n).LT.1.E-08)Qx(I,K,ens4n)=1.E-08
- enddo
- enddo
- enddo !n
- enddo !nn
- 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
- 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.
- ENDDO
- do n=1,ens4
- DO K=kts,ktf-1
- DO I = its,itf
- dq=(q2d(i,k+1)-q2d(i,k))
- mconv(i,n)=mconv(i,n)+omeg(i,k,n)*dq/g
- enddo
- ENDDO
- ENDDO
- do n=1,ens4
- DO I = its,itf
- if(mconv(i,n).lt.0.)mconv(i,n)=0.
- ENDDO
- ENDDO
- !
- !---- CALL CUMULUS PARAMETERIZATION
- !
- #if ( WRF_DFI_RADAR == 1 )
- if(do_capsuppress == 1 ) then
- DO I= its,itf
- cap_suppress_j(i)=cap_suppress_loc(i,j)
- ENDDO
- endif
- #endif
- CALL CUP_enss_3d(outqc,j,AAEQ,T2d,Q2d,TER11,subm,TN,QO,PO,PRET,&
- P2d,OUTT,OUTQ,DT,itimestep,tkm,PSUR,US,VS,tcrit,iens,tx,qx, &
- tshall,qshall,kpbli,DHDT,outts,outqs,tscl_kf, &
- k22s,kbcons,ktops,xmbs, &
- mconv,massflni,iact_old_gr,omeg,direction,MASSi_FLX, &
- maxiens,maxens,maxens2,maxens3,ensdim, &
- APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
- APRi_CAPMA,APRi_CAPME,APRi_CAPMI,kbcon,ktop,cupclw, &
- xfi_ens,pri_ens,XLANDi,gswi,edti_out,subt,subq, &
- ! ruc lv_p,rv_p,cpd_p,g0_p,ichoice,ipr,jpr, &
- xlv,r_v,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
- ishallow_g3,itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte &
- #if ( WRF_DFI_RADAR == 1 )
- ,do_capsuppress,cap_suppress_j &
- #endif
- )
- if(j.lt.jbegc.or.j.gt.jendc)go to 100
- DO I=ibegc,iendc
- xmb_shallow(i,j)=xmbs(i)
- k22_shallow(i,j)=k22s(i)
- kbcon_shallow(i,j)=kbcons(i)
- ktop_shallow(i,j)=ktops(i)
- cuten(i)=0.
- if(pret(i).gt.0.)then
- cuten(i)=1.
- ! raincv(i,j)=pret(i)*dt
- endif
- ENDDO
- ! if(j.eq.jpr)write(0,*)'precip,ktop,kbcon = ',pret(ipr),ktop(ipr),kbcon(ipr)
- DO I=ibegc,iendc
- DO K=kts,ktf
- cugd_ttens(I,K,J)=subt(i,k)*cuten(i)*sub_spread
- cugd_qvtens(I,K,J)=subq(i,k)*cuten(i)*sub_spread
- ! cugd_tten(I,K,J)=outt(i,k)*cuten(i)
- ! cugd_qvten(I,K,J)=outq(i,k)*cuten(i)
- cugd_tten(I,K,J)=outts(i,k)+outt(i,k)*cuten(i)
- cugd_qvten(I,K,J)=outqs(i,k)+outq(i,k)*cuten(i)
- cugd_qcten(I,K,J)=outqc(i,k)*cuten(i)
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! write(0,*)subt(i,k)+outt(i,k),subq(i,k)+outq(i,k),outts(i,k),outqs(i,k)
- ! endif
- ENDDO
- ENDDO
- DO I=ibegc,iendc
- if(pret(i).gt.0.)then
- raincv(i,j)=pret(i)*dt
- pratec(i,j)=pret(i)
- 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
- ENDDO
- DO n= 1,ensdim
- DO I= ibegc,iendc
- xf_ens(i,j,n)=xfi_ens(i,j,n)
- pr_ens(i,j,n)=pri_ens(i,j,n)
- ENDDO
- ENDDO
- DO I= ibegc,iendc
- APR_GR(i,j)=apri_gr(i,j)
- APR_w(i,j)=apri_w(i,j)
- APR_mc(i,j)=apri_mc(i,j)
- APR_st(i,j)=apri_st(i,j)
- APR_as(i,j)=apri_as(i,j)
- APR_capma(i,j)=apri_capma(i,j)
- APR_capme(i,j)=apri_capme(i,j)
- APR_capmi(i,j)=apri_capmi(i,j)
- mass_flux(i,j)=massi_flx(i,j)
- edt_out(i,j)=edti_out(i,j)
- ENDDO
- IF(PRESENT(RQCCUTEN)) THEN
- IF ( F_QC ) THEN
- DO K=kts,ktf
- DO I=ibegc,iendc
- 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=ibegc,iendc
- if(t2d(i,k).lt.258.)then
- RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
- cugd_qcten(i,k,j)=0.
- 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
- 100 continue
- END SUBROUTINE G3DRV
- SUBROUTINE CUP_enss_3d(OUTQC,J,AAEQ,T,Q,Z1,sub_mas, &
- TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,ktau,tkmax,PSUR,US,VS, &
- TCRIT,iens,tx,qx, &
- tshall,qshall,kpbl,dhdt,outts,outqs,tscl_kf, &
- k23,kbcon3,ktop3,xmb3, &
- 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,cupclw, & !-lxz
- xf_ens,pr_ens,xland,gsw,edt_out,subt,subq, &
- xl,rv,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
- ishallow_g3,itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte &
- #if ( WRF_DFI_RADAR == 1 )
- ! Optional CAP suppress option
- ,do_capsuppress,cap_suppress_j &
- #endif
- )
- IMPLICIT NONE
- integer &
- ,intent (in ) :: &
- itf,jtf,ktf,ktau, &
- its,ite, jts,jte, kts,kte,ipr,jpr,ens4,high_resolution
- integer, intent (in ) :: &
- j,ensdim,maxiens,ishallow_g3,maxens,maxens2,maxens3,ichoice,iens
- !
- !
- !
- real, dimension (its:ite,jts:jte,1:ensdim) &
- ,intent (inout) :: &
- massfln,xf_ens,pr_ens
- real, dimension (its:ite,jts:jte) &
- ,intent (inout ) :: &
- APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
- APR_CAPME,APR_CAPMI,massflx,edt_out
- real, dimension (its:ite,jts:jte) &
- ,intent (in ) :: &
- 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 (inout ) :: &
- DHDT,OUTT,OUTQ,OUTQC,subt,subq,sub_mas,cupclw,outts,outqs
- real, dimension (its:ite) &
- ,intent (out ) :: &
- pre,xmb3
- integer, dimension (its:ite) &
- ,intent (out ) :: &
- kbcon,ktop,k23,kbcon3,ktop3
- integer, dimension (its:ite) &
- ,intent (in ) :: &
- kpbl
- !
- ! 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,PO,P,US,VS,tn,tshall,qshall
- real, dimension (its:ite,kts:kte,1:ens4) &
- ,intent (inout ) :: &
- omeg,tx,qx
- real, dimension (its:ite,kts:kte) &
- ,intent (inout) :: &
- Q,QO
- real, dimension (its:ite) &
- ,intent (in ) :: &
- Z1,PSUR,AAEQ,direction,tkmax,xland
- real, dimension (its:ite,1:ens4) &
- ,intent (in ) :: &
- mconv
-
- real &
- ,intent (in ) :: &
- dtime,tcrit,xl,cp,rv,g,tscl_kf
- #if ( WRF_DFI_RADAR == 1 )
- !
- ! option of cap suppress:
- ! do_capsuppress = 1 do
- ! do_capsuppress = other don't
- !
- !
- INTEGER, INTENT(IN ) ,OPTIONAL :: do_capsuppress
- REAL, DIMENSION( its:ite ),INTENT(IN ) ,OPTIONAL :: cap_suppress_j
- #endif
- !
- ! 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,subt_ens,subq_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) :: &
- he3,hes3,qes3,z3,zdo3,zu3_0,hc3_0,dby3_0,gamma3_0_cup, &
- qes3_cup,q3_cup,he3_cup,hes3_cup,z3_cup,gamma3_cup,t3_cup, &
- xhe3,xhes3,xqes3,xz3,xt3,xq3, &
- xqes3_cup,xq3_cup,xhe3_cup,xhes3_cup,xz3_cup,xgamma3_cup, &
- xt3_cup, &
- xdby3,xqc3,xhc3,xqrc3,xzu3, &
- dby3,qc3,pw3,hc3,qrc3,zu3,cd3,DELLAH3,DELLAQ3, &
- dsubt3,dsubq3,DELLAT3,DELLAQC3
- 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,dsubt,dsubq
- ! 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) :: &
- aa3_0,aa3,hkb3,qkb3,pwav3,bu3,xaa3,xhkb3, &
- hkb3_0,edt,edto,edtx,AA1,AA0,XAA0,HKB, &
- HKBO,aad,XHKB,QKB,QKBO,edt3, &
- XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO, &
- PWEVO,BU,BUO,cap_max,xland1, &
- cap_max_increment,closure_n,cap_max3
- real, dimension (its:ite,1:ens4) :: &
- axx
- integer, dimension (its:ite) :: &
- kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x,jmin3,kdet3, & !-lxz
- KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX,ierr5,ierr5_0
- integer :: &
- nall,iedt,nens,nens3,ki,I,K,KK,iresult
- real :: &
- day,dz,mbdt,mbdt_s,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
- zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
- massfld,dh,cap_maxs,trash,entr_rate3,mentr_rate3
- integer :: jmini
- logical :: keep_going
- real xff_shal(9),blqe,xkshal
- day=86400.
- do i=its,itf
- xmb3(i)=0.
- closure_n(i)=16.
- xland1(i)=1.
- if(xland(i).gt.1.5)xland1(i)=0.
- ! cap_max_increment(i)=50.
- 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
- entr_rate3=.2/200.
- !
- !--- entrainment of mass
- !
- mentrd_rate=0.
- mentr_rate=entr_rate
- mentr_rate3=entr_rate3
- !
- !--- initial detrainmentrates
- !
- do k=kts,ktf
- do i=its,itf
- cupclw(i,k)=0.
- cd(i,k)=0.01*entr_rate
- cd3(i,k)=entr_rate3
- cdd(i,k)=0.
- zdo3(i,k)=0.
- hcdo(i,k)=0.
- qrcdo(i,k)=0.
- dellaqc(i,k)=0.
- enddo
- enddo
- !
- !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
- ! base mass flux
- !
- edtmax=1.
- edtmin=.2
- !
- !--- minimum depth (m), clouds must have
- !
- depth_min=500.
- !
- !--- maximum depth (mb) of capping
- !--- inversion (larger cap = no convection)
- !
- ! cap_maxs=125.
- cap_maxs=75.
- DO i=its,itf
- kbmax(i)=1
- jmin3(i)=0
- kdet3(i)=0
- aa0(i)=0.
- aa3_0(i)=0.
- aa1(i)=0.
- aa3(i)=0.
- aad(i)=0.
- edt(i)=0.
- edt3(i)=0.
- kstabm(i)=ktf-1
- IERR(i)=0
- IERR2(i)=0
- IERR3(i)=0
- IERR5(i)=0
- IERR5_0(i)=0
- enddo
- !
- !--- first check for upstream convection
- !
- #if ( WRF_DFI_RADAR == 1 )
- if(do_capsuppress == 1) then
- do i=its,itf
- cap_max(i)=cap_maxs
- cap_max3(i)=25.
- if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
- if (abs(cap_suppress_j(i) - 1.0 ) < 0.1 ) then
- cap_max(i)=cap_maxs+75.
- elseif (abs(cap_suppress_j(i) - 0.0 ) < 0.1 ) then
- cap_max(i)=10.0
- endif
- iresult=0
- enddo
- else
- do i=its,itf
- cap_max(i)=cap_maxs
- cap_max3(i)=25.
- if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
- iresult=0
- enddo
- endif
- do i=its,itf
- edt_out(i,j)=cap_max(i)
- enddo
- #else
- do i=its,itf
- cap_max(i)=cap_maxs
- cap_max3(i)=25.
- if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
- iresult=0
- enddo
- #endif
- !
- !--- 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
- !
- !--- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
- psur,ierr,tcrit,0,xl,cp, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- do i=its,itf
- if(aaeq(i).lt.-0.1)then
- ierr(i)=20
- endif
- ! if(ierr(i).eq.0)then
- !
- do k=kts,ktf
- 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, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !--- increase detrainment in stable layers
- !
- CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
- itf,jtf,ktf, &
- 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)+.15*entr_rate
- if(cd(i,k).gt.1.0*entr_rate)cd(i,k)=1.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,'deep', &
- itf,jtf,ktf, &
- 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,'deep', &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !--- DETERMINE CLOUD TOP - KTOP
- !
- call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- 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
- !
- 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
- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
- 1,kdet,z1, &
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- 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,&
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !--- calculate moisture properties of downdraft
- !
- call cup_dd_moisture_3d(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,high_resolution, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_dd_moisture_3d(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,high_resolution, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !--- calculate moisture properties of updraft
- !
- call cup_up_moisture('deep',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, &
- itf,jtf,ktf, &
- 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('deep',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,&
- itf,jtf,ktf, &
- 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, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
- kbcon,ktop,ierr, &
- itf,jtf,ktf, &
- 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
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! NEXT section for shallow convection
- !
- if(ishallow_g3.eq.1)then
- ! write(0,*)'now do shallow for j = ',j
- call cup_env(z3,qes3,he3,hes3,tshall,qshall,po,z1, &
- psur,ierr5,tcrit,0,xl,cp, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_env_clev(tshall,qes3,qshall,he3,hes3,z3,po,qes3_cup,q3_cup, &
- he3_cup,hes3_cup,z3_cup,po_cup,gamma3_cup,t3_cup,psur, &
- ierr5,z1,xl,rv,cp, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- CALL cup_MAXIMI(HE3_CUP,1,kbmax,K23,ierr5, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- DO i=its,itf
- if(kpbl(i).gt.5)cap_max3(i)=po_cup(i,kpbl(i))
- IF(ierr5(I).eq.0.)THEN
- IF(K23(I).Gt.Kbmax(i))ierr5(i)=2
- if(kpbl(i).gt.5)k23(i)=kpbl(i)
- endif
- ierr5_0(i)=ierr5(i)
- ENDDO
- call cup_kbcon(cap_max_increment,5,k23,kbcon3,he3_cup,hes3_cup, &
- ierr5,kbmax,po_cup,cap_max3, &
- ! ierr5,kpbl,po_cup,cap_max3, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_he(k23,hkb3,z3_cup,cd3,mentr_rate3,he3_cup,hc3, &
- kbcon3,ierr5,dby3,he3,hes3_cup,'shallow', &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_he(k23,hkb3_0,z_cup,cd3,mentr_rate3,he_cup,hc3_0, &
- kbcon3,ierr5,dby3_0,he,hes_cup,'shallow', &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_ktop(1,dby3,kbcon3,ktop3,ierr5, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_nms(zu3,z3_cup,mentr_rate3,cd3,kbcon3,ktop3, &
- ierr5,k23, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_nms(zu3_0,z_cup,mentr_rate3,cd3,kbcon3,ktop3, &
- ierr5,k23, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- ! first calculate aa3_0_cup
- !
- call cup_up_aa0(aa3_0,z,zu3_0,dby3_0,GAMMA3_0_CUP,t_cup, &
- kbcon3,ktop3,ierr5, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- ! now what is necessary for aa3 and feedbacks
- !
- call cup_up_moisture('shallow',ierr5,z3_cup,qc3,qrc3,pw3,pwav3, &
- kbcon3,ktop3,cd3,dby3,mentr_rate3,clw_all, &
- qshall,GAMMA3_cup,zu3,qes3_cup,k23,q3_cup,xl,&
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_aa0(aa3,z3,zu3,dby3,GAMMA3_CUP,t3_cup, &
- kbcon3,ktop3,ierr5, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- ! do i=its,itf
- ! if(ierr5(i).eq.0)then
- ! if(aa3(i).eq.0.)then
- ! ierr5(i)=17
- ! endif
- ! endif
- ! enddo
- ! call cup_dellabot('shallow',ipr,jpr,q3_cup,ierr5,z3_cup,po,qrcdo,edto, &
- ! zdo,cdd,q3,dellaq3,dsubq,j,mentrd_rate,z3,g,&
- ! itf,jtf,ktf, &
- ! its,ite, jts,jte, kts,kte)
- call cup_dellas_3d(ierr5,z3_cup,po_cup,hcdo,edt3,zdo3,cdd, &
- he3,dellah3,dsubt3,j,mentrd_rate,zu3,g, &
- cd3,hc3,ktop3,k23,kbcon3,mentr_rate3,jmin,he3_cup,kdet, &
- k23,ipr,jpr,'shallow',0, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_dellas_3d(ierr5,z3_cup,po_cup,qrcdo,edt3,zdo3,cdd, &
- qshall,dellaq3,dsubq3,j,mentrd_rate,zu3,g, &
- cd3,qc3,ktop3,k23,kbcon3,mentr_rate3,jmin,q3_cup,kdet, &
- k23,ipr,jpr,'shallow',0, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte )
- mbdt_s=1.e-1*mbdt_ens(1)
- do k=kts,ktf
- do i=its,itf
- dellat3(i,k)=0.
- if(ierr5(i).eq.0)then
- trash=dsubt3(i,k)
- XHE3(I,K)=(dsubt3(i,k)+DELLAH3(I,K))*MBDT_S+HE3(I,K)
- XQ3(I,K)=(dsubq3(i,k)+DELLAQ3(I,K))*MBDT_S+QSHALL(I,K)
- DELLAT3(I,K)=(1./cp)*(DELLAH3(I,K)-xl*DELLAQ3(I,K))
- dSUBT3(I,K)=(1./cp)*(dsubt3(i,k)-xl*dsubq3(i,k))
- XT3(I,K)= (DELLAT3(I,K)+dsubt3(i,k))*MBDT_S+TSHALL(I,K)
- IF(XQ3(I,K).LE.0.)XQ3(I,K)=1.E-08
- ! if(i.eq.ipr.and.j.eq.jpr)then
- ! write(0,*)k,trash,DELLAQ3(I,K),dsubq3(I,K),dsubt3(i,k)
- ! endif
- ENDIF
- enddo
- enddo
- do i=its,itf
- if(ierr5(i).eq.0)then
- XHE3(I,ktf)=HE3(I,ktf)
- XQ3(I,ktf)=QSHALL(I,ktf)
- XT3(I,ktf)=TSHALL(I,ktf)
- IF(XQ3(I,ktf).LE.0.)XQ3(I,ktf)=1.E-08
- endif
- enddo
- !
- !--- calculate moist static energy, heights, qes
- !
- call cup_env(xz3,xqes3,xhe3,xhes3,xt3,xq3,po,z1, &
- psur,ierr5,tcrit,2,xl,cp, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !--- environmental values on cloud levels
- !
- call cup_env_clev(xt3,xqes3,xq3,xhe3,xhes3,xz3,po,xqes3_cup,xq3_cup, &
- xhe3_cup,xhes3_cup,xz3_cup,po_cup,gamma3_cup,xt3_cup,psur, &
- ierr5,z1,xl,rv,cp, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !
- !**************************** static control
- !
- !--- moist static energy inside cloud
- !
- do i=its,itf
- if(ierr5(i).eq.0)then
- xhkb3(i)=xhe3(i,k23(i))
- endif
- enddo
- call cup_up_he(k23,xhkb3,xz3_cup,cd3,mentr_rate3,xhe3_cup,xhc3, &
- kbcon3,ierr5,xdby3,xhe3,xhes3_cup,'shallow', &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- !c--- normalized mass flux profile and CWF
- !
- call cup_up_nms(xzu3,xz3_cup,mentr_rate3,cd3,kbcon3,ktop3,ierr5,k23, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- call cup_up_aa0(xaa3,xz3,xzu3,xdby3,GAMMA3_CUP,xt3_cup, &
- kbcon3,ktop3,ierr5, &
- itf,jtf,ktf, &
- its,ite, jts,jte, kts,kte)
- !
- ! now for shallow forcing
- !
- do i=its,itf
- xmb3(i)=0.
- xff_shal(1:9)=0.
- if(ierr5(i).eq.0)then
- xkshal=(xaa3(i)-aa3(i))/mbdt_s
- if(xkshal.ge.0.)xkshal=+1.e6
- if(xkshal.gt.-1.e-4 .and. xkshal.lt.0.)xkshal=-1.e-4
- xff_shal(1)=max(0.,-(aa3(i)-aa3_0(i))/(xkshal*dtime))
- xff_shal(2)=max(0.,-(aa3(i)-aa3_…
Large files files are truncated, but you can click here to view the full file