/wrfv2_fire/phys/module_cu_kfeta.F
FORTRAN Legacy | 3174 lines | 2002 code | 116 blank | 1056 comment | 19 complexity | 83c977654d6529edf366d4c5881efc42 MD5 | raw file
Possible License(s): AGPL-1.0
- MODULE module_cu_kfeta
- USE module_wrf_error
- !
- ! V3.3: A new trigger function is added based Ma and Tan (2009):
- ! Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of
- ! the cumulus parameterization for tropical cyclone prediction:
- ! Convection trigger. Atmospheric Research, 92, 190 - 211.
- !
- !--------------------------------------------------------------------
- ! Lookup table variables:
- INTEGER, PARAMETER :: KFNT=250,KFNP=220
- REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
- REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
- REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
- REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
- ! Note: KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
- ! TPMIX2DD, ENVIRTHT
- ! End of Lookup table variables:
- CONTAINS
- SUBROUTINE KF_eta_CPS( &
- ids,ide, jds,jde, kds,kde &
- ,ims,ime, jms,jme, kms,kme &
- ,its,ite, jts,jte, kts,kte &
- ,trigger &
- ,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
- ,CUDTACTTIME &
- ,rho,RAINCV,PRATEC,NCA &
- ,U,V,TH,T,W,dz8w,Pcps,pi &
- ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
- ,EP2,SVP1,SVP2,SVP3,SVPT0 &
- ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
- ,QV &
- ! optionals
- ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
- ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
- ,RQICUTEN,RQSCUTEN, RQVFTEN &
- )
- !
- !-------------------------------------------------------------
- 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 ) :: trigger
- INTEGER, INTENT(IN ) :: STEPCU
- LOGICAL, INTENT(IN ) :: warm_rain
- REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
- REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
- REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
- INTEGER, INTENT(IN ) :: KTAU
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
- INTENT(IN ) :: &
- U, &
- V, &
- W, &
- TH, &
- T, &
- QV, &
- dz8w, &
- Pcps, &
- rho, &
- pi
- !
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
- INTENT(INOUT) :: &
- W0AVG
- REAL, INTENT(IN ) :: DT, DX
- REAL, INTENT(IN ) :: CUDT
- REAL, INTENT(IN ) :: CURR_SECS
- LOGICAL,OPTIONAL,INTENT(IN ) :: ADAPT_STEP_FLAG
- REAL, INTENT (INOUT) :: CUDTACTTIME
- !
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: RAINCV
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: PRATEC
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: NCA
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(OUT) :: CUBOT, &
- CUTOP
- LOGICAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: CU_ACT_FLAG
- !
- ! Optional arguments
- !
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
- OPTIONAL, &
- INTENT(INOUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQCCUTEN, &
- RQRCUTEN, &
- RQICUTEN, &
- RQSCUTEN, &
- RQVFTEN
- !
- ! 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
- ! LOCAL VARS
- LOGICAL :: flag_qr, flag_qi, flag_qs
- REAL, DIMENSION( kts:kte ) :: &
- U1D, &
- V1D, &
- T1D, &
- DZ1D, &
- QV1D, &
- P1D, &
- RHO1D, &
- tpart_v1D, &
- tpart_h1D, &
- W0AVG1D
- REAL, DIMENSION( kts:kte ):: &
- DQDT, &
- DQIDT, &
- DQCDT, &
- DQRDT, &
- DQSDT, &
- DTDT
- REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) :: aveh_t, aveh_q
- REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
- REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_t, avev_q
- REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
- REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: coef_v, coef_h, tpart_h, tpart_v
- INTEGER :: ii,jj,kk
- REAL :: ttop
- REAL, DIMENSION (kts:kte) :: z0
- REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
- integer :: ibegh,iendh,jbegh,jendh
- integer :: istart,iend,jstart,jend
- INTEGER :: i,j,k,NTST
- REAL :: lastdt = -1.0
- REAL :: W0AVGfctr, W0fctr, W0den
- LOGICAL :: run_param , doing_adapt_dt , decided
-
- !
- DXSQ=DX*DX
- !----------------------
- NTST=STEPCU
- TST=float(NTST*2)
- flag_qr = .FALSE.
- flag_qi = .FALSE.
- flag_qs = .FALSE.
- IF ( PRESENT(F_QR) ) flag_qr = F_QR
- IF ( PRESENT(F_QI) ) flag_qi = F_QI
- IF ( PRESENT(F_QS) ) flag_qs = F_QS
- !
- if (lastdt < 0) then
- lastdt = dt
- endif
-
- if (ADAPT_STEP_FLAG) then
- W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
- W0fctr = dt
- W0den = 2 * MAX(CUDT*60,dt)
- else
- W0AVGfctr = (TST-1.)
- W0fctr = 1.
- W0den = TST
- endif
- DO J = jts,jte
- DO K=kts,kte
- DO I= its,ite
- ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
- ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
- ! RHOE=Pcps(I,K,J)/(R*TV)
- ! W0=-101.9368*SCR1/RHOE
- W0=0.5*(w(I,K,J)+w(I,K+1,J))
- ! Old:
- !
- ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
- !
- ! New, to support adaptive time step:
- !
- W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
- ENDDO
- ENDDO
- ENDDO
- lastdt = dt
- ! Initialization for adaptive time step.
- doing_adapt_dt = .FALSE.
- IF ( PRESENT(adapt_step_flag) ) THEN
- IF ( adapt_step_flag ) THEN
- doing_adapt_dt = .TRUE.
- IF ( cudtacttime .EQ. 0. ) THEN
- cudtacttime = curr_secs + cudt*60.
- END IF
- END IF
- END IF
- ! Do we run through this scheme or not?
- ! Test 1: If this is the initial model time, then yes.
- ! KTAU=1
- ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
- ! CUDT=0 or STEPCU=1
- ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
- ! MOD(KTAU,NST)=0
- ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
- ! CURR_SECS >= CUDTACTTIME
- ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
- ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
- ! We only proceed to other tests if the previous tests all have left decided as FALSE.
- ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
- ! cumulus run.
- decided = .FALSE.
- run_param = .FALSE.
- IF ( ( .NOT. decided ) .AND. &
- ( ktau .EQ. 1 ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( .NOT. doing_adapt_dt ) .AND. &
- ( MOD(ktau,ntst) .EQ. 0 ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- END IF
- IF ( ( .NOT. decided ) .AND. &
- ( doing_adapt_dt ) .AND. &
- ( curr_secs .GE. cudtacttime ) ) THEN
- run_param = .TRUE.
- decided = .TRUE.
- cudtacttime = curr_secs + cudt*60
- END IF
- if (run_param) then
- ! New trigger function
- IF (trigger.eq.2) THEN
- !
- ! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
- !
- aveh_t=-999 ! horizontal 9-point ave
- aveh_q=-999
- avev_t=0 ! vertical 3-level ave
- avev_q=0
- avev_qmax=0
- avev_qmin=0
- aveh_qmax=0
- aveh_qmin=0
- tpart_h=0
- tpart_v=0
- coef_h=0
- coef_v=0
- ibegh=max(its-1, ids+1) ! start from 2
- jbegh=max(jts-1, jds+1)
- iendh=min(ite+1, ide-2) ! end at ide-2
- jendh=min(jte+1, jde-2)
- DO J = jbegh,jendh
- DO K = kts,kte
- DO I = ibegh,iendh
- aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j) +T(i-1,k,j+1)+ &
- T(i,k,j-1) +T(i,k,j) +T(i,k,j+1)+ &
- T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
- aveh_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.
- ENDDO
- ENDDO
- ENDDO
- ! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
- DO K = kts,kte
- DO J = jts-1,jte+1
- DO I = its-1,ite+1
- if(i.eq.ids) then
- aveh_t(i,k,j)=aveh_t(i+1,k,j)
- aveh_q(i,k,j)=aveh_q(i+1,k,j)
- elseif(i.eq.ide-1) then
- aveh_t(i,k,j)=aveh_t(i-1,k,j)
- aveh_q(i,k,j)=aveh_q(i-1,k,j)
- endif
- if(j.eq.jds) then
- aveh_t(i,k,j)=aveh_t(i,k,j+1)
- aveh_q(i,k,j)=aveh_q(i,k,j+1)
- elseif(j.eq.jde-1) then
- aveh_t(i,k,j)=aveh_t(i,k,j-1)
- aveh_q(i,k,j)=aveh_q(i,k,j-1)
- endif
- if(j.eq.jds.and.i.eq.ids) then
- aveh_q(i,k,j)=aveh_q(i+1,k,j+1)
- aveh_t(i,k,j)=aveh_t(i+1,k,j+1)
- endif
- if(j.eq.jde-1.and.i.eq.ids) then
- aveh_q(i,k,j)=aveh_q(i+1,k,j-1)
- aveh_t(i,k,j)=aveh_t(i+1,k,j-1)
- endif
- if(j.eq.jde-1.and.i.eq.ide-1) then
- aveh_q(i,k,j)=aveh_q(i-1,k,j-1)
- aveh_t(i,k,j)=aveh_t(i-1,k,j-1)
- endif
- if(j.eq.jds.and.i.eq.ide-1) then
- aveh_q(i,k,j)=aveh_q(i-1,k,j+1)
- aveh_t(i,k,j)=aveh_t(i-1,k,j+1)
- endif
- ENDDO
- ENDDO
- ENDDO
- ! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
- istart=max(its, ids+1) ! start from 2
- jstart=max(jts, jds+1)
- iend=min(ite, ide-2) ! end at ide-2
- jend=min(jte, jde-2)
- DO K = kts,kte
- DO J = jstart,jend
- DO I = istart,iend
- aveh_qmax(i,k,j)=aveh_q(i,k,j)
- aveh_qmin(i,k,j)=aveh_q(i,k,j)
- DO ii=-1, 1
- DO jj=-1,1
- if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
- if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
- ENDDO
- ENDDO
- if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
- coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
- else
- coef_h(i,k,j)=0.
- endif
- coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
- coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
- tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
- ENDDO
- ENDDO
- ENDDO
- 89 continue
- ! vertical 3-layer calculation
- DO J = jts, jte
- DO I = its, ite
- z0(1) = 0.5 * dz8w(i,1,j)
- DO K = 2, kte
- Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
- ENDDO
- DO K = kts+1,kte-1
- ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
- avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
- ! avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
- avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
- ENDDO
- avev_t(i,kts,j)=avev_t(i,kts+1,j) ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
- avev_q(i,kts,j)=avev_q(i,kts+1,j)
- avev_t(i,kte,j)=avev_t(i,kte-1,j) ! highest level value
- avev_q(i,kte,j)=avev_q(i,kte-1,j)
- ENDDO
- ENDDO
- ! max /min value
- DO J = jts, jte
- DO I = its, ite
- DO K = kts+1,kte-1
- avev_qmax(i,k,j)=avev_q(i,k,j)
- avev_qmin(i,k,j)=avev_q(i,k,j)
- DO kk=-1,1
- if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j)
- if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j)
- ENDDO
- if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
- coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
- else
- coef_v(i,k,j)=0
- endif
- tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
- ENDDO
- tpart_v(i,kts,j)= tpart_v(i,kts+1,j) ! lowest level
- tpart_v(i,kte,j)= tpart_v(i,kte-1,j) ! highest level
- ENDDO
- ENDDO
- ENDIF ! endif (trigger.eq.2)
- !
- DO J = jts,jte
- DO I= its,ite
- CU_ACT_FLAG(i,j) = .true.
- ENDDO
- ENDDO
- DO J = jts,jte
- DO I=its,ite
-
- IF ( NCA(I,J) .ge. 0.5*DT ) then
- CU_ACT_FLAG(i,j) = .false.
- ELSE
- DO k=kts,kte
- DQDT(k)=0.
- DQIDT(k)=0.
- DQCDT(k)=0.
- DQRDT(k)=0.
- DQSDT(k)=0.
- DTDT(k)=0.
- ENDDO
- RAINCV(I,J)=0.
- CUTOP(I,J)=KTS
- CUBOT(I,J)=KTE+1
- PRATEC(I,J)=0.
- !
- ! assign vars from 3D to 1D
- DO K=kts,kte
- U1D(K) =U(I,K,J)
- V1D(K) =V(I,K,J)
- T1D(K) =T(I,K,J)
- RHO1D(K) =rho(I,K,J)
- QV1D(K)=QV(I,K,J)
- P1D(K) =Pcps(I,K,J)
- W0AVG1D(K) =W0AVG(I,K,J)
- DZ1D(k)=dz8w(I,K,J)
- IF (trigger.eq.2) THEN
- tpart_h1D(K) =tpart_h(I,K,J)
- tpart_v1D(K) =tpart_v(I,K,J)
- ELSE
- tpart_h1D(K) = 0.
- tpart_v1D(K) = 0.
- ENDIF
- ENDDO
- CALL KF_eta_PARA(I, J, &
- U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
- tpart_h1D,tpart_v1D, &
- trigger, &
- DT,DX,DXSQ,RHO1D, &
- XLV0,XLV1,XLS0,XLS1,CP,R,G, &
- EP2,SVP1,SVP2,SVP3,SVPT0, &
- DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
- RAINCV,PRATEC,NCA, &
- flag_QI,flag_QS,warm_rain, &
- CUTOP,CUBOT,CUDT, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
- DO K=kts,kte
- RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
- RQVCUTEN(I,K,J)=DQDT(K)
- ENDDO
- ENDIF
- IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
- IF( F_QR )THEN
- DO K=kts,kte
- RQRCUTEN(I,K,J)=DQRDT(K)
- RQCCUTEN(I,K,J)=DQCDT(K)
- ENDDO
- ELSE
- ! This is the case for Eta microphysics without 3d rain field
- DO K=kts,kte
- RQRCUTEN(I,K,J)=0.
- RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
- ENDDO
- ENDIF
- ENDIF
- !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
- IF(PRESENT( rqicuten )) THEN
- IF ( F_QI ) THEN
- DO K=kts,kte
- RQICUTEN(I,K,J)=DQIDT(K)
- ENDDO
- ENDIF
- ENDIF
- IF(PRESENT( rqscuten )) THEN
- IF ( F_QS ) THEN
- DO K=kts,kte
- RQSCUTEN(I,K,J)=DQSDT(K)
- ENDDO
- ENDIF
- ENDIF
- !
- ENDIF
- ENDDO ! i-loop
- ENDDO ! j-loop
- ENDIF ! run_param
- !
- END SUBROUTINE KF_eta_CPS
- ! ****************************************************************************
- !-----------------------------------------------------------
- SUBROUTINE KF_eta_PARA (I, J, &
- U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
- TPART_H0,TPART_V0, &
- trigger, &
- DT,DX,DXSQ,rhoe, &
- XLV0,XLV1,XLS0,XLS1,CP,R,G, &
- EP2,SVP1,SVP2,SVP3,SVPT0, &
- DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
- RAINCV,PRATEC,NCA, &
- F_QI,F_QS,warm_rain, &
- CUTOP,CUBOT,CUDT, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !-----------------------------------------------------------
- !***** The KF scheme that is currently used in experimental runs of EMCs
- !***** Eta model....jsk 8/00
- !
- IMPLICIT NONE
- !-----------------------------------------------------------
- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte, &
- I,J
- ! ,P_QI,P_QS,P_FIRST_SCALAR
- INTEGER, INTENT(IN ) :: trigger
- LOGICAL, INTENT(IN ) :: F_QI, F_QS
- LOGICAL, INTENT(IN ) :: warm_rain
- !
- REAL, DIMENSION( kts:kte ), &
- INTENT(IN ) :: U0, &
- V0, &
- TPART_H0, &
- TPART_V0, &
- T0, &
- QV0, &
- P0, &
- rhoe, &
- DZQ, &
- W0AVG1D
- !
- REAL, INTENT(IN ) :: DT,DX,DXSQ
- !
- REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
- REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
- !
- REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
- DQDT, &
- DQIDT, &
- DQCDT, &
- DQRDT, &
- DQSDT, &
- DTDT
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: NCA
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: RAINCV
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: PRATEC
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(OUT) :: CUBOT, &
- CUTOP
- REAL, INTENT(IN ) :: CUDT
- !
- !...DEFINE LOCAL VARIABLES...
- !
- REAL, DIMENSION( kts:kte ) :: &
- Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
- QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
- UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
- UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
- THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
- QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
- DETLQ2,DETIC2,RATIO,RATIO2
- REAL, DIMENSION( kts:kte ) :: &
- DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
- QDT,FXM,THTAG,THPA,THFXOUT, &
- THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
- QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
- QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
- QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
- REAL, DIMENSION( kts:kte+1 ) :: OMG
- REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
- REAL, DIMENSION( kts:kte ) :: &
- CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
- ! LOCAL VARS
- REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
- TTFRZ,TBFRZ,C5,RATE
- REAL :: GDRY,ROCP,ALIQ,BLIQ, &
- CLIQ,DLIQ
- REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
- ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
- CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
- ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
- TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
- UPNEW,ABE,WKLCL,TTEMP,FRC1, &
- QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
- DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
- THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
- UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
- THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
- CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
- DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
- DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
- UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
- DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
- AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
- DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
- TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
- UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
- RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
- DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
- REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
- QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
- !
- INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
- REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
- REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
- INTEGER :: KX,K,KL
- !
- INTEGER :: NCHECK
- INTEGER, DIMENSION (kts:kte) :: KCHECK
- INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
- LC,MXLAYR,LLFC,NLAYRS,NK, &
- KPBL,KLCL,LCL,LET,IFLAG, &
- NK1,LTOP,NJ,LTOP1, &
- LTOPM1,LVF,KSTART,KMIN,LFS, &
- ND,NIC,LDB,LDT,ND1,NDK, &
- NM,LMAX,NCOUNT,NOITR, &
- NSTEP,NTC,NCHM,ISHALL,NSHALL
- LOGICAL :: IPRNT
- REAL :: u00,qslcl,rhlcl,dqssdt !jfb
- CHARACTER*1024 message
- !
- DATA P00,T00/1.E5,273.16/
- DATA RLF/3.339E5/
- DATA RHIC,RHBC/1.,0.90/
- DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
- DATA RATE/0.03/ ! wrf default
- ! DATA RATE/0.01/ ! value used in NRCM
- ! DATA RATE/0.001/ ! effectively turn off autoconversion
- !-----------------------------------------------------------
- IPRNT=.FALSE.
- GDRY=-G/CP
- ROCP=R/CP
- NSHALL = 0
- KL=kte
- KX=kte
- !
- ! ALIQ = 613.3
- ! BLIQ = 17.502
- ! CLIQ = 4780.8
- ! DLIQ = 32.19
- ALIQ = SVP1*1000.
- BLIQ = SVP2
- CLIQ = SVP2*SVPT0
- DLIQ = SVP3
- !
- !
- !****************************************************************************
- ! ! PPT FB MODS
- !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
- !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
- !...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
- !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
- FBFRC=0.0 ! PPT FB MODS
- !...mods to allow shallow convection...
- NCHM = 0
- ISHALL = 0
- DPMIN = 5.E3
- !...
- P300=P0(1)-30000.
- !
- !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
- !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
- !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
- !
- !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
- !...FROM BOTTOM-UP IN THE KF SCHEME...
- !
- ML=0
- !SUE tmprpsb=1./PSB(I,J)
- !SUE CELL=PTOP*tmprpsb
- !
- DO K=1,KX
- !
- ! Saturation vapor pressure (ES) is calculated following Buck (1981)
- !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
- !
- ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
- QES(K)=0.622*ES/(P0(K)-ES)
- Q0(K)=AMIN1(QES(K),QV0(K))
- Q0(K)=AMAX1(0.000001,Q0(K))
- QL0(K)=0.
- QI0(K)=0.
- QR0(K)=0.
- QS0(K)=0.
- RH(K) = Q0(K)/QES(K)
- DILFRC(K) = 1.
- TV0(K)=T0(K)*(1.+0.608*Q0(K))
- ! RHOE(K)=P0(K)/(R*TV0(K))
- ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
- DP(K)=rhoe(k)*g*DZQ(k)
- ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
- ! use it for shallow convection...For now, assume it is not available....
- ! TKE(K) = Q2(I,J,NK)
- TKE(K) = 0.
- CLDHGT(K) = 0.
- ! IF(P0(K).GE.500E2)L5=K
- IF(P0(K).GE.0.5*P0(1))L5=K
- IF(P0(K).GE.P300)LLFC=K
- ENDDO
- !
- !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
- Z0(1)=.5*DZQ(1)
- !cdir novector
- DO K=2,KL
- Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
- DZA(K-1)=Z0(K)-Z0(K-1)
- ENDDO
- DZA(KL)=0.
- !
- !
- ! To save time, specify a pressure interval to move up in sequential
- ! check of different ~50 mb deep groups of adjacent model layers in
- ! the process of identifying updraft source layer (USL). Note that
- ! this search is terminated as soon as a buoyant parcel is found and
- ! this parcel can produce a cloud greater than specifed minimum depth
- ! (CHMIN)...For now, set interval at 15 mb...
- !
- NCHECK = 1
- KCHECK(NCHECK)=1
- PM15 = P0(1)-15.E2
- DO K=2,LLFC
- IF(P0(K).LT.PM15)THEN
- NCHECK = NCHECK+1
- KCHECK(NCHECK) = K
- PM15 = PM15-15.E2
- ENDIF
- ENDDO
- !
- NU=0
- NUCHM=0
- usl: DO
- NU = NU+1
- IF(NU.GT.NCHECK)THEN
- IF(ISHALL.EQ.1)THEN
- CHMAX = 0.
- NCHM = 0
- DO NK = 1,NCHECK
- NNN=KCHECK(NK)
- IF(CLDHGT(NNN).GT.CHMAX)THEN
- NCHM = NNN
- NUCHM = NK
- CHMAX = CLDHGT(NNN)
- ENDIF
- ENDDO
- NU = NUCHM-1
- FBFRC=1.
- CYCLE usl
- ELSE
- RETURN
- ENDIF
- ENDIF
- KMIX = KCHECK(NU)
- LOW=KMIX
- !...
- LC = LOW
- !
- !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
- !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
- !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
- !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
- !
- NLAYRS=0
- DPTHMX=0.
- NK=LC-1
- IF ( NK+1 .LT. KTS ) THEN
- WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
- CALL wrf_message (TRIM(message))
- ELSE
- DO
- NK=NK+1
- IF ( NK .GT. KTE ) THEN
- WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
- CALL wrf_message (TRIM(message))
- EXIT
- ENDIF
- DPTHMX=DPTHMX+DP(NK)
- NLAYRS=NLAYRS+1
- IF(DPTHMX.GT.DPMIN)THEN
- EXIT
- ENDIF
- END DO
- ENDIF
- IF(DPTHMX.LT.DPMIN)THEN
- RETURN
- ENDIF
- KPBL=LC+NLAYRS-1
- !
- !...********************************************************
- !...for computational simplicity without much loss in accuracy,
- !...mix temperature instead of theta for evaluating convective
- !...initiation (triggering) potential...
- ! THMIX=0.
- TMIX=0.
- QMIX=0.
- ZMIX=0.
- PMIX=0.
- !
- !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
- !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
- !...LAYERS...
- !
- !cdir novector
- DO NK=LC,KPBL
- TMIX=TMIX+DP(NK)*T0(NK)
- QMIX=QMIX+DP(NK)*Q0(NK)
- ZMIX=ZMIX+DP(NK)*Z0(NK)
- PMIX=PMIX+DP(NK)*P0(NK)
- ENDDO
- ! THMIX=THMIX/DPTHMX
- TMIX=TMIX/DPTHMX
- QMIX=QMIX/DPTHMX
- ZMIX=ZMIX/DPTHMX
- PMIX=PMIX/DPTHMX
- EMIX=QMIX*PMIX/(0.622+QMIX)
- !
- !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
- !
- ! TLOG=ALOG(EMIX/ALIQ)
- ! ...calculate dewpoint using lookup table...
- !
- astrt=1.e-3
- ainc=0.075
- a1=emix/aliq
- tp=(a1-astrt)/ainc
- indlu=int(tp)+1
- value=(indlu-1)*ainc+astrt
- aintrp=(a1-value)/ainc
- tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
- TLCL=AMIN1(TLCL,TMIX)
- TVLCL=TLCL*(1.+0.608*QMIX)
- ZLCL = ZMIX+(TLCL-TMIX)/GDRY
- NK = LC-1
- DO
- NK = NK+1
- KLCL=NK
- IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
- EXIT
- ENDIF
- ENDDO
- IF(NK.GT.KL)THEN
- RETURN
- ENDIF
- K=KLCL-1
- ! calculate DLP using Z instead of log(P)
- DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
- !
- !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
- !
- TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
- QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
- TVEN=TENV*(1.+0.608*QENV)
- !
- !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
- !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
- !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
- !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
- !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
- !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
- !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
- !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
- !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
- !
- IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
- WKLCL=0.02*ZLCL/2.E3
- ELSE
- WKLCL=0.02 ! units of m/s
- ENDIF
- WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
- IF(WKL.LT.0.0001)THEN
- DTLCL=0.
- ELSE
- DTLCL=4.64*WKL**0.33 ! Kain (2004) Eq. 1
- ENDIF
- ! New trigger function
- IF(trigger.eq.2) then
- DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
- ENDIF
- ! end for trigger function
- !
- dtrh = 0.
- if (trigger .eq. 3) then
- !...for ETA model, give parcel an extra temperature perturbation based
- !...the threshold RH for condensation (U00)...
- ! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
- !
- !...for now, just assume U00=0.75...
- !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
- U00 = 0.75
- IF(U00.lt.1.)THEN
- QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
- RHLCL = QENV/QSLCL
- DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
- IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
- DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
- ELSEIF(RHLCL.GT.0.95)THEN
- DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
- ELSE
- DTRH = 0.
- ENDIF
- ENDIF
- endif ! trigger 3
- ! IF(ISHALL.EQ.1)IPRNT=.TRUE.
- ! IPRNT=.TRUE.
- ! IF(TLCL+DTLCL.GT.TENV)GOTO 45
-
- trigger2: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
- !
- ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
- !
- CYCLE usl
- !
- ELSE ! Parcel is buoyant, determine updraft
- !
- !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
- !...EQUIVALENT POTENTIAL TEMPERATURE
- !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
- !
- CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
- !
- !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
- !
- DTTOT = DTLCL+DTRH
- IF(DTTOT.GT.1.E-4)THEN
- GDT=2.*G*DTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
- WLCL=1.+0.5*SQRT(GDT)
- WLCL = AMIN1(WLCL,3.)
- ELSE
- WLCL=1.
- ENDIF
- PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
- WTW=WLCL*WLCL
- !
- TVLCL=TLCL*(1.+0.608*QMIX)
- RHOLCL=PLCL/(R*TVLCL)
- !
- LCL=KLCL
- LET=LCL
- ! make RAD a function of background vertical velocity... (Kain (2004) Eq. 6)
- IF(WKL.LT.0.)THEN
- RAD = 1000.
- ELSEIF(WKL.GT.0.1)THEN
- RAD = 2000.
- ELSE
- RAD = 1000.+1000*WKL/0.1
- ENDIF
- !
- !*******************************************************************
- ! *
- ! COMPUTE UPDRAFT PROPERTIES *
- ! *
- !*******************************************************************
- !
- !
- !...
- !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
- !
- WU(K)=WLCL
- AU0=0.01*DXSQ
- UMF(K)=RHOLCL*AU0
- VMFLCL=UMF(K)
- UPOLD=VMFLCL
- UPNEW=UPOLD
- !
- !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
- !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
- !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
- !...PRODUCTION...
- !
- RATIO2(K)=0.
- UER(K)=0.
- ABE=0.
- TRPPT=0.
- TU(K)=TLCL
- TVU(K)=TVLCL
- QU(K)=QMIX
- EQFRC(K)=1.
- QLIQ(K)=0.
- QICE(K)=0.
- QLQOUT(K)=0.
- QICOUT(K)=0.
- DETLQ(K)=0.
- DETIC(K)=0.
- PPTLIQ(K)=0.
- PPTICE(K)=0.
- IFLAG=0
- !
- !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
- !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
- !...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
- !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
- !...PREVIOUS MODEL LEVEL...
- !
- TTEMP=TTFRZ
- !
- !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
- !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
- !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
- !
- ! **1 variables indicate the bottom of a model layer
- ! **2 variables indicate the top of a model layer
- !
- EE1=1.
- UD1=0.
- REI = 0.
- DILBE = 0.
- updraft: DO NK=K,KL-1
- NK1=NK+1
- RATIO2(NK1)=RATIO2(NK)
- FRC1=0.
- TU(NK1)=T0(NK1)
- THETEU(NK1)=THETEU(NK)
- QU(NK1)=QU(NK)
- QLIQ(NK1)=QLIQ(NK)
- QICE(NK1)=QICE(NK)
- call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
- qice(nk1),qnewlq,qnewic,XLV1,XLV0)
- !
- !
- !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
- !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
- !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
- !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
- !...LIQUID WATER IS FROZEN AT EACH LEVEL...
- !
- IF(TU(NK1).LE.TTFRZ)THEN
- IF(TU(NK1).GT.TBFRZ)THEN
- IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
- FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
- ELSE
- FRC1=1.
- IFLAG=1
- ENDIF
- TTEMP=TU(NK1)
- !
- ! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
- !...IS BELOW TTFRZ...
- !
- QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
- QNEWIC=QNEWIC+QNEWLQ*FRC1
- QNEWLQ=QNEWLQ-QNEWLQ*FRC1
- QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
- QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
- CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
- QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
- ENDIF
- TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
- !
- ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
- !
- IF(NK.EQ.K)THEN
- BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
- BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
- DZZ=Z0(NK1)-ZLCL
- ELSE
- BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
- BOTERM=2.*DZA(NK)*G*BE/1.5
- DZZ=DZA(NK)
- ENDIF
- ENTERM=2.*REI*WTW/UPOLD
- CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
- RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
- !
- !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
- !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
- !
- IF(WTW.LT.1.E-3)THEN
- EXIT
- ELSE
- WU(NK1)=SQRT(WTW)
- ENDIF
- !...Calculate value of THETA-E in environment to entrain into updraft...
- !
- CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
- !
- !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
- !
- REI=VMFLCL*DP(NK1)*0.03/RAD ! KF (1990) Eq. 1; Kain (2004) Eq. 5
- TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
- IF(NK.EQ.K)THEN
- DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
- ELSE
- DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
- ENDIF
- IF(DILBE.GT.0.)ABE=ABE+DILBE*G
- !
- !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
- !...ENTRAINMENT (0.5*REI) IS IMPOSED...
- !
- IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
- EE2=0.5 ! Kain (2004) Eq. 4
- UD2=1.
- EQFRC(NK1)=0.
- ELSE
- LET=NK1
- TTMP=TVQU(NK1)
- !
- !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
- !
- F1=0.95
- F2=1.-F1
- THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
- QTMP=F1*Q0(NK1)+F2*QU(NK1)
- TMPLIQ=F2*QLIQ(NK1)
- TMPICE=F2*QICE(NK1)
- call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
- qnewlq,qnewic,XLV1,XLV0)
- TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
- IF(TU95.GT.TV0(NK1))THEN
- EE2=1.
- UD2=0.
- EQFRC(NK1)=1.0
- ELSE
- F1=0.10
- F2=1.-F1
- THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
- QTMP=F1*Q0(NK1)+F2*QU(NK1)
- TMPLIQ=F2*QLIQ(NK1)
- TMPICE=F2*QICE(NK1)
- call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
- qnewlq,qnewic,XLV1,XLV0)
- TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
- TVDIFF = ABS(TU10-TVQU(NK1))
- IF(TVDIFF.LT.1.e-3)THEN
- EE2=1.
- UD2=0.
- EQFRC(NK1)=1.0
- ELSE
- EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
- EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
- EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
- IF(EQFRC(NK1).EQ.1)THEN
- EE2=1.
- UD2=0.
- ELSEIF(EQFRC(NK1).EQ.0.)THEN
- EE2=0.
- UD2=1.
- ELSE
- !
- !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
- ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
- !
- CALL PROF5(EQFRC(NK1),EE2,UD2)
- ENDIF
- ENDIF
- ENDIF
- ENDIF ! End of Entrain/Detrain IF BLOCK
- !
- !
- !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
- ! VALUES IN THE LAYER...
- !
- EE2 = AMAX1(EE2,0.5)
- UD2 = 1.5*UD2
- UER(NK1)=0.5*REI*(EE1+EE2)
- UDR(NK1)=0.5*REI*(UD1+UD2)
- !
- !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
- ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
- !
- IF(UMF(NK)-UDR(NK1).LT.10.)THEN
- !
- !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
- ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
- ! First, correct ABE calculation if needed...
- !
- IF(DILBE.GT.0.)THEN
- ABE=ABE-DILBE*G
- ENDIF
- LET=NK
- ! WRITE(98,1015)P0(NK1)/100.
- EXIT
- ELSE
- EE1=EE2
- UD1=UD2
- UPOLD=UMF(NK)-UDR(NK1)
- UPNEW=UPOLD+UER(NK1)
- UMF(NK1)=UPNEW
- DILFRC(NK1) = UPNEW/UPOLD
- !
- !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
- !...ICE IN THE DETRAINING UPDRAFT MASS...
- !
- DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
- DETIC(NK1)=QICE(NK1)*UDR(NK1)
- QDT(NK1)=QU(NK1)
- QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
- THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
- QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
- QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
- !
- !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
- !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
- !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
- !...CURRENT MODEL LEVEL...
- !
- PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
- PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
- !
- TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
- IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
- ENDIF
- !
- END DO updraft
- !
- !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
- ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
- ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
- ! THE LET AND CLOUD TOP...
- !
- !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
- ! FIRST BECOMES NEGATIVE...
- !
- LTOP=NK
- CLDHGT(LC)=Z0(LTOP)-ZLCL
- !
- !...Instead of using the same minimum cloud height (for deep convection)
- !...everywhere, try specifying minimum cloud depth as a function of TLCL...
- !
- ! Kain (2004) Eq. 7
- !
- IF(TLCL.GT.293.)THEN
- CHMIN = 4.E3
- ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
- CHMIN = 2.E3 + 100.*(TLCL-273.)
- ELSEIF(TLCL.LT.273.)THEN
- CHMIN = 2.E3
- ENDIF
- !
- !...If cloud top height is less than the specified minimum for deep
- !...convection, save value to consider this level as source for
- !...shallow convection, go back up to check next level...
- !
- !...Try specifying minimum cloud depth as a function of TLCL...
- !
- !
- !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
- !
- !... 1.) if there is no CAPE, or
- !... 2.) cloud top is at model level just above LCL, or
- !... 3.) cloud top is within updraft source layer, or
- !... 4.) cloud-top detrainment layer begins within
- !... updraft source layer.
- !
- IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
- CLDHGT(LC)=0.
- DO NK=K,LTOP
- UMF(NK)=0.
- UDR(NK)=0.
- UER(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- ENDDO
- !
- ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
- ISHALL=0
- EXIT usl
- ELSE
- !
- !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
- ISHALL = 1
- IF(NU.EQ.NUCHM)THEN
- EXIT usl ! Shallow Convection from this layer
- ELSE
- ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
- DO NK=K,LTOP
- UMF(NK)=0.
- UDR(NK)=0.
- UER(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- ENDDO
- ENDIF
- ENDIF
- ENDIF trigger2
- END DO usl
- IF(ISHALL.EQ.1)THEN
- KSTART=MAX0(KPBL,KLCL)
- LET=KSTART
- endif
- !
- !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
- ! THIS LEVEL...
- !
- IF(LET.EQ.LTOP)THEN
- UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
- DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
- DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
- UER(LTOP)=0.
- UMF(LTOP)=0.
- ELSE
- !
- ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
- !
- DPTT=0.
- DO NJ=LET+1,LTOP
- DPTT=DPTT+DP(NJ)
- ENDDO
- DUMFDP=UMF(LET)/DPTT
- !
- !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
- ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
- !
- DO NK=LET+1,LTOP
- !
- !...entrainment is allowed at every level except for LTOP, so disallow
- !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
- !...so the the dilution factor due to entyrianment is not changed but
- !...the actual entrainment rate will change due due forced total
- !...detrainment in this layer...
- !
- IF(NK.EQ.LTOP)THEN
- UDR(NK) = UMF(NK-1)
- UER(NK) = 0.
- DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
- DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
- ELSE
- UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
- UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
- UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
- DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
- DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
- ENDIF
- IF(NK.GE.LET+2)THEN
- TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
- PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
- PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
- TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
- ENDIF
- ENDDO
- ENDIF
- !
- ! Initialize some arrays below cloud base and above cloud top...
- !
- DO NK=1,LTOP
- IF(T0(NK).GT.T00)ML=NK
- ENDDO
- DO NK=1,K
- IF(NK.GE.LC)THEN
- IF(NK.EQ.LC)THEN
- UMF(NK)=VMFLCL*DP(NK)/DPTHMX
- UER(NK)=VMFLCL*DP(NK)/DPTHMX
- ELSEIF(NK.LE.KPBL)THEN
- UER(NK)=VMFLCL*DP(NK)/DPTHMX
- UMF(NK)=UMF(NK-1)+UER(NK)
- ELSE
- UMF(NK)=VMFLCL
- UER(NK)=0.
- ENDIF
- TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
- QU(NK)=QMIX
- WU(NK)=WLCL
- ELSE
- TU(NK)=0.
- QU(NK)=0.
- UMF(NK)=0.
- WU(NK)=0.
- UER(NK)=0.
- ENDIF
- UDR(NK)=0.
- QDT(NK)=0.
- QLIQ(NK)=0.
- QICE(NK)=0.
- QLQOUT(NK)=0.
- QICOUT(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- RATIO2(NK)=0.
- CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
- EQFRC(NK)=1.0
- ENDDO
- !
- LTOP1=LTOP+1
- LTOPM1=LTOP-1
- !
- !...DEFINE VARIABLES ABOVE CLOUD TOP...
- !
- DO NK=LTOP1,KX
- UMF(NK)=0.
- UDR(NK)=0.
- UER(NK)=0.
- QDT(NK)=0.
- QLIQ(NK)=0.
- QICE(NK)=0.
- QLQOUT(NK)=0.
- QICOUT(NK)=0.
- DETLQ(NK)=0.
- DETIC(NK)=0.
- PPTLIQ(NK)=0.
- PPTICE(NK)=0.
- IF(NK.GT.LTOP1)THEN
- TU(NK)=0.
- QU(NK)=0.
- WU(NK)=0.
- ENDIF
- THTA0(NK)=0.
- THTAU(NK)=0.
- EMS(NK)=0.
- EMSD(NK)=0.
- TG(NK)=T0(NK)
- QG(NK)=Q0(NK)
- QLG(NK)=0.
- QIG(NK)=0.
- QRG(NK)=0.
- QSG(NK)=0.
- OMG(NK)=0.
- ENDDO
- OMG(KX+1)=0.
- DO NK=1,LTOP
- EMS(NK)=DP(NK)*DXSQ/G
- EMSD(NK)=1./EMS(NK)
- !
- !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
- !
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
- THTAU(NK)=TU(NK)*EXN(NK)
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
- THTA0(NK)=T0(NK)*EXN(NK)
- DDILFRC(NK) = 1./DILFRC(NK)
- OMG(NK)=0.
- ENDDO
- ! IF (XTIME.LT.10.)THEN
- ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
- ! * TMIX-T00,PMIX,QMIX,ABE
- ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
- ! * WLCL,CLDHGT
- ! ENDIF
- !
- !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
- !...AND MIDTROPOSPHERE IS USED.
- !
- WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
- WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
- WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
- VCONV=.5*(WSPD(KLCL)+WSPD(L5))
- !...for ETA model, DX is a function of location...
- ! TIMEC=DX(I,J)/VCONV
- TIMEC=DX/VCONV
- TADVEC=TIMEC
- TIMEC=AMAX1(1800.,TIMEC) ! 30 minutes >= TIMEC <= 60 minutes
- TIMEC=AMIN1(3600.,TIMEC)
- IF(ISHALL.EQ.1)TIMEC=2400. ! shallow convection TIMEC = 40 minutes
- NIC=NINT(TIMEC/DT)
- TIMEC=FLOAT(NIC)*DT
- !
- !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
- !
- IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
- SHSIGN=1.
- ELSE
- SHSIGN=-1.
- ENDIF
- VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
- (V0(LTOP)-V0(KLCL))
- VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
- PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
- PEF=AMAX1(PEF,.2)
- PEF=AMIN1(PEF,.9)
- !
- !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
- !
- CBH=(ZLCL-Z0(1))*3.281E-3
- IF(CBH.LT.3.)THEN
- RCBH=.02
- ELSE
- RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
- 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
- ENDIF
- IF(CBH.GT.25)RCBH=2.4
- PEFCBH=1./(1.+RCBH)
- PEFCBH=AMIN1(PEFCBH,.9)
- !
- !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
- !
- PEFF=.5*(PEF+PEFCBH)
- PEFF2 = PEFF ! JSK MODS
- IF(IPRNT)THEN
- ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
- WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
- CALL wrf_message( message )
- ! call flush(98)
- endif
- ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
- !*****************************************************************
- ! *
- ! COMPUTE DOWNDRAFT PROPERTIES *
- ! *
- !*****************************************************************
- !
- !
- TDER=0.
- devap:IF(ISHALL.EQ.1)THEN
- LFS = 1
- ELSE
- !
- !...start downdraft about 150 mb above cloud base...
- !
- ! KSTART=MAX0(KPBL,KLCL)
- ! KSTART=KPBL ! Changed 7/23/99
- KSTART=KPBL+1 ! Changed 7/23/99
- KLFS = LET-1
- DO NK = KSTART+1,KL
- DPPP = P0(KSTART)-P0(NK)
- ! IF(DPPP.GT.200.E2)THEN
- IF(DPPP.GT.150.E2)THEN
- KLFS = NK
- EXIT
- ENDIF
- ENDDO
- KLFS = MIN0(KLFS,LET-1)
- LFS = KLFS
- !
- !...if LFS is not at least 50 mb above cloud base (implying that the
- !...level of equil temp, LET, is just above cloud base) do not allow a
- !...downdraft...
- !
- IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
- THETED(LFS) = THETEE(LFS)
- QD(LFS) = Q0(LFS)
- !
- !...call tpmix2dd to find wet-bulb temp, qv...
- !
- call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
- THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
- !
- !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
- !
- TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
- RDD=P0(LFS)/(R*TVD(LFS))
- A1=(1.-PEFF)*AU0
- DMF(LFS)=-A1*RDD
- DER(LFS)=DMF(LFS)
- DDR(LFS)=0.
- RHBAR = RH(LFS)*DP(LFS)
- DPTT = DP(LFS)
- DO ND = LFS-1,KSTART,-1
- ND1 = ND+1
- DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
- DDR(ND)=0.
- DMF(ND)=DMF(ND1)+DER(ND)
- THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
- QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
- DPTT = DPTT+DP(ND)
- RHBAR = RHBAR+RH(ND)*DP(ND)
- ENDDO
- RHBAR = RHBAR/DPTT
- DMFFRC = 2.*(1.-RHBAR) ! Kain (2004) eq. 11
- DPDD = 0.
- !...Calculate melting effect
- !... first, compute total frozen precipitation generated...
- !
- pptmlt = 0.
- DO NK = KLCL,LTOP
- PPTMLT = PPTMLT+PPTICE(NK)
- ENDDO
- if(lc.lt.ml)then
- !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
- !...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
- !...12/14/98 jsk...
- DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
- else
- DTMELT = 0.
- endif
- LDT = MIN0(LFS-1,KSTART-1)
- !
- call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
- !
- tz(kstart) = tz(kstart)-dtmelt
- ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
- QSS=0.622*ES/(P0(KSTART)-ES)
- THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))* &
- EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
- !....
- LDT = MIN0(LFS-1,KSTART-1)
- DO ND = LDT,1,-1
- DPDD = DPDD+DP(ND)
- THETED(ND) = THETED(KSTART)
- QD(ND) = QD(KSTART)
- !
- !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
- !
- call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
- qsd(nd) = qss
- !
- !...specify RH decrease of 20%/km in downdraft...
- !
- RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
- !
- !...adjust downdraft TEMP, Q to specified RH:
- !
- IF(RHH.LT.1.)THEN
- DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
- RL=XLV0-XLV1*TZ(ND)
- DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
- T1RH=TZ(ND)+DTMP
- ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
- QSRH=0.622*ES/(P0(ND)-ES)
- !
- !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
- !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
- !
- IF(QSRH.LT.QD(ND))THEN
- QSRH=QD(ND)
- T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
- ENDIF
- TZ(ND)=T1RH
- QSS=QSRH
- QSD(ND) = QSS
- ENDIF
- TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
- IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
- LDB=ND
- EXIT
- ENDIF
- ENDDO
- IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
- DO ND=LDT,LDB,-1
- ND1 = ND+1
- DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
- DER(ND) = 0.
- DMF(ND) = DMF(ND1)+DDR(ND)
- TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
- QD(ND)=QSD(nd)
- THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
- ENDDO
- ENDIF
- ENDIF
- ENDIF devap
- !
- !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
- !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
- !
- d_mf: IF(TDER.LT.1.)THEN
- ! WRITE(98,3004)I,J
- !3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
- PPTFLX=TRPPT
- CPR=TRPPT
- TDER=0.
- CNDTNF=0.
- UPDINC=1.
- LDB=LFS
- DO NDK=1,LTOP
- DMF(NDK)=0.
- DER(NDK)=0.
- DDR(NDK)=0.
- THTAD(NDK)=0.
- WD(NDK)=0.
- TZ(NDK)=0.
- QD(NDK)=0.
- ENDDO
- AINCM2=100.
- ELSE
- DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
- UPDINC=1.
- IF(TDER*DDINC.GT.TRPPT)THEN
- DDINC = TRPPT/TDER
- ENDIF
- TDER = TDER*DDINC
- DO NK=LDB,LFS
- DMF(NK)=DMF(NK)*DDINC
- DER(NK)=DER(NK)*DDINC
- DDR(NK)=DDR(NK)*DDINC
- ENDDO
- CPR=TRPPT
- PPTFLX = TRPPT-TDER
- PEFF=PPTFLX/TRPPT
- IF(IPRNT)THEN
- ! write(98,*)'PRECIP EFFICIENCY =',PEFF
- write(message,*)'PRECIP EFFICIENCY =',PEFF
- CALL wrf_message(message)
- ! call flush(98)
- ENDIF
- !
- !
- !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
- ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
- ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
- !
- ! DO NK=LC,LFS
- ! UMF(NK)=UMF(NK)*UPDINC
- ! UDR(NK)=UDR(NK)*UPDINC
- ! UER(NK)=UER(NK)*UPDINC
- ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
- ! PPTICE(NK)=PPTICE(NK)*UPDINC
- ! DETLQ(NK)=DETLQ(NK)*UPDINC
- ! DETIC(NK)=DETIC(NK)*UPDINC
- ! ENDDO
- !
- !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
- !...DOWNDRAFT...
- !
- IF(LDB.GT.1)THEN
- DO NK=1,LDB-1
- DMF(NK)=0.
- DER(NK)=0.
- DDR(NK)=0.
- WD(NK)=0.
- TZ(NK)=0.
- QD(NK)=0.
- THTAD(NK)=0.
- ENDDO
- ENDIF
- DO NK=LFS+1,KX
- DMF(NK)=0.
- DER(NK)=0.
- DDR(NK)=0.
- WD(NK)=0.
- TZ(NK)=0.
- QD(NK)=0.
- THTAD(NK)=0.
- ENDDO
- DO NK=LDT+1,LFS-1
- TZ(NK)=0.
- QD(NK)=0.
- THTAD(NK)=0.
- ENDDO
- ENDIF d_mf
- !
- !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
- ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
- ! IN THAT LAYER INITIALLY...
- !
- AINCMX=1000.
- LMAX=MAX0(KLCL,LFS)
- DO NK=LC,LMAX
- IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
- AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
- AINCMX=AMIN1(AINCMX,AINCM1)
- ENDIF
- ENDDO
- AINC=1.
- IF(AINCMX.LT.AINC)AINC=AINCMX
- !
- !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
- !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
- !...CLOSURE...
- !
- TDER2=TDER
- PPTFL2=PPTFLX
- DO NK=1,LTOP
- DETLQ2(NK)=DETLQ(NK)
- DETIC2(NK)=DETIC(NK)
- UDR2(NK)=UDR(NK)
- UER2(NK)=UER(NK)
- DDR2(NK)=DDR(NK)
- DER2(NK)=DER(NK)
- UMF2(NK)=UMF(NK)
- DMF2(NK)=DMF(NK)
- ENDDO
- FABE=1.
- STAB=0.95
- NOITR=0
- ISTOP=0
- !
- IF(ISHALL.EQ.1)THEN ! First for shallow convection
- !
- ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
- ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
- ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
- !
- !...find the maximum TKE value between LC and KLCL...
- ! TKEMAX = 0.
- TKEMAX = 5.
- ! DO 173 K = LC,KLCL
- ! NK = KX-K+1
- ! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
- ! 173 CONTINUE
- ! TKEMAX = AMIN1(TKEMAX,10.)
- ! TKEMAX = AMAX1(TKEMAX,5.)
- !c TKEMAX = 10.
- !c...3_24_99...DPMIN was changed for shallow convection so that it is the
- !c... the same as for deep convection (5.E3). Since this doubles
- !c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
- !c... lation of EVAC...
- !c EVAC = TKEMAX*0.1
- EVAC = 0.5*TKEMAX*0.1
- ! AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
- ! AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
- AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
- TDER=TDER2*AINC
- PPTFLX=PPTFL2*AINC
- DO NK=1,LTOP
- UMF(NK)=UMF2(NK)*AINC
- DMF(NK)=DMF2(NK)*AINC
- DETLQ(NK)=DETLQ2(NK)*AINC
- DETIC(NK)=DETIC2(NK)*AINC
- UDR(NK)=UDR2(NK)*AINC
- UER(NK)=UER2(NK)*AINC
- DER(NK)=DER2(NK)*AINC
- DDR(NK)=DDR2(NK)*AINC
- ENDDO
- ENDIF ! Otherwise for deep convection
- ! use iterative procedure to find mass fluxes...
- iter: DO NCOUNT=1,10
- !
- !*****************************************************************
- ! *
- ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
- ! *
- !*****************************************************************
- !
- !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
- !...SATISFY MASS CONTINUITY...
- !
- DTT=TIMEC
- DO NK=1,LTOP
- DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
- IF(NK.GT.1)THEN
- OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
- ABSOMG = ABS(OMG(NK))
- ABSOMGTC = ABSOMG*TIMEC
- FRDP = 0.75*DP(NK-1)
- IF(ABSOMGTC.GT.FRDP)THEN
- DTT1 = FRDP/ABSOMG
- DTT=AMIN1(DTT,DTT1)
- ENDIF
- ENDIF
- ENDDO
- DO NK=1,LTOP
- THPA(NK)=THTA0(NK)
- QPA(NK)=Q0(NK)
- NSTEP=NINT(TIMEC/DTT+1)
- DTIME=TIMEC/FLOAT(NSTEP)
- FXM(NK)=OMG(NK)*DXSQ/G
- ENDDO
- !
- !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
- !
- DO NTC=1,NSTEP
- !
- !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
- !...SIGN OF OMEGA...
- !
- DO NK=1,LTOP
- THFXIN(NK)=0.
- THFXOUT(NK)=0.
- QFXIN(NK)=0.
- QFXOUT(NK)=0.
- ENDDO
- DO NK=2,LTOP
- IF(OMG(NK).LE.0.)THEN
- THFXIN(NK)=-FXM(NK)*THPA(NK-1)
- QFXIN(NK)=-FXM(NK)*QPA(NK-1)
- THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
- QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
- ELSE
- THFXOUT(NK)=FXM(NK)*THPA(NK)
- QFXOUT(NK)=FXM(NK)*QPA(NK)
- THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
- QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
- ENDIF
- ENDDO
- !
- !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
- !
- DO NK=1,LTOP
- THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
- THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
- DTIME*EMSD(NK)
- QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)- &
- QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
- ENDDO
- ENDDO
- DO NK=1,LTOP
- THTAG(NK)=THPA(NK)
- QG(NK)=QPA(NK)
- ENDDO
- !
- !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
- !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
- !
- DO NK=1,LTOP
- IF(QG(NK).LT.0.)THEN
- IF(NK.EQ.1)THEN ! JSK MODS
- ! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
- ! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
- CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
- ENDIF ! JSK MODS
- NK1=NK+1
- IF(NK.EQ.LTOP)THEN
- NK1=KLCL
- ENDIF
- TMA=QG(NK1)*EMS(NK1)
- TMB=QG(NK-1)*EMS(NK-1)
- TMM=(QG(NK)-1.E-9)*EMS(NK )
- BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
- ACOEFF=BCOEFF*TMA/TMB
- TMB=TMB*(1.-BCOEFF)
- TMA=TMA*(1.-ACOEFF)
- IF(NK.EQ.LTOP)THEN
- QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
- ! IF(ABS(QVDIFF).GT.1.)THEN
- ! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
- ! QVDIFF, &
- ! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
- ! 'VALUES IN KAIN-FRITSCH'
- ! ENDIF
- ENDIF
- QG(NK)=1.E-9
- QG(NK1)=TMA*EMSD(NK1)
- QG(NK-1)=TMB*EMSD(NK-1)
- ENDIF
- ENDDO
- TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
- IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
- ! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
- ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
- ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
- ISTOP=1
- IPRNT=.TRUE.
- EXIT iter
- ENDIF
- !
- !...CONVERT THETA TO T...
- !
- DO NK=1,LTOP
- EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
- TG(NK)=THTAG(NK)/EXN(NK)
- TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
- ENDDO
- IF(ISHALL.EQ.1)THEN
- EXIT iter
- ENDIF
- !
- !*******************************************************************
- ! *
- ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
- ! *
- !*******************************************************************
- !
- !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
- !
- ! THMIX=0.
- TMIX=0.
- QMIX=0.
- !
- !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
- !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
- !...LAYERS...
- !
- DO NK=LC,KPBL
- TMIX=TMIX+DP(NK)*TG(NK)
- QMIX=QMIX+DP(NK)*QG(NK)
- ENDDO
- TMIX=TMIX/DPTHMX
- QMIX=QMIX/DPTHMX
- ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
- QSS=0.622*ES/(PMIX-ES)
- !
- !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
- !
- IF(QMIX.GT.QSS)THEN
- RL=XLV0-XLV1*TMIX
- CPM=CP*(1.+0.887*QMIX)
- DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
- DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
- TMIX=TMIX+RL/CP*DQ
- QMIX=QMIX-DQ
- TLCL=TMIX
- ELSE
- QMIX=AMAX1(QMIX,0.)
- EMIX=QMIX*PMIX/(0.622+QMIX)
- astrt=1.e-3
- binc=0.075
- a1=emix/aliq
- tp=(a1-astrt)/binc
- indlu=int(tp)+1
- value=(indlu-1)*binc+astrt
- aintrp=(a1-value)/binc
- tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
- TLCL=AMIN1(TLCL,TMIX)
- ENDIF
- TVLCL=TLCL*(1.+0.608*QMIX)
- ZLCL = ZMIX+(TLCL-TMIX)/GDRY
- DO NK = LC,KL
- KLCL=NK
- IF(ZLCL.LE.Z0(NK))THEN
- EXIT
- ENDIF
- ENDDO
- K=KLCL-1
- DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
- !
- !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
- !
- TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
- QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
- TVEN=TENV*(1.+0.608*QENV)
- PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
- THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
- EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
- !
- !...COMPUTE ADJUSTED ABE(ABEG).
- !
- ABEG=0.
- DO NK=K,LTOPM1
- NK1=NK+1
- THETEU(NK1) = THETEU(NK)
- !
- call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
- !
- TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
- IF(NK.EQ.K)THEN
- DZZ=Z0(KLCL)-ZLCL
- DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
- ELSE
- DZZ=DZA(NK)
- DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
- ENDIF
- IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
- !
- !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
- !
- CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
- THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
- ENDDO
- !
- !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
- !...THE PERIOD TIMEC...
- !
- IF(NOITR.EQ.1)THEN
- ! write(98,*)' '
- ! write(98,*)'TAU, I, J, =',NTSD,I,J
- ! WRITE(98,1060)FABE
- ! GOTO 265
- EXIT iter
- ENDIF
- DABE=AMAX1(ABE-ABEG,0.1*ABE)
- FABE=ABEG/ABE
- IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
- ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
- ! *GRID POINT; NO CONVECTION ALLOWED!'
- RETURN
- ENDIF
- IF(NCOUNT.NE.1)THEN
- IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
- NOITR=1
- AINC=AINCOLD
- CYCLE iter
- ENDIF
- DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
- IF(DFDA.GT.0.)THEN
- NOITR=1
- AINC=AINCOLD
- CYCLE iter
- ENDIF
- ENDIF
- AINCOLD=AINC
- FABEOLD=FABE
- IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
- ! write(98,*)' '
- ! write(98,*)'TAU, I, J, =',NTSD,I,J
- ! WRITE(98,1055)FABE
- ! GOTO 265
- EXIT
- ENDIF
- IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
- EXIT iter
- ELSE
- IF(NCOUNT.GT.10)THEN
- ! write(98,*)' '
- ! write(98,*)'TAU, I, J, =',NTSD,I,J
- ! WRITE(98,1060)FABE
- ! GOTO 265
- EXIT
- ENDIF
- !
- !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
- !...MASS FLUX BY THE FACTOR AINC:
- !
- IF(FABE.EQ.0.)THEN
- AINC=AINC*0.5
- ELSE
- IF(DABE.LT.1.e-4)THEN
- NOITR=1
- AINC=AINCOLD
- CYCLE iter
- ELSE
- AINC=AINC*STAB*ABE/DABE
- ENDIF
- ENDIF
- ! AINC=AMIN1(AINCMX,AINC)
- AINC=AMIN1(AINCMX,AINC)
- !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
- !...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
- IF(AINC.LT.0.05)then
- RETURN ! JSK MODS
- ENDIF
- ! AINC=AMAX1(AINC,0.05) ! JSK MODS
- TDER=TDER2*AINC
- PPTFLX=PPTFL2*AINC
- ! IF (XTIME.LT.10.)THEN
- ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
- ! * FABEOLD,AINCOLD
- ! ENDIF
- DO NK=1,LTOP
- UMF(NK)=UMF2(NK)*AINC
- DMF(NK)=DMF2(NK)*AINC
- DETLQ(NK)=DETLQ2(NK)*AINC
- DETIC(NK)=DETIC2(NK)*AINC
- UDR(NK)=UDR2(NK)*AINC
- UER(NK)=UER2(NK)*AINC
- DER(NK)=DER2(NK)*AINC
- DDR(NK)=DDR2(NK)*AINC
- ENDDO
- !
- !...GO BACK UP FOR ANOTHER ITERATION...
- !
- ENDIF
- ENDDO iter
- !
- !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
- !
- !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
- !...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
- !
- ! Redistribute hydormeteors according to the final mass-flux values:
- !
- IF(CPR.GT.0.)THEN
- FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
- ELSE
- FRC2=0.
- ENDIF
- DO NK=1,LTOP
- QLPA(NK)=QL0(NK)
- QIPA(NK)=QI0(NK)
- QRPA(NK)=QR0(NK)
- QSPA(NK)=QS0(NK)
- RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
- SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
- ENDDO
- DO NTC=1,NSTEP
- !
- !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
- !...BASED ON THE SIGN OF OMEGA...
- !
- DO NK=1,LTOP
- QLFXIN(NK)=0.
- QLFXOUT(NK)=0.
- QIFXIN(NK)=0.
- QIFXOUT(NK)=0.
- QRFXIN(NK)=0.
- QRFXOUT(NK)=0.
- QSFXIN(NK)=0.
- QSFXOUT(NK)=0.
- ENDDO
- DO NK=2,LTOP
- IF(OMG(NK).LE.0.)THEN
- QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
- QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
- QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
- QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
- QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
- QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
- QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
- QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
- ELSE
- QLFXOUT(NK)=FXM(NK)*QLPA(NK)
- QIFXOUT(NK)=FXM(NK)*QIPA(NK)
- QRFXOUT(NK)=FXM(NK)*QRPA(NK)
- QSFXOUT(NK)=FXM(NK)*QSPA(NK)
- QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
- QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
- QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
- QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
- ENDIF
- ENDDO
- !
- !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
- !
- DO NK=1,LTOP
- QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
- QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
- QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
- QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
- ENDDO
- ENDDO
- DO NK=1,LTOP
- QLG(NK)=QLPA(NK)
- QIG(NK)=QIPA(NK)
- QRG(NK)=QRPA(NK)
- QSG(NK)=QSPA(NK)
- ENDDO
- !
- !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
- !...GRID POINT...
- !
- ! IF (XTIME.LT.10.)THEN
- ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
- ! ENDIF
- IF(IPRNT)THEN
- ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
- WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
- CALL wrf_message(message)
- ! call flush(98)
- endif
- !
- !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
- !
- !297 IF(IPRNT)then
- IF(IPRNT)then
- ! if(I.eq.16 .and. J.eq.41)then
- ! IF(ISTOP.EQ.1)THEN
- write(98,*)
- ! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
- write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
- TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
- call wrf_message(message)
- write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
- DTRH,TENV
- call wrf_message(message)
- WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
- TMIX-T00,PMIX,QMIX,ABE
- call wrf_message(message)
- WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
- WLCL,CLDHGT(LC)
- call wrf_message(message)
- WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
- call wrf_message(message)
- write(message,*)'PRECIP EFFICIENCY =',PEFF
- call wrf_message(message)
- WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
- call wrf_message(message)
- ! ENDIF
- !!!!! HERE !!!!!!!
- WRITE(message,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
- ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
- ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
- call wrf_message(message)
- write(message,*)'just before DO 300...'
- call wrf_message(message)
- ! call flush(98)
- DO NK=1,LTOP
- K=LTOP-NK+1
- DTT=(TG(K)-T0(K))*86400./TIMEC
- RL=XLV0-XLV1*TG(K)
- DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
- UDFRC=UDR(K)*TIMEC*EMSD(K)
- UEFRC=UER(K)*TIMEC*EMSD(K)
- DDFRC=DDR(K)*TIMEC*EMSD(K)
- DEFRC=-DER(K)*TIMEC*EMSD(K)
- WRITE(message,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
- UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
- W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
- TIMEC*EMSD(K)*1.E3
- call wrf_message(message)
- ENDDO
- WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
- 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
- call wrf_message(message)
- DO NK=1,KL
- K=KX-NK+1
- DTT=TG(K)-T0(K)
- TUC=TU(K)-T00
- IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
- TDC=TZ(K)-T00
- IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
- IF(T0(K).LT.T00)THEN
- ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
- ELSE
- ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
- ENDIF
- QGS=ES*0.622/(P0(K)-ES)
- RH0=Q0(K)/QES(K)
- RHG=QG(K)/QGS
- WRITE(message,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
- TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
- 1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
- QSG(K)*1000.,RH0,RHG
- call wrf_message(message)
- ENDDO
- !
- !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
- !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
- !
- ! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
- ! IF(ISHALL.NE.1)THEN
- ! write(98,4421)i,j,iyr,imo,idy,ihr,imn
- ! write(98)i,j,iyr,imo,idy,ihr,imn,kl
- ! 4421 format(7i4)
- ! write(98,4422)kl
- ! 4422 format(i6)
- DO 310 NK = 1,KL
- k = kl - nk + 1
- write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
- u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
- ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
- ! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
- ! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
- 310 CONTINUE
- IF(ISTOP.EQ.1)THEN
- CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
- ENDIF
- ! ENDIF
- 4455 format(8f11.3)
- ENDIF
- CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
- PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
- RAINCV(I,J)=DT*PRATEC(I,J) ! PPT FB MODS
- ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
- ! RNC=0.1*TIMEC*PPTFLX/DXSQ
- RNC=RAINCV(I,J)*NIC
- IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
- ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
- !
- ! EVALUATE MOISTURE BUDGET...
- !
- QINIT=0.
- QFNL=0.
- DPT=0.
- DO 315 NK=1,LTOP
- DPT=DPT+DP(NK)
- QINIT=QINIT+Q0(NK)*EMS(NK)
- QFNL=QFNL+QG(NK)*EMS(NK)
- QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
- 315 CONTINUE
- QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) ! PPT FB MODS
- ! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
- ERR2=(QFNL-QINIT)*100./QINIT
- IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
- IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
- ! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
- ! WRITE(99,1110)QINIT,QFNL,ERR2
- IPRNT=.TRUE.
- ISTOP=1
- write(98,4422)kl
- 4422 format(i6)
- DO 311 NK = 1,KL
- k = kl - nk + 1
- ! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
- ! u0(k),v0(k),W0AVG1D(K),dp(k)
- ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
- ! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
- ! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
- WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
- U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
- 311 CONTINUE
- ! call flush(98)
- ! GOTO 297
- ! STOP 'QVERR'
- ENDIF
- 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
- 4456 format(8f12.3)
- IF(PPTFLX.GT.0.)THEN
- RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
- ELSE
- RELERR=0.
- ENDIF
- IF(IPRNT)THEN
- WRITE(98,1120)RELERR
- WRITE(98,*)'TDER, CPR, TRPPT =', &
- TDER,CPR*AINC,TRPPT*AINC
- ENDIF
- !
- !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
- !
- !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
- !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
- !
- IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
- NCA(I,J) = REAL(NIC)*DT
- IF(ISHALL.EQ.1)THEN
- TIMEC = 2400.
- NCA(I,J) = CUDT*60.
- NSHALL = NSHALL+1
- ENDIF
- DO K=1,KX
- ! IF(IMOIST(INEST).NE.2)THEN
- !
- !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
- !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
- !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
- !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
- !...OF QG...
- !
- ! RLC=XLV0-XLV1*TG(K)
- ! RLS=XLS0-XLS1*TG(K)
- ! CPM=CP*(1.+0.887*QG(K))
- ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
- ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
- ! DQLDT(I,J,NK)=0.
- ! DQIDT(I,J,NK)=0.
- ! DQRDT(I,J,NK)=0.
- ! DQSDT(I,J,NK)=0.
- ! ELSE
- !
- !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
- !
- IF(.NOT. F_QI .and. warm_rain)THEN
- CPM=CP*(1.+0.887*QG(K))
- TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
- DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
- DQIDT(K)=0.
- DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
- DQSDT(K)=0.
- ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
- !
- !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
- !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
- !
- CPM=CP*(1.+0.887*QG(K))
- IF(K.LE.ML)THEN
- TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
- ELSEIF(K.GT.ML)THEN
- TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
- ENDIF
- DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
- DQIDT(K)=0.
- DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
- DQSDT(K)=0.
- ELSEIF(F_QI) THEN
- !
- !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
- !...OF HYDROMETEORS DIRECTLY...
- !
- DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
- DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
- DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
- IF (F_QS) THEN
- DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
- ELSE
- DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
- ENDIF
- ELSE
- ! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
- CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
- ENDIF
- DTDT(K)=(TG(K)-T0(K))/TIMEC
- DQDT(K)=(QG(K)-Q0(K))/TIMEC
- ENDDO
- PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
- RAINCV(I,J)=DT*PRATEC(I,J)
- ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
- ! RNC=0.1*TIMEC*PPTFLX/DXSQ
- RNC=RAINCV(I,J)*NIC
- 909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
- ! write (98,909)I,J,RNC
- ! write (6,909)I,J,RNC
- ! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
- ! * NCCNT
- ! call flush(98)
- 1000 FORMAT(' ',10A8)
- 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
- 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
- 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
- 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
- ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
- I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
- ' CAPE=',0PF7.1)
- 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
- E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
- F8.1)
- 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
- ,F6.3,'VWS=',F5.2)
- !1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
- ! ', NO MORE MASS FLUX IS ALLOWED!')
- !1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
- ! &DEGREE OF STABILIZATION! FABE= ',F6.4)
- 1070 FORMAT (16A8)
- 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
- 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
- 2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
- 1085 FORMAT (A3,16A7,2A8)
- 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
- 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
- 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
- E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
- 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
- ' TOTAL WATER CHANGE =',F8.2,'%')
- ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
- 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
- !
- !-----------------------------------------------------------------------
- !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
- !-----------------------------------------------------------------------
- !
- CUTOP(I,J)=REAL(LTOP)
- CUBOT(I,J)=REAL(LCL)
- !
- !-----------------------------------------------------------------------
- END SUBROUTINE KF_eta_PARA
- !********************************************************************
- ! ***********************************************************************
- SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
- !
- ! Lookup table variables:
- ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
- ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
- ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
- ! REAL, SAVE, DIMENSION(1:200) :: ALU
- ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
- ! End of Lookup table variables:
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
- REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
- REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
- REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
- TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
- INTEGER :: IPTB,ITHTB
- !-----------------------------------------------------------------------
- !c******** LOOKUP TABLE VARIABLES... ****************************
- ! parameter(kfnt=250,kfnp=220)
- !c
- ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
- ! * alu(200),rdpr,rdthk,plutop
- !C***************************************************************
- !c
- !c***********************************************************************
- !c scaling pressure and tt table index
- !c***********************************************************************
- !c
- tp=(p-plutop)*rdpr
- qq=tp-aint(tp)
- iptb=int(tp)+1
- !
- !***********************************************************************
- ! base and scaling factor for the
- !***********************************************************************
- !
- ! scaling the and tt table index
- bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
- tth=(thes-bth)*rdthk
- pp =tth-aint(tth)
- ithtb=int(tth)+1
- IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
- write(98,*)'**** OUT OF BOUNDS *********'
- ! call flush(98)
- ENDIF
- !
- t00=ttab(ithtb ,iptb )
- t10=ttab(ithtb+1,iptb )
- t01=ttab(ithtb ,iptb+1)
- t11=ttab(ithtb+1,iptb+1)
- !
- q00=qstab(ithtb ,iptb )
- q10=qstab(ithtb+1,iptb )
- q01=qstab(ithtb ,iptb+1)
- q11=qstab(ithtb+1,iptb+1)
- !
- !***********************************************************************
- ! parcel temperature
- !***********************************************************************
- !
- temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
- !
- qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
- !
- DQ=QS-QU
- IF(DQ.LE.0.)THEN
- QNEW=QU-QS
- QU=QS
- ELSE
- !
- ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
- ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
- !
- QNEW=0.
- QTOT=QLIQ+QICE
- !
- ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
- ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
- ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
- ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
- ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
- !
- !...subsaturated values only occur in calculations involving various mixtures of
- !...updraft and environmental air for estimation of entrainment and detrainment.
- !...For these purposes, assume that reasonable estimates can be given using
- !...liquid water saturation calculations only - i.e., ignore the effect of the
- !...ice phase in this process only...will not affect conservative properties...
- !
- IF(QTOT.GE.DQ)THEN
- qliq=qliq-dq*qliq/(qtot+1.e-10)
- qice=qice-dq*qice/(qtot+1.e-10)
- QU=QS
- ELSE
- RLL=XLV0-XLV1*TEMP
- CPP=1004.5*(1.+0.89*QU)
- IF(QTOT.LT.1.E-10)THEN
- !
- !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
- TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
- ELSE
- !
- !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
- ! THE TEMPERATURE IS GIVEN BY:
- !
- TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
- QU=QU+QTOT
- QTOT=0.
- QLIQ=0.
- QICE=0.
- ENDIF
- ENDIF
- ENDIF
- TU=TEMP
- qnewlq=qnew
- qnewic=0.
- !
- END SUBROUTINE TPMIX2
- !******************************************************************************
- SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
- REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
- REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
- !-----------------------------------------------------------------------
- !
- !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
- !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
- !...TTFRZ TO TBFRZ...
- !...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
- !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
- !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
- !
- RLC=2.5E6-2369.276*(TU-273.16)
- RLS=2833922.-259.532*(TU-273.16)
- RLF=RLS-RLC
- CPP=1004.5*(1.+0.89*QU)
- !
- ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
- ! FOR SATURATION VAPOR PRESSURE...
- !
- A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
- DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
- TU = TU+DTFRZ
-
- ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
- QS = ES*0.622/(P-ES)
- !
- !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
- !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
- !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
- !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
- !...TEMPERATURE TO THE SATURATION VALUE...
- !
- DQEVAP = QS-QU
- QICE = QICE-DQEVAP
- QU = QU+DQEVAP
- PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
- THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
- !
- END SUBROUTINE DTFRZNEW
- ! --------------------------------------------------------------------------------
- SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
- QNEWIC,QLQOUT,QICOUT,G)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
- ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
- ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
- ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
- ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
- REAL, INTENT(IN ) :: G
- REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
- REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
- REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
- !
- QTOT=QLIQ+QICE
- QNEW=QNEWLQ+QNEWIC
- !
- ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
- ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
- ! LEVELS...
- !
- QEST=0.5*(QTOT+QNEW)
- G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
- IF(G1.LT.0.0)G1=0.
- WAVG=0.5*(SQRT(WTW)+SQRT(G1))
- CONV=RATE*DZ/WAVG ! KF90 Eq. 9
- !
- ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
- ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
- ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
- ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
- !
- RATIO3=QNEWLQ/(QNEW+1.E-8)
- ! OLDQ=QTOT
- QTOT=QTOT+0.6*QNEW
- OLDQ=QTOT
- RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
- QTOT=QTOT*EXP(-CONV) ! KF90 Eq. 9
- !
- ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
- ! PARCEL AT THIS LEVEL...
- !
- DQ=OLDQ-QTOT
- QLQOUT=RATIO4*DQ
- QICOUT=(1.-RATIO4)*DQ
- !
- ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
- ! LATE VERTICAL VELOCITY
- !
- PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
- WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
- IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
- !
- ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
- ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
- !
- QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
- QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
- QNEWLQ=0.
- QNEWIC=0.
- END SUBROUTINE CONDLOAD
- ! ----------------------------------------------------------------------
- SUBROUTINE PROF5(EQ,EE,UD)
- !
- !***********************************************************************
- !***** GAUSSIAN TYPE MIXING PROFILE....******************************
- !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
- ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
- ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
- ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
- ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
- ! JACK KAIN
- ! 7/6/89
- ! Solves for KF90 Eq. 2
- !
- !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: EQ
- REAL, INTENT(INOUT) :: EE,UD
- REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
- DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
- 0.9372980,0.33267,0.166666667,0.202765151/
- X=(EQ-0.5)/SIGMA
- Y=6.*EQ-3.
- EY=EXP(Y*Y/(-2))
- E45=EXP(-4.5)
- T2=1./(1.+P*ABS(Y))
- T1=0.500498
- C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
- C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
- IF(Y.GE.0.)THEN
- EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
- UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
- EQ)
- ELSE
- EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
- UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
- EQ/2.-EQ)
- ENDIF
- EE=EE/FE
- UD=UD/FE
- END SUBROUTINE PROF5
- ! ------------------------------------------------------------------------
- SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
- !
- ! Lookup table variables:
- ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
- ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
- ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
- ! REAL, SAVE, DIMENSION(1:200) :: ALU
- ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
- ! End of Lookup table variables:
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: P,THES
- REAL, INTENT(INOUT) :: TS,QS
- INTEGER, INTENT(IN ) :: i,j ! avail for debugging
- REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
- INTEGER :: IPTB,ITHTB
- CHARACTER*256 :: MESS
- !-----------------------------------------------------------------------
- !
- !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
- ! parameter(kfnt=250,kfnp=220)
- !
- ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
- ! alu(200),rdpr,rdthk,plutop
- !***************************************************************
- !
- !***********************************************************************
- ! scaling pressure and tt table index
- !***********************************************************************
- !
- tp=(p-plutop)*rdpr
- qq=tp-aint(tp)
- iptb=int(tp)+1
- !
- !***********************************************************************
- ! base and scaling factor for the
- !***********************************************************************
- !
- ! scaling the and tt table index
- bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
- tth=(thes-bth)*rdthk
- pp =tth-aint(tth)
- ithtb=int(tth)+1
- !
- t00=ttab(ithtb ,iptb )
- t10=ttab(ithtb+1,iptb )
- t01=ttab(ithtb ,iptb+1)
- t11=ttab(ithtb+1,iptb+1)
- !
- q00=qstab(ithtb ,iptb )
- q10=qstab(ithtb+1,iptb )
- q01=qstab(ithtb ,iptb+1)
- q11=qstab(ithtb+1,iptb+1)
- !
- !***********************************************************************
- ! parcel temperature and saturation mixing ratio
- !***********************************************************************
- !
- ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
- !
- qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
- !
- END SUBROUTINE TPMIX2DD
- ! -----------------------------------------------------------------------
- SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
- !
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
- REAL, INTENT(INOUT) :: THT1
- REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
- T00,P00,C1,C2,C3,C4,C5
- INTEGER :: INDLU
- !-----------------------------------------------------------------------
- DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
- 0.278296,1.0723E-3/
- !
- ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
- !
- ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
- ! For example, KF90 Eq. 10 no longer used
- !
- EE=Q1*P1/(0.622+Q1)
- ! TLOG=ALOG(EE/ALIQ)
- ! ...calculate LOG term using lookup table...
- !
- astrt=1.e-3
- ainc=0.075
- a1=ee/aliq
- tp=(a1-astrt)/ainc
- indlu=int(tp)+1
- value=(indlu-1)*ainc+astrt
- aintrp=(a1-value)/ainc
- tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
- !
- TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
- TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
- THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
- THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
- !
- END SUBROUTINE ENVIRTHT
- ! ***********************************************************************
- !====================================================================
- SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
- RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
- SVP1,SVP2,SVP3,SVPT0, &
- P_FIRST_SCALAR,restart,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_QI,P_QS,P_FIRST_SCALAR
- REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
- RTHCUTEN, &
- RQVCUTEN, &
- RQCCUTEN, &
- RQRCUTEN, &
- RQICUTEN, &
- RQSCUTEN
- REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
- REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
- INTEGER :: i, j, k, itf, jtf, ktf
- REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
- 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.
- RQCCUTEN(i,k,j)=0.
- RQRCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- 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
- IF (P_QS .ge. P_FIRST_SCALAR) THEN
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- RQSCUTEN(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- DO j=jts,jtf
- DO i=its,itf
- NCA(i,j)=-100.
- ENDDO
- ENDDO
- DO j=jts,jtf
- DO k=kts,ktf
- DO i=its,itf
- W0AVG(i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- endif
-
- CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
- END SUBROUTINE kf_eta_init
- !-------------------------------------------------------
- subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
- !
- ! This subroutine is a lookup table.
- ! Given a series of series of saturation equivalent potential
- ! temperatures, the temperature is calculated.
- !
- !--------------------------------------------------------------------
- IMPLICIT NONE
- !--------------------------------------------------------------------
- ! Lookup table variables
- ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
- ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
- ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
- ! REAL, SAVE, DIMENSION(1:200) :: ALU
- ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
- ! End of Lookup table variables
- INTEGER :: KP,IT,ITCNT,I
- REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
- TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
- ASTRT,AINC,A1,THTGS
- ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
- REAL :: ALIQ,BLIQ,CLIQ,DLIQ
- REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
- !
- ! equivalent potential temperature increment
- data dth/1./
- ! minimum starting temp
- data tmin/150./
- ! tolerance for accuracy of temperature
- data toler/0.001/
- ! top pressure (pascals)
- plutop=5000.0
- ! bottom pressure (pascals)
- pbot=110000.0
- ALIQ = SVP1*1000.
- BLIQ = SVP2
- CLIQ = SVP2*SVPT0
- DLIQ = SVP3
- !
- ! compute parameters
- !
- ! 1._over_(sat. equiv. theta increment)
- rdthk=1./dth
- ! pressure increment
- !
- DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
- ! dpr=(pbot-plutop)/REAL(kfnp-1)
- ! 1._over_(pressure increment)
- rdpr=1./dpr
- ! compute the spread of thes
- ! thespd=dth*(kfnt-1)
- !
- ! calculate the starting sat. equiv. theta
- !
- temp=tmin
- p=plutop-dpr
- do kp=1,kfnp
- p=p+dpr
- es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
- qs=0.622*es/(p-es)
- pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
- the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
- (1.+0.81*qs))
- enddo
- !
- ! compute temperatures for each sat. equiv. potential temp.
- !
- p=plutop-dpr
- do kp=1,kfnp
- thes=the0k(kp)-dth
- p=p+dpr
- do it=1,kfnt
- ! define sat. equiv. pot. temp.
- thes=thes+dth
- ! iterate to find temperature
- ! find initial guess
- if(it.eq.1) then
- tgues=tmin
- else
- tgues=ttab(it-1,kp)
- endif
- es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
- qs=0.622*es/(p-es)
- pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
- thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
- (1.+0.81*qs))
- f0=thgues-thes
- t1=tgues-0.5*f0
- t0=tgues
- itcnt=0
- ! iteration loop
- do itcnt=1,11
- es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
- qs=0.622*es/(p-es)
- pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
- thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
- f1=thtgs-thes
- if(abs(f1).lt.toler)then
- exit
- endif
- ! itcnt=itcnt+1
- dt=f1*(t1-t0)/(f1-f0)
- t0=t1
- f0=f1
- t1=t1-dt
- enddo
- ttab(it,kp)=t1
- qstab(it,kp)=qs
- enddo
- enddo
- !
- ! lookup table for tlog(emix/aliq)
- !
- ! set up intial values for lookup tables
- !
- astrt=1.e-3
- ainc=0.075
- !
- a1=astrt-ainc
- do i=1,200
- a1=a1+ainc
- alu(i)=alog(a1)
- enddo
- !
- END SUBROUTINE KF_LUTAB
- END MODULE module_cu_kfeta