UPP  001
 All Data Structures Files Functions Pages
UPP_PHYSICS.f
Go to the documentation of this file.
1 
5 
27  module upp_physics
28 
29  implicit none
30 
31  private
32 
33  public :: calcape, calcape2
34  public :: calrh
35  public :: calrh_gfs, calrh_gsd, calrh_nam
36  public :: calrh_pw
37  public :: fpvsnew
38  public :: tvirtual
39 
40  contains
41 !
42 !-------------------------------------------------------------------------------------
43 !
44  SUBROUTINE calrh(P1,T1,Q1,RH)
45 
46  use ctlblk_mod, only: im, jsta, jend, modelname
47  implicit none
48 
49  REAL,dimension(IM,jsta:jend),intent(in) :: p1,t1
50  REAL,dimension(IM,jsta:jend),intent(inout) :: q1
51  REAL,dimension(IM,jsta:jend),intent(out) :: rh
52 
53  IF(modelname == 'RAPR')THEN
54  CALL calrh_gsd(p1,t1,q1,rh)
55  ELSE
56  CALL calrh_nam(p1,t1,q1,rh)
57  END IF
58 
59  END SUBROUTINE calrh
60 !
61 !-------------------------------------------------------------------------------------
62 !
91  SUBROUTINE calrh_nam(P1,T1,Q1,RH)
92  use params_mod, only: pq0, a2, a3, a4, rhmin
93  use ctlblk_mod, only: jsta, jend, spval, im
94 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95  implicit none
96 !
97 ! SET PARAMETER.
98 !
99 ! DECLARE VARIABLES.
100 !
101  REAL,dimension(IM,jsta:jend),intent(in) :: p1,t1
102  REAL,dimension(IM,jsta:jend),intent(inout) :: q1
103  REAL,dimension(IM,jsta:jend),intent(out) :: rh
104  REAL qc
105  integer i,j
106 !***************************************************************
107 !
108 ! START CALRH.
109 !
110  DO j=jsta,jend
111  DO i=1,im
112  IF (t1(i,j) < spval) THEN
113  IF (abs(p1(i,j)) >= 1) THEN
114  qc = pq0/p1(i,j)*exp(a2*(t1(i,j)-a3)/(t1(i,j)-a4))
115 !
116  rh(i,j) = q1(i,j)/qc
117 !
118 ! BOUNDS CHECK
119 !
120  IF (rh(i,j) > 1.0) THEN
121  rh(i,j) = 1.0
122  q1(i,j) = rh(i,j)*qc
123  ENDIF
124  IF (rh(i,j) < rhmin) THEN !use smaller RH limit for stratosphere
125  rh(i,j) = rhmin
126  q1(i,j) = rh(i,j)*qc
127  ENDIF
128 !
129  ENDIF
130  ELSE
131  rh(i,j) = spval
132  ENDIF
133  ENDDO
134  ENDDO
135 !
136 !
137 ! END SUBROUTINE CALRH
138  END SUBROUTINE calrh_nam
139 !
140 !-------------------------------------------------------------------------------------
141 !
171  SUBROUTINE calrh_gfs(P1,T1,Q1,RH)
172  use params_mod, only: rhmin
173  use ctlblk_mod, only: jsta, jend, spval, im
174 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175  implicit none
176 !
177  real,parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K)
178  real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
179  real,parameter:: con_eps =con_rd/con_rv
180  real,parameter:: con_epsm1 =con_rd/con_rv-1
181 ! real,external::FPVSNEW
182 
183 ! INTERFACE
184 ! ELEMENTAL FUNCTION FPVSNEW (t)
185 ! REAL FPVSNEW
186 ! REAL, INTENT(IN) :: t
187 ! END FUNCTION FPVSNEW
188 ! END INTERFACE
189 !
190  REAL,dimension(IM,jsta:jend),intent(in) :: p1,t1
191  REAL,dimension(IM,jsta:jend),intent(inout):: q1,rh
192  REAL es,qc
193  integer :: i,j
194 !***************************************************************
195 !
196 ! START CALRH.
197 !
198 !$omp parallel do private(i,j,es,qc)
199  DO j=jsta,jend
200  DO i=1,im
201  IF (t1(i,j) < spval .AND. p1(i,j) < spval.AND.q1(i,j)/=spval) THEN
202 ! IF (ABS(P1(I,J)) > 1.0) THEN
203 ! IF (P1(I,J) > 1.0) THEN
204  IF (p1(i,j) >= 1.0) THEN
205  es = min(fpvsnew(t1(i,j)),p1(i,j))
206  qc = con_eps*es/(p1(i,j)+con_epsm1*es)
207 
208 ! QC=PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4))
209 
210  rh(i,j) = min(1.0,max(q1(i,j)/qc,rhmin))
211  q1(i,j) = rh(i,j)*qc
212 
213 ! BOUNDS CHECK
214 !
215 ! IF (RH(I,J) > 1.0) THEN
216 ! RH(I,J) = 1.0
217 ! Q1(I,J) = RH(I,J)*QC
218 ! ELSEIF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere
219 ! RH(I,J) = RHmin
220 ! Q1(I,J) = RH(I,J)*QC
221 ! ENDIF
222 
223  ENDIF
224  ELSE
225  rh(i,j) = spval
226  ENDIF
227  ENDDO
228  ENDDO
229 
230  END SUBROUTINE calrh_gfs
231 !
232 !-------------------------------------------------------------------------------------
233 !
234  SUBROUTINE calrh_gsd(P1,T1,Q1,RHB)
235 !
236 ! Algorithm use at GSD for RUC and Rapid Refresh
237 !------------------------------------------------------------------
238 !
239 
240  use ctlblk_mod, only: jsta, jend, im, spval
241 
242  implicit none
243 
244  integer :: j, i
245  real :: tx, pol, esx, es, e
246  real, dimension(im,jsta:jend) :: p1, t1, q1, rhb
247 
248 
249  DO j=jsta,jend
250  DO i=1,im
251  IF (t1(i,j) < spval .AND. p1(i,j) < spval .AND. q1(i,j) < spval) THEN
252 ! - compute relative humidity
253  tx=t1(i,j)-273.15
254  pol = 0.99999683 + tx*(-0.90826951e-02 + &
255  tx*(0.78736169e-04 + tx*(-0.61117958e-06 + &
256  tx*(0.43884187e-08 + tx*(-0.29883885e-10 + &
257  tx*(0.21874425e-12 + tx*(-0.17892321e-14 + &
258  tx*(0.11112018e-16 + tx*(-0.30994571e-19)))))))))
259  esx = 6.1078/pol**8
260 
261  es = esx
262  e = p1(i,j)/100.*q1(i,j)/(0.62197+q1(i,j)*0.37803)
263  rhb(i,j) = min(1.,e/es)
264  ELSE
265  rhb(i,j) = spval
266  ENDIF
267  ENDDO
268  ENDDO
269 
270  END SUBROUTINE calrh_gsd
271 !
272 !-------------------------------------------------------------------------------------
273 !
274  SUBROUTINE calrh_pw(RHPW)
275 !
276 ! Algorithm use at GSD for RUC and Rapid Refresh
277 !------------------------------------------------------------------
278 !
279 
280  use vrbls3d, only: q, pmid, t
281  use params_mod, only: g
282  use ctlblk_mod, only: lm, jsta, jend, lm, im, spval
283 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284  implicit none
285 
286  real,PARAMETER :: svp1=6.1153,svp2=17.67,svp3=29.65
287 
288  REAL, dimension(im,jsta:jend):: pw, pw_sat, rhpw
289  REAL deltp,sh,qv,temp,es,qs,qv_sat
290  integer i,j,l,k,ka,kb
291 
292  pw = 0.
293  pw_sat = 0.
294  rhpw = 0.
295 
296  DO l=1,lm
297  k=lm-l+1
298  DO j=jsta,jend
299  DO i=1,im
300 ! -- use specific humidity for PW calculation
301  if(t(i,j,k)<spval.and.q(i,j,k)<spval) then
302  sh = q(i,j,k)
303  qv = sh/(1.-sh)
304  ka = max(1,k-1)
305  kb = min(lm,k+1)
306 
307 ! assumes that P is in mb at this point - be careful!
308  deltp = 0.5*(pmid(i,j,kb)-pmid(i,j,ka))
309  pw(i,j) = pw(i,j) + sh *deltp/g
310 
311 !Csgb -- Add more for RH w.r.t. PW-sat
312 
313  temp = t(i,j,k)
314 ! --- use saturation mixing ratio w.r.t. water here
315 ! for this check.
316  es = svp1*exp(svp2*(temp-273.15)/(temp-svp3))
317 ! -- get saturation specific humidity (w.r.t. total air)
318  qs = 0.62198*es/(pmid(i,j,k)*1.e-2-0.37802*es)
319 ! -- get saturation mixing ratio (w.r.t. dry air)
320  qv_sat = qs/(1.-qs)
321 
322  pw_sat(i,j) = pw_sat(i,j) + max(sh,qs)*deltp/g
323 
324  if (i==120 .and. j==120 ) &
325  write (6,*)'pw-sat', temp, sh, qs, pmid(i,j,kb) &
326  ,pmid(i,j,ka),pw(i,j),pw_sat(i,j)
327 
328 !sgb - This IS RH w.r.t. PW-sat.
329  rhpw(i,j) = min(1.,pw(i,j) / pw_sat(i,j)) * 100.
330  else
331  rhpw(i,j) = spval
332  endif
333  ENDDO
334  ENDDO
335  ENDDO
336 
337  END SUBROUTINE calrh_pw
338 !
339 !-------------------------------------------------------------------------------------
340 !
341  elemental function fpvsnew(t)
342 
364  implicit none
365  integer,parameter:: nxpvs=7501
366  real,parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt
367  real,parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt
368  real,parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K)
369  real,parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq
370  real,parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond
371  real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
372  real,parameter:: con_csol =2.1060e+3 ! spec heat H2O ice
373  real,parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion
374  real,parameter:: tliq=con_ttp
375  real,parameter:: tice=con_ttp-20.0
376  real,parameter:: dldtl=con_cvap-con_cliq
377  real,parameter:: heatl=con_hvap
378  real,parameter:: xponal=-dldtl/con_rv
379  real,parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
380  real,parameter:: dldti=con_cvap-con_csol
381  real,parameter:: heati=con_hvap+con_hfus
382  real,parameter:: xponai=-dldti/con_rv
383  real,parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
384  real tr,w,pvl,pvi
385  real fpvsnew
386  real,intent(in):: t
387  integer jx
388  real xj,x,tbpvs(nxpvs),xp1
389  real xmin,xmax,xinc,c2xpvs,c1xpvs
390 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391  xmin=180.0
392  xmax=330.0
393  xinc=(xmax-xmin)/(nxpvs-1)
394 ! c1xpvs=1.-xmin/xinc
395  c2xpvs=1./xinc
396  c1xpvs=1.-xmin*c2xpvs
397 ! xj=min(max(c1xpvs+c2xpvs*t,1.0),real(nxpvs,krealfp))
398  xj=min(max(c1xpvs+c2xpvs*t,1.0),float(nxpvs))
399  jx=min(xj,float(nxpvs)-1.0)
400  x=xmin+(jx-1)*xinc
401 
402  tr=con_ttp/x
403  if(x>=tliq) then
404  tbpvs(jx)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
405  elseif(x<tice) then
406  tbpvs(jx)=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
407  else
408  w=(t-tice)/(tliq-tice)
409  pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
410  pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
411  tbpvs(jx)=w*pvl+(1.-w)*pvi
412  endif
413 
414  xp1=xmin+(jx-1+1)*xinc
415 
416  tr=con_ttp/xp1
417  if(xp1>=tliq) then
418  tbpvs(jx+1)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
419  elseif(xp1<tice) then
420  tbpvs(jx+1)=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
421  else
422  w=(t-tice)/(tliq-tice)
423  pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
424  pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
425  tbpvs(jx+1)=w*pvl+(1.-w)*pvi
426  endif
427 
428  fpvsnew=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
429 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430  end function fpvsnew
431 !
432 !-------------------------------------------------------------------------------------
523  SUBROUTINE calcape(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, &
524  cins,pparc,zeql,thund)
525  use vrbls3d, only: pmid, t, q, zint
526  use vrbls2d, only: teql,ieql
527  use masks, only: lmh
528  use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, &
529  oneps, g
530  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
531  plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, &
532  itbq, jtbq, rdpq, the0q, stheq, rdtheq
533  use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, me, spval
534 !
535 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
536  implicit none
537 !
538 ! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980).
539  real,PARAMETER :: ismthp=2,ismtht=2,ismthq=2
540 !
541 ! DECLARE VARIABLES.
542 !
543  integer,intent(in) :: itype
544  real,intent(in) :: dpbnd
545  integer, dimension(IM,Jsta:jend),intent(in) :: l1d
546  real, dimension(IM,Jsta:jend),intent(in) :: p1d,t1d
547  real, dimension(IM,jsta:jend),intent(inout) :: q1d,cape,cins,pparc,zeql
548 !
549  integer, dimension(im,jsta:jend) :: iptb, ithtb, parcel, klres, khres, lcl, idx
550 !
551  real, dimension(im,jsta:jend) :: thesp, psp, cape20, qq, pp, thund
552  REAL, ALLOCATABLE :: tpar(:,:,:)
553 
554  LOGICAL thunder(im,jsta:jend), needthun
555  real psfck,pkl,tbtk,qbtk,apebtk,tthbtk,tthk,apespk,tpspk, &
556  bqs00k,sqs00k,bqs10k,sqs10k,bqk,sqk,tqk,presk,gdzkl,thetap, &
557  thetaa,p00k,p10k,p01k,p11k,tthesk,esatp,qsatp,tvp,tv
558 ! real,external :: fpvsnew
559  integer i,j,l,knuml,knumh,lbeg,lend,iq, kb,ittbk
560 
561 ! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK
562 !
563 !**************************************************************
564 ! START CALCAPE HERE.
565 !
566  ALLOCATE(tpar(im,jsta_2l:jend_2u,lm))
567 !
568 ! COMPUTE CAPE/CINS
569 !
570 ! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
571 ! G * (LN(THETAP) - LN(THETAA)) * DZ
572 !
573 ! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
574 !
575 ! WHERE:
576 ! THETAP IS THE PARCEL THETA
577 ! THETAA IS THE AMBIENT THETA
578 ! DZ IS THE THICKNESS OF THE LAYER
579 !
580 ! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
581 ! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
582 !
583 ! IEQL = EQ LEVEL
584 ! P_thetaemax - real pressure of theta-e max parcel (Pa)
585 !
586 ! INITIALIZE CAPE AND CINS ARRAYS
587 !
588 !$omp parallel do
589  DO j=jsta,jend
590  DO i=1,im
591  cape(i,j) = d00
592  cape20(i,j) = d00
593  cins(i,j) = d00
594  lcl(i,j) = 0
595  thesp(i,j) = d00
596  ieql(i,j) = lm
597  parcel(i,j) = lm
598  psp(i,j) = d00
599  pparc(i,j) = d00
600  thunder(i,j) = .true.
601  ENDDO
602  ENDDO
603 !
604 !$omp parallel do
605  DO l=1,lm
606  DO j=jsta,jend
607  DO i=1,im
608  tpar(i,j,l) = d00
609  ENDDO
610  ENDDO
611  ENDDO
612 !
613 ! TYPE 2 CAPE/CINS:
614 ! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
615 ! ARE DUMMY ARRAYS.
616 !
617  IF (itype == 2) THEN
618 !$omp parallel do private(i,j)
619  DO j=jsta,jend
620  DO i=1,im
621  q1d(i,j) = min(max(h1m12,q1d(i,j)),h99999)
622  ENDDO
623  ENDDO
624  ENDIF
625 !-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
626 !-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
627 !--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
628 
629  DO kb=1,lm
630 !hc IF (ITYPE==2.AND.KB>1) cycle
631  IF (itype == 1 .OR. (itype == 2 .AND. kb == 1)) THEN
632 
633 !$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, &
634 !$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, &
635 !$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk)
636  DO j=jsta,jend
637  DO i=1,im
638  psfck = pmid(i,j,nint(lmh(i,j)))
639  pkl = pmid(i,j,kb)
640  IF(psfck<spval.and.pkl<spval)THEN
641 
642 !hc IF (ITYPE==1.AND.(PKL<PSFCK-DPBND.OR.PKL>PSFCK)) cycle
643  IF (itype ==2 .OR. &
644  (itype == 1 .AND. (pkl >= psfck-dpbnd .AND. pkl <= psfck)))THEN
645  IF (itype == 1) THEN
646  tbtk = t(i,j,kb)
647  qbtk = max(0.0, q(i,j,kb))
648  apebtk = (h10e5/pkl)**capa
649  ELSE
650  pkl = p1d(i,j)
651  tbtk = t1d(i,j)
652  qbtk = max(0.0, q1d(i,j))
653  apebtk = (h10e5/pkl)**capa
654  ENDIF
655 
656 !----------Breogan Gomez - 2009-02-06
657 ! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK
658 ! and a floating invalid.
659 
660 ! if(QBTK < 0) QBTK = 0
661 
662 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
663  tthbtk = tbtk*apebtk
664  tthk = (tthbtk-thl)*rdth
665  qq(i,j) = tthk - aint(tthk)
666  ittbk = int(tthk) + 1
667 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
668  IF(ittbk < 1) THEN
669  ittbk = 1
670  qq(i,j) = d00
671  ENDIF
672  IF(ittbk >= jtb) THEN
673  ittbk = jtb-1
674  qq(i,j) = d00
675  ENDIF
676 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
677  bqs00k = qs0(ittbk)
678  sqs00k = sqs(ittbk)
679  bqs10k = qs0(ittbk+1)
680  sqs10k = sqs(ittbk+1)
681 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
682  bqk = (bqs10k-bqs00k)*qq(i,j) + bqs00k
683  sqk = (sqs10k-sqs00k)*qq(i,j) + sqs00k
684  tqk = (qbtk-bqk)/sqk*rdq
685  pp(i,j) = tqk-aint(tqk)
686  iq = int(tqk)+1
687 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
688  IF(iq < 1) THEN
689  iq = 1
690  pp(i,j) = d00
691  ENDIF
692  IF(iq >= itb) THEN
693  iq = itb-1
694  pp(i,j) = d00
695  ENDIF
696 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
697  p00k = ptbl(iq ,ittbk )
698  p10k = ptbl(iq+1,ittbk )
699  p01k = ptbl(iq ,ittbk+1)
700  p11k = ptbl(iq+1,ittbk+1)
701 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
702  tpspk = p00k + (p10k-p00k)*pp(i,j) + (p01k-p00k)*qq(i,j) &
703  + (p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
704 
705 !!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA
706  if (tpspk > 1.0e-3) then
707  apespk = (max(0.,h10e5/ tpspk))**capa
708  else
709  apespk = 0.0
710  endif
711 
712  tthesk = tthbtk * exp(elocp*qbtk*apespk/tthbtk)
713 !--------------CHECK FOR MAXIMUM THETA E--------------------------------
714  IF(tthesk > thesp(i,j)) THEN
715  psp(i,j) = tpspk
716  thesp(i,j) = tthesk
717  parcel(i,j) = kb
718  ENDIF
719  END IF
720  ENDIF !end PSFCK<spval.and.PKL<spval
721  ENDDO ! I loop
722  ENDDO ! J loop
723  END IF
724  ENDDO ! KB loop
725 
726 !----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED
727 !$omp parallel do private(i,j)
728  DO j=jsta,jend
729  DO i=1,im
730  pparc(i,j) = pmid(i,j,parcel(i,j))
731  ENDDO
732  ENDDO
733 !
734 !-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
735 !-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
736 !-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
737  DO l=1,lm
738 !$omp parallel do private(i,j)
739  DO j=jsta,jend
740  DO i=1,im
741  IF (pmid(i,j,l) < psp(i,j)) lcl(i,j) = l+1
742  ENDDO
743  ENDDO
744  ENDDO
745 !$omp parallel do private(i,j)
746  DO j=jsta,jend
747  DO i=1,im
748  IF (lcl(i,j) > nint(lmh(i,j))) lcl(i,j) = nint(lmh(i,j))
749  IF (itype > 2) THEN
750  IF (t(i,j,lcl(i,j)) < 263.15) THEN
751  thunder(i,j) = .false.
752  ENDIF
753  ENDIF
754  ENDDO
755  ENDDO
756 !-----------------------------------------------------------------------
757 !---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
758 !-----------------------------------------------------------------------
759 
760  DO l=lm,1,-1
761 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
762  knuml = 0
763  knumh = 0
764  DO j=jsta,jend
765  DO i=1,im
766  klres(i,j) = 0
767  khres(i,j) = 0
768  IF(l <= lcl(i,j)) THEN
769  IF(pmid(i,j,l) < plq)THEN
770  knuml = knuml + 1
771  klres(i,j) = 1
772  ELSE
773  knumh = knumh + 1
774  khres(i,j) = 1
775  ENDIF
776  ENDIF
777  ENDDO
778  ENDDO
779 !***
780 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
781 !**
782  IF(knuml > 0) THEN
783  CALL ttblex(tpar(1,jsta_2l,l),ttbl,itb,jtb,klres &
784  , pmid(1,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
785  , rdthe,thesp,iptb,ithtb)
786  ENDIF
787 !***
788 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
789 !**
790  IF(knumh > 0) THEN
791  CALL ttblex(tpar(1,jsta_2l,l),ttblq,itbq,jtbq,khres &
792  , pmid(1,jsta_2l,l),plq,qq,pp,rdpq &
793  ,the0q,stheq,rdtheq,thesp,iptb,ithtb)
794  ENDIF
795 
796 !------------SEARCH FOR EQ LEVEL----------------------------------------
797 !$omp parallel do private(i,j)
798  DO j=jsta,jend
799  DO i=1,im
800  IF(khres(i,j) > 0) THEN
801  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
802  ENDIF
803  ENDDO
804  ENDDO
805 !
806 !$omp parallel do private(i,j)
807  DO j=jsta,jend
808  DO i=1,im
809  IF(klres(i,j) > 0) THEN
810  IF(tpar(i,j,l) > t(i,j,l) .AND. &
811  pmid(i,j,l)>100.) ieql(i,j) = l
812  ENDIF
813  ENDDO
814  ENDDO
815 !-----------------------------------------------------------------------
816  ENDDO ! end of do l=lm,1,-1 loop
817 !------------COMPUTE CAPE AND CINS--------------------------------------
818  lbeg = 1000
819  lend = 0
820  DO j=jsta,jend
821  DO i=1,im
822  lbeg = min(ieql(i,j),lbeg)
823  lend = max(lcl(i,j),lend)
824  ENDDO
825  ENDDO
826 !
827 !$omp parallel do private(i,j)
828  DO j=jsta,jend
829  DO i=1,im
830  IF(t(i,j,ieql(i,j)) > 255.65) THEN
831  thunder(i,j) = .false.
832  ENDIF
833  ENDDO
834  ENDDO
835 !
836  DO l=lbeg,lend
837 
838 !$omp parallel do private(i,j)
839  DO j=jsta,jend
840  DO i=1,im
841  idx(i,j) = 0
842  IF(l >= ieql(i,j).AND.l <= lcl(i,j)) THEN
843  idx(i,j) = 1
844  ENDIF
845  ENDDO
846  ENDDO
847 !
848 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv)
849  DO j=jsta,jend
850  DO i=1,im
851  IF(idx(i,j) > 0) THEN
852  presk = pmid(i,j,l)
853  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
854  esatp = min(fpvsnew(tpar(i,j,l)),presk)
855  qsatp = eps*esatp/(presk-esatp*oneps)
856 ! TVP = TPAR(I,J,L)*(1+0.608*QSATP)
857  tvp = tvirtual(tpar(i,j,l),qsatp)
858  thetap = tvp*(h10e5/presk)**capa
859 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
860  tv = tvirtual(t(i,j,l),q(i,j,l))
861  thetaa = tv*(h10e5/presk)**capa
862  IF(thetap < thetaa) THEN
863  cins(i,j) = cins(i,j) + (log(thetap)-log(thetaa))*gdzkl
864  ELSEIF(thetap > thetaa) THEN
865  cape(i,j) = cape(i,j) + (log(thetap)-log(thetaa))*gdzkl
866  IF (thunder(i,j) .AND. t(i,j,l) < 273.15 &
867  .AND. t(i,j,l) > 253.15) THEN
868  cape20(i,j) = cape20(i,j) + (log(thetap)-log(thetaa))*gdzkl
869  ENDIF
870  ENDIF
871  ENDIF
872  ENDDO
873  ENDDO
874  ENDDO
875 !
876 ! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
877 ! LIMIT OF 0.0 ON CINS.
878 !
879 !$omp parallel do private(i,j)
880  DO j=jsta,jend
881  DO i=1,im
882  cape(i,j) = max(d00,cape(i,j))
883  cins(i,j) = min(cins(i,j),d00)
884 ! add equillibrium height
885  zeql(i,j) = zint(i,j,ieql(i,j))
886  teql(i,j) = t(i,j,ieql(i,j))
887  IF (cape20(i,j) < 75.) THEN
888  thunder(i,j) = .false.
889  ENDIF
890  IF (thunder(i,j)) THEN
891  thund(i,j) = 1.0
892  ELSE
893  thund(i,j) = 0.0
894  ENDIF
895  ENDDO
896  ENDDO
897 !
898  DEALLOCATE(tpar)
899 !
900  END SUBROUTINE calcape
901 !
902 !-------------------------------------------------------------------------------------
998  SUBROUTINE calcape2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, &
999  cape,cins,lfc,esrhl,esrhh, &
1000  dcape,dgld,esp)
1001  use vrbls3d, only: pmid, t, q, zint
1002  use vrbls2d, only: fis,ieql
1003  use gridspec_mod, only: gridtype
1004  use masks, only: lmh
1005  use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, &
1006  oneps, g, tfrz
1007  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
1008  plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, &
1009  itbq, jtbq, rdpq, the0q, stheq, rdtheq
1010  use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, jm, me, jsta_m, jend_m, spval
1011 !
1012 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1013  implicit none
1014 !
1015 ! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980).
1016  real,PARAMETER :: ismthp=2,ismtht=2,ismthq=2
1017 !
1018 ! DECLARE VARIABLES.
1019 !
1020  integer,intent(in) :: itype
1021  real,intent(in) :: dpbnd
1022  integer, dimension(IM,Jsta:jend),intent(in) :: l1d
1023  real, dimension(IM,Jsta:jend),intent(in) :: p1d,t1d
1024 ! real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL
1025  real, dimension(IM,jsta:jend),intent(inout) :: q1d,cape,cins
1026  real, dimension(IM,jsta:jend) :: pparc,zeql
1027  real, dimension(IM,jsta:jend),intent(inout) :: lfc,esrhl,esrhh
1028  real, dimension(IM,jsta:jend),intent(inout) :: dcape,dgld,esp
1029  integer, dimension(im,jsta:jend) ::l12,l17,l3km
1030 !
1031  integer, dimension(im,jsta:jend) :: iptb, ithtb, parcel, klres, khres, lcl, idx
1032 !
1033  real, dimension(im,jsta:jend) :: thesp, psp, cape20, qq, pp, thund
1034  integer, dimension(im,jsta:jend) :: parcel2
1035  real, dimension(im,jsta:jend) :: thesp2,psp2
1036  real, dimension(im,jsta:jend) :: cape4,cins4
1037  REAL, ALLOCATABLE :: tpar(:,:,:)
1038  REAL, ALLOCATABLE :: tpar2(:,:,:)
1039 
1040  LOGICAL thunder(im,jsta:jend), needthun
1041  real psfck,pkl,tbtk,qbtk,apebtk,tthbtk,tthk,apespk,tpspk, &
1042  bqs00k,sqs00k,bqs10k,sqs10k,bqk,sqk,tqk,presk,gdzkl,thetap, &
1043  thetaa,p00k,p10k,p01k,p11k,tthesk,esatp,qsatp,tvp,tv
1044  real presk2, esatp2, qsatp2, tvp2, thetap2, tv2, thetaa2
1045 ! real,external :: fpvsnew
1046  integer i,j,l,knuml,knumh,lbeg,lend,iq, kb,ittbk
1047  integer ie,iw,jn,js,ive(jm),ivw(jm),jvn,jvs
1048  integer istart,istop,jstart,jstop
1049  real, dimension(IM,jsta:jend) :: htsfc
1050 
1051 ! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK
1052 !
1053 !**************************************************************
1054 ! START CALCAPE HERE.
1055 !
1056  ALLOCATE(tpar(im,jsta_2l:jend_2u,lm))
1057  ALLOCATE(tpar2(im,jsta_2l:jend_2u,lm))
1058 !
1059 ! COMPUTE CAPE/CINS
1060 !
1061 ! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
1062 ! G * (LN(THETAP) - LN(THETAA)) * DZ
1063 !
1064 ! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
1065 !
1066 ! WHERE:
1067 ! THETAP IS THE PARCEL THETA
1068 ! THETAA IS THE AMBIENT THETA
1069 ! DZ IS THE THICKNESS OF THE LAYER
1070 !
1071 ! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
1072 ! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
1073 !
1074 ! IEQL = EQ LEVEL
1075 ! P_thetaemax - real pressure of theta-e max parcel (Pa)
1076 !
1077 ! INITIALIZE CAPE AND CINS ARRAYS
1078 !
1079 !$omp parallel do
1080  DO j=jsta,jend
1081  DO i=1,im
1082  cape(i,j) = d00
1083  cape20(i,j) = d00
1084  cape4(i,j) = d00
1085  cins(i,j) = d00
1086  cins4(i,j) = d00
1087  lcl(i,j) = 0
1088  thesp(i,j) = d00
1089  ieql(i,j) = lm
1090  parcel(i,j) = lm
1091  psp(i,j) = d00
1092  pparc(i,j) = d00
1093  thunder(i,j) = .true.
1094  lfc(i,j) = d00
1095  esrhl(i,j) = d00
1096  esrhh(i,j) = d00
1097  dcape(i,j) = d00
1098  dgld(i,j) = d00
1099  esp(i,j) = d00
1100  thesp2(i,j) = 1000.
1101  psp2(i,j) = d00
1102  parcel2(i,j) = lm
1103  ENDDO
1104  ENDDO
1105 !
1106 !$omp parallel do
1107  DO l=1,lm
1108  DO j=jsta,jend
1109  DO i=1,im
1110  tpar(i,j,l) = d00
1111  tpar2(i,j,l) = d00
1112  ENDDO
1113  ENDDO
1114  ENDDO
1115 !
1116 ! FIND SURFACE HEIGHT
1117 !
1118  IF(gridtype == 'E')THEN
1119  jvn = 1
1120  jvs = -1
1121  do j=jsta,jend
1122  ive(j) = mod(j,2)
1123  ivw(j) = ive(j)-1
1124  enddo
1125  istart = 2
1126  istop = im-1
1127  jstart = jsta_m
1128  jstop = jend_m
1129  ELSE IF(gridtype == 'B')THEN
1130  jvn = 1
1131  jvs = 0
1132  do j=jsta,jend
1133  ive(j)=1
1134  ivw(j)=0
1135  enddo
1136  istart = 2
1137  istop = im-1
1138  jstart = jsta_m
1139  jstop = jend_m
1140  ELSE
1141  jvn = 0
1142  jvs = 0
1143  do j=jsta,jend
1144  ive(j) = 0
1145  ivw(j) = 0
1146  enddo
1147  istart = 1
1148  istop = im
1149  jstart = jsta
1150  jstop = jend
1151  END IF
1152 !!$omp parallel do private(htsfc,ie,iw)
1153  IF(gridtype /= 'A') CALL exch(fis(1:im,jsta:jend))
1154  DO j=jstart,jstop
1155  DO i=istart,istop
1156  ie = i+ive(j)
1157  iw = i+ivw(j)
1158  jn = j+jvn
1159  js = j+jvs
1160 !mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+
1161 !mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25
1162 !mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT
1163  IF (gridtype=='B')THEN
1164  htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(ie,jn))
1165  ELSE
1166  htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(i,js))
1167  ENDIF
1168  ENDDO
1169  ENDDO
1170 !
1171 ! TYPE 2 CAPE/CINS:
1172 ! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
1173 ! ARE DUMMY ARRAYS.
1174 !
1175  IF (itype == 2) THEN
1176 !$omp parallel do private(i,j)
1177  DO j=jsta,jend
1178  DO i=1,im
1179  q1d(i,j) = min(max(h1m12,q1d(i,j)),h99999)
1180  ENDDO
1181  ENDDO
1182  ENDIF
1183 !-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
1184 !-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
1185 !--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
1186 
1187  DO kb=1,lm
1188 !hc IF (ITYPE==2.AND.KB>1) cycle
1189  IF (itype == 1 .OR. (itype == 2 .AND. kb == 1)) THEN
1190 
1191 !$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, &
1192 !$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, &
1193 !$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk)
1194  DO j=jsta,jend
1195  DO i=1,im
1196  psfck = pmid(i,j,nint(lmh(i,j)))
1197  pkl = pmid(i,j,kb)
1198 
1199 !hc IF (ITYPE==1.AND.(PKL<PSFCK-DPBND.OR.PKL>PSFCK)) cycle
1200  IF (itype ==2 .OR. &
1201  (itype == 1 .AND. (pkl >= psfck-dpbnd .AND. pkl <= psfck)))THEN
1202  IF (itype == 1) THEN
1203  tbtk = t(i,j,kb)
1204  qbtk = max(0.0, q(i,j,kb))
1205  apebtk = (h10e5/pkl)**capa
1206  ELSE
1207  pkl = p1d(i,j)
1208  tbtk = t1d(i,j)
1209  qbtk = max(0.0, q1d(i,j))
1210  apebtk = (h10e5/pkl)**capa
1211  ENDIF
1212 
1213 !----------Breogan Gomez - 2009-02-06
1214 ! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK
1215 ! and a floating invalid.
1216 
1217 ! if(QBTK < 0) QBTK = 0
1218 
1219 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
1220  tthbtk = tbtk*apebtk
1221  tthk = (tthbtk-thl)*rdth
1222  qq(i,j) = tthk - aint(tthk)
1223  ittbk = int(tthk) + 1
1224 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
1225  IF(ittbk < 1) THEN
1226  ittbk = 1
1227  qq(i,j) = d00
1228  ENDIF
1229  IF(ittbk >= jtb) THEN
1230  ittbk = jtb-1
1231  qq(i,j) = d00
1232  ENDIF
1233 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
1234  bqs00k = qs0(ittbk)
1235  sqs00k = sqs(ittbk)
1236  bqs10k = qs0(ittbk+1)
1237  sqs10k = sqs(ittbk+1)
1238 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
1239  bqk = (bqs10k-bqs00k)*qq(i,j) + bqs00k
1240  sqk = (sqs10k-sqs00k)*qq(i,j) + sqs00k
1241  tqk = (qbtk-bqk)/sqk*rdq
1242  pp(i,j) = tqk-aint(tqk)
1243  iq = int(tqk)+1
1244 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
1245  IF(iq < 1) THEN
1246  iq = 1
1247  pp(i,j) = d00
1248  ENDIF
1249  IF(iq >= itb) THEN
1250  iq = itb-1
1251  pp(i,j) = d00
1252  ENDIF
1253 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
1254  p00k = ptbl(iq ,ittbk )
1255  p10k = ptbl(iq+1,ittbk )
1256  p01k = ptbl(iq ,ittbk+1)
1257  p11k = ptbl(iq+1,ittbk+1)
1258 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
1259  tpspk = p00k + (p10k-p00k)*pp(i,j) + (p01k-p00k)*qq(i,j) &
1260  + (p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
1261 
1262 !!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA
1263  if (tpspk > 1.0e-3) then
1264  apespk = (max(0.,h10e5/ tpspk))**capa
1265  else
1266  apespk = 0.0
1267  endif
1268 
1269  tthesk = tthbtk * exp(elocp*qbtk*apespk/tthbtk)
1270 !--------------CHECK FOR MAXIMUM THETA E--------------------------------
1271  IF(tthesk > thesp(i,j)) THEN
1272  psp(i,j) = tpspk
1273  thesp(i,j) = tthesk
1274  parcel(i,j) = kb
1275  ENDIF
1276 !--------------CHECK FOR MINIMUM THETA E--------------------------------
1277  IF(tthesk < thesp2(i,j)) THEN
1278  psp2(i,j) = tpspk
1279  thesp2(i,j) = tthesk
1280  parcel2(i,j) = kb
1281  ENDIF
1282  END IF
1283  ENDDO ! I loop
1284  ENDDO ! J loop
1285  END IF
1286  ENDDO ! KB loop
1287 
1288 !----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED
1289 !$omp parallel do private(i,j)
1290  DO j=jsta,jend
1291  DO i=1,im
1292  pparc(i,j) = pmid(i,j,parcel(i,j))
1293  ENDDO
1294  ENDDO
1295 !
1296 !-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
1297 !-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
1298 !-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
1299  DO l=1,lm
1300 !$omp parallel do private(i,j)
1301  DO j=jsta,jend
1302  DO i=1,im
1303  IF (pmid(i,j,l) < psp(i,j)) lcl(i,j) = l+1
1304  ENDDO
1305  ENDDO
1306  ENDDO
1307 !$omp parallel do private(i,j)
1308  DO j=jsta,jend
1309  DO i=1,im
1310  IF (lcl(i,j) > nint(lmh(i,j))) lcl(i,j) = nint(lmh(i,j))
1311  IF (itype > 2) THEN
1312  IF (t(i,j,lcl(i,j)) < 263.15) THEN
1313  thunder(i,j) = .false.
1314  ENDIF
1315  ENDIF
1316  ENDDO
1317  ENDDO
1318 !-----------------------------------------------------------------------
1319 !---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
1320 !-----------------------------------------------------------------------
1321  DO l=lm,1,-1
1322 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
1323  knuml = 0
1324  knumh = 0
1325  DO j=jsta,jend
1326  DO i=1,im
1327  klres(i,j) = 0
1328  khres(i,j) = 0
1329  IF(l <= lcl(i,j)) THEN
1330  IF(pmid(i,j,l) < plq)THEN
1331  knuml = knuml + 1
1332  klres(i,j) = 1
1333  ELSE
1334  knumh = knumh + 1
1335  khres(i,j) = 1
1336  ENDIF
1337  ENDIF
1338  ENDDO
1339  ENDDO
1340 !***
1341 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
1342 !**
1343  IF(knuml > 0) THEN
1344  CALL ttblex(tpar(1,jsta_2l,l),ttbl,itb,jtb,klres &
1345  , pmid(1,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
1346  , rdthe,thesp,iptb,ithtb)
1347  ENDIF
1348 !***
1349 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
1350 !**
1351  IF(knumh > 0) THEN
1352  CALL ttblex(tpar(1,jsta_2l,l),ttblq,itbq,jtbq,khres &
1353  , pmid(1,jsta_2l,l),plq,qq,pp,rdpq &
1354  ,the0q,stheq,rdtheq,thesp,iptb,ithtb)
1355  ENDIF
1356 
1357 !------------SEARCH FOR EQ LEVEL----------------------------------------
1358 !$omp parallel do private(i,j)
1359  DO j=jsta,jend
1360  DO i=1,im
1361  IF(khres(i,j) > 0) THEN
1362  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
1363  ENDIF
1364  ENDDO
1365  ENDDO
1366 !
1367 !$omp parallel do private(i,j)
1368  DO j=jsta,jend
1369  DO i=1,im
1370  IF(klres(i,j) > 0) THEN
1371  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
1372  ENDIF
1373  ENDDO
1374  ENDDO
1375 !-----------------------------------------------------------------------
1376  ENDDO ! end of do l=lm,1,-1 loop
1377 !------------COMPUTE CAPE AND CINS--------------------------------------
1378  lbeg = 1000
1379  lend = 0
1380  DO j=jsta,jend
1381  DO i=1,im
1382  lbeg = min(ieql(i,j),lbeg)
1383  lend = max(lcl(i,j),lend)
1384  ENDDO
1385  ENDDO
1386 !
1387 !$omp parallel do private(i,j)
1388  DO j=jsta,jend
1389  DO i=1,im
1390  IF(t(i,j,ieql(i,j)) > 255.65) THEN
1391  thunder(i,j) = .false.
1392  ENDIF
1393  ENDDO
1394  ENDDO
1395 !
1396 !reverse L order from bottom up for ESRH calculation
1397 !
1398  esrhh = lcl
1399  esrhl = lcl
1400 ! DO L=LBEG,LEND
1401  DO l=lend,lbeg,-1
1402 
1403 !$omp parallel do private(i,j)
1404  DO j=jsta,jend
1405  DO i=1,im
1406  idx(i,j) = 0
1407  IF(l >= ieql(i,j).AND.l <= lcl(i,j)) THEN
1408  idx(i,j) = 1
1409  ENDIF
1410  ENDDO
1411  ENDDO
1412 !
1413 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv,&
1414 !$omp & presk2,esatp2,qsatp2,tvp2,thetap2,tv2,thetaa2)
1415  DO j=jsta,jend
1416  DO i=1,im
1417  IF(idx(i,j) > 0) THEN
1418  presk = pmid(i,j,l)
1419  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
1420  esatp = min(fpvsnew(tpar(i,j,l)),presk)
1421  qsatp = eps*esatp/(presk-esatp*oneps)
1422 ! TVP = TPAR(I,J,L)*(1+0.608*QSATP)
1423  tvp = tvirtual(tpar(i,j,l),qsatp)
1424  thetap = tvp*(h10e5/presk)**capa
1425 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
1426  tv = tvirtual(t(i,j,l),q(i,j,l))
1427  thetaa = tv*(h10e5/presk)**capa
1428  IF(thetap < thetaa) THEN
1429  cins4(i,j) = cins4(i,j) + (log(thetap)-log(thetaa))*gdzkl
1430  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1431  cins(i,j) = cins(i,j) + (log(thetap)-log(thetaa))*gdzkl
1432  ENDIF
1433  ELSEIF(thetap > thetaa) THEN
1434  cape4(i,j) = cape4(i,j) + (log(thetap)-log(thetaa))*gdzkl
1435  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1436  cape(i,j) = cape(i,j) + (log(thetap)-log(thetaa))*gdzkl
1437  ENDIF
1438  IF (thunder(i,j) .AND. t(i,j,l) < 273.15 &
1439  .AND. t(i,j,l) > 253.15) THEN
1440  cape20(i,j) = cape20(i,j) + (log(thetap)-log(thetaa))*gdzkl
1441  ENDIF
1442  ENDIF
1443 
1444 ! LFC
1445  IF (itype /= 1) THEN
1446  presk2 = pmid(i,j,l+1)
1447  esatp2 = min(fpvsnew(tpar(i,j,l+1)),presk2)
1448  qsatp2 = eps*esatp2/(presk2-esatp2*oneps)
1449 ! TVP2 = TPAR(I,J,L+1)*(1+0.608*QSATP2)
1450  tvp2 = tvirtual(tpar(i,j,l+1),qsatp2)
1451  thetap2 = tvp2*(h10e5/presk2)**capa
1452 ! TV2 = T(I,J,L+1)*(1+0.608*Q(I,J,L+1))
1453  tv2 = tvirtual(t(i,j,l+1),q(i,j,l+1))
1454  thetaa2 = tv2*(h10e5/presk2)**capa
1455  IF(thetap >= thetaa .AND. thetap2 <= thetaa2) THEN
1456  IF(lfc(i,j) == d00)THEN
1457  lfc(i,j) = zint(i,j,l)
1458  ENDIF
1459  ENDIF
1460  ENDIF
1461 !
1462 ! ESRH/CAPE threshold check
1463  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1464  IF(cape4(i,j) >= 100. .AND. cins4(i,j) >= -250.) THEN
1465  IF(esrhl(i,j) == lcl(i,j)) esrhl(i,j)=l
1466  ENDIF
1467  esrhh(i,j)=l
1468  ENDIF
1469 
1470  ENDIF !(IDX(I,J) > 0)
1471  ENDDO
1472  ENDDO
1473  ENDDO
1474 
1475 !$omp parallel do private(i,j)
1476  DO j=jsta,jend
1477  DO i=1,im
1478  IF(esrhh(i,j) > esrhl(i,j)) esrhh(i,j)=ieql(i,j)
1479  ENDDO
1480  ENDDO
1481 !
1482 ! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
1483 ! LIMIT OF 0.0 ON CINS.
1484 ! ENFORCE LFC ABOVE LCL AND BELOW EL
1485 !
1486 !$omp parallel do private(i,j)
1487  DO j=jsta,jend
1488  DO i=1,im
1489  cape(i,j) = max(d00,cape(i,j))
1490  cins(i,j) = min(cins(i,j),d00)
1491 ! equillibrium height
1492  zeql(i,j) = zint(i,j,ieql(i,j))
1493  lfc(i,j) = min(lfc(i,j),zint(i,j,ieql(i,j)))
1494  lfc(i,j) = max(zint(i,j, lcl(i,j)),lfc(i,j))
1495  IF (cape20(i,j) < 75.) THEN
1496  thunder(i,j) = .false.
1497  ENDIF
1498  IF (thunder(i,j)) THEN
1499  thund(i,j) = 1.0
1500  ELSE
1501  thund(i,j) = 0.0
1502  ENDIF
1503  ENDDO
1504  ENDDO
1505 !------------COMPUTE DCAPE--------------------------------------
1506 !-----------------------------------------------------------------------
1507 !---------FIND TEMP OF PARCEL DESCENDED ALONG MOIST ADIABAT (TPAR)---------
1508 !-----------------------------------------------------------------------
1509  IF (itype == 1) THEN
1510 
1511  DO l=lm,1,-1
1512 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
1513  knuml = 0
1514  knumh = 0
1515  DO j=jsta,jend
1516  DO i=1,im
1517  klres(i,j) = 0
1518  khres(i,j) = 0
1519  psfck = pmid(i,j,nint(lmh(i,j)))
1520  pkl = pmid(i,j,l)
1521  IF(pkl >= psfck-dpbnd) THEN
1522  IF(pmid(i,j,l) < plq)THEN
1523  knuml = knuml + 1
1524  klres(i,j) = 1
1525  ELSE
1526  knumh = knumh + 1
1527  khres(i,j) = 1
1528  ENDIF
1529  ENDIF
1530  ENDDO
1531  ENDDO
1532 !***
1533 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
1534 !**
1535  IF(knuml > 0) THEN
1536  CALL ttblex(tpar2(1,jsta_2l,l),ttbl,itb,jtb,klres &
1537  , pmid(1,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
1538  , rdthe,thesp2,iptb,ithtb)
1539  ENDIF
1540 !***
1541 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
1542 !**
1543  IF(knumh > 0) THEN
1544  CALL ttblex(tpar2(1,jsta_2l,l),ttblq,itbq,jtbq,khres &
1545  , pmid(1,jsta_2l,l),plq,qq,pp,rdpq &
1546  , the0q,stheq,rdtheq,thesp2,iptb,ithtb)
1547  ENDIF
1548  ENDDO ! end of do l=lm,1,-1 loop
1549 
1550  lbeg = 1
1551  lend = lm
1552 
1553  DO l=lbeg,lend
1554 !$omp parallel do private(i,j)
1555  DO j=jsta,jend
1556  DO i=1,im
1557  idx(i,j) = 0
1558  IF(l >= parcel2(i,j).AND.l < nint(lmh(i,j))) THEN
1559  idx(i,j) = 1
1560  ENDIF
1561  ENDDO
1562  ENDDO
1563 !
1564 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv)
1565  DO j=jsta,jend
1566  DO i=1,im
1567  IF(idx(i,j) > 0) THEN
1568  presk = pmid(i,j,l)
1569  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
1570  esatp = min(fpvsnew(tpar2(i,j,l)),presk)
1571  qsatp = eps*esatp/(presk-esatp*oneps)
1572 ! TVP = TPAR2(I,J,L)*(1+0.608*QSATP)
1573  tvp = tvirtual(tpar2(i,j,l),qsatp)
1574  thetap = tvp*(h10e5/presk)**capa
1575 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
1576  tv = tvirtual(t(i,j,l),q(i,j,l))
1577  thetaa = tv*(h10e5/presk)**capa
1578  !IF(THETAP < THETAA) THEN
1579  dcape(i,j) = dcape(i,j) + (log(thetap)-log(thetaa))*gdzkl
1580  !ENDIF
1581  ENDIF
1582  ENDDO
1583  ENDDO
1584  ENDDO
1585 
1586 !$omp parallel do private(i,j)
1587  DO j=jsta,jend
1588  DO i=1,im
1589  dcape(i,j) = min(d00,dcape(i,j))
1590  ENDDO
1591  ENDDO
1592 
1593  ENDIF !ITYPE=1 FOR DCAPE
1594 
1595 !
1596 ! Dendritic Growth Layer depth
1597 ! the layer with temperatures from -12 to -17 C in meters
1598 !
1599  l12=lm
1600  l17=lm
1601  DO l=lm,1,-1
1602 !$omp parallel do private(i,j)
1603  DO j=jsta,jend
1604  DO i=1,im
1605  IF(t(i,j,l) <= tfrz-12. .AND. l12(i,j)==lm) l12(i,j)=l
1606  IF(t(i,j,l) <= tfrz-17. .AND. l17(i,j)==lm) l17(i,j)=l
1607  ENDDO
1608  ENDDO
1609  ENDDO
1610 !$omp parallel do private(i,j)
1611  DO j=jsta,jend
1612  DO i=1,im
1613  IF(l12(i,j)/=lm .AND. l17(i,j)/=lm) THEN
1614  dgld(i,j)=zint(i,j,l17(i,j))-zint(i,j,l12(i,j))
1615  dgld(i,j)=max(dgld(i,j),0.)
1616  ENDIF
1617  ENDDO
1618  ENDDO
1619 !
1620 ! Enhanced Stretching Potential
1621 ! ESP = (0-3 km MLCAPE / 50 J kg-1) * ((0-3 km lapse rate - 7.0) / 1.0 C (km-1)
1622 ! https://www.spc.noaa.gov/exper/soundings/help/params4.html
1623 !
1624  l3km=lm
1625  DO l=lm,1,-1
1626 !$omp parallel do private(i,j)
1627  DO j=jsta,jend
1628  DO i=1,im
1629  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) l3km(i,j)=l
1630  ENDDO
1631  ENDDO
1632  ENDDO
1633 !$omp parallel do private(i,j)
1634  DO j=jsta,jend
1635  DO i=1,im
1636  esp(i,j) = (cape(i,j) / 50.) * (t(i,j,lm) - t(i,j,l3km(i,j)) - 7.0)
1637  IF((t(i,j,lm) - t(i,j,l3km(i,j))) < 7.0) esp(i,j) = 0.
1638 ! IF(CAPE(I,J) < 250.) ESP(I,J) = 0.
1639  ENDDO
1640  ENDDO
1641 !
1642  DEALLOCATE(tpar)
1643  DEALLOCATE(tpar2)
1644 !
1645  END SUBROUTINE calcape2
1646 !
1647 !-------------------------------------------------------------------------------------
1648 !
1649 !
1650 !-------------------------------------------------------------------------------------
1651 !
1652  elemental function tvirtual(T,Q)
1653 !
1654 ! COMPUTE VIRTUAL TEMPERATURE
1655 !
1656  IMPLICIT NONE
1657  REAL tvirtual
1658  REAL, INTENT(IN) :: t, q
1659 
1660  tvirtual = t*(1+0.608*q)
1661 
1662  end function tvirtual
1663 !
1664 !-------------------------------------------------------------------------------------
1665 !
1666  end module upp_physics
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:341
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27
lookup_mod
Definition: LOOKUP.f:1