UPP  V11.0.0
 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 :: caldiv
35  public :: calgradps
36  public :: calrh
37  public :: calrh_gfs, calrh_gsd, calrh_nam
38  public :: calrh_pw
39  public :: calvor
40 
41  public :: fpvsnew
42  public :: tvirtual
43 
44  contains
45 !
46 !-------------------------------------------------------------------------------------
47 !
48  SUBROUTINE calrh(P1,T1,Q1,RH)
49 
50  use ctlblk_mod, only: ista, iend, jsta, jend, modelname
51  implicit none
52 
53  REAL,dimension(ista:iend,jsta:jend),intent(in) :: p1,t1
54  REAL,dimension(ista:iend,jsta:jend),intent(inout) :: q1
55  REAL,dimension(ista:iend,jsta:jend),intent(out) :: rh
56 
57  IF(modelname == 'RAPR')THEN
58  CALL calrh_gsd(p1,t1,q1,rh)
59  ELSE
60  CALL calrh_nam(p1,t1,q1,rh)
61  END IF
62 
63  END SUBROUTINE calrh
64 !
65 !-------------------------------------------------------------------------------------
66 !
95  SUBROUTINE calrh_nam(P1,T1,Q1,RH)
96  use params_mod, only: pq0, a2, a3, a4, rhmin
97  use ctlblk_mod, only: ista, iend, jsta, jend, spval
98 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99  implicit none
100 !
101 ! SET PARAMETER.
102 !
103 ! DECLARE VARIABLES.
104 !
105  REAL,dimension(ista:iend,jsta:jend),intent(in) :: p1,t1
106  REAL,dimension(ista:iend,jsta:jend),intent(inout) :: q1
107  REAL,dimension(ista:iend,jsta:jend),intent(out) :: rh
108  REAL qc
109  integer i,j
110 !***************************************************************
111 !
112 ! START CALRH.
113 !
114  DO j=jsta,jend
115  DO i=ista,iend
116  IF (t1(i,j) < spval) THEN
117  IF (abs(p1(i,j)) >= 1) THEN
118  qc = pq0/p1(i,j)*exp(a2*(t1(i,j)-a3)/(t1(i,j)-a4))
119 !
120  rh(i,j) = q1(i,j)/qc
121 !
122 ! BOUNDS CHECK
123 !
124  IF (rh(i,j) > 1.0) THEN
125  rh(i,j) = 1.0
126  q1(i,j) = rh(i,j)*qc
127  ENDIF
128  IF (rh(i,j) < rhmin) THEN !use smaller RH limit for stratosphere
129  rh(i,j) = rhmin
130  q1(i,j) = rh(i,j)*qc
131  ENDIF
132 !
133  ENDIF
134  ELSE
135  rh(i,j) = spval
136  ENDIF
137  ENDDO
138  ENDDO
139 !
140 !
141 ! END SUBROUTINE CALRH
142  END SUBROUTINE calrh_nam
143 !
144 !-------------------------------------------------------------------------------------
145 !
175  SUBROUTINE calrh_gfs(P1,T1,Q1,RH)
176  use params_mod, only: rhmin
177  use ctlblk_mod, only: ista, iend, jsta, jend, spval
178 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179  implicit none
180 !
181  real,parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K)
182  real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
183  real,parameter:: con_eps =con_rd/con_rv
184  real,parameter:: con_epsm1 =con_rd/con_rv-1
185 ! real,external::FPVSNEW
186 
187 ! INTERFACE
188 ! ELEMENTAL FUNCTION FPVSNEW (t)
189 ! REAL FPVSNEW
190 ! REAL, INTENT(IN) :: t
191 ! END FUNCTION FPVSNEW
192 ! END INTERFACE
193 !
194  REAL,dimension(ista:iend,jsta:jend),intent(in) :: p1,t1
195  REAL,dimension(ista:iend,jsta:jend),intent(inout):: q1,rh
196  REAL es,qc
197  integer :: i,j
198 !***************************************************************
199 !
200 ! START CALRH.
201 !
202 !$omp parallel do private(i,j,es,qc)
203  DO j=jsta,jend
204  DO i=ista,iend
205  IF (t1(i,j) < spval .AND. p1(i,j) < spval.AND.q1(i,j)/=spval) THEN
206 ! IF (ABS(P1(I,J)) > 1.0) THEN
207 ! IF (P1(I,J) > 1.0) THEN
208  IF (p1(i,j) >= 1.0) THEN
209  es = min(fpvsnew(t1(i,j)),p1(i,j))
210  qc = con_eps*es/(p1(i,j)+con_epsm1*es)
211 
212 ! QC=PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4))
213 
214  rh(i,j) = min(1.0,max(q1(i,j)/qc,rhmin))
215  q1(i,j) = rh(i,j)*qc
216 
217 ! BOUNDS CHECK
218 !
219 ! IF (RH(I,J) > 1.0) THEN
220 ! RH(I,J) = 1.0
221 ! Q1(I,J) = RH(I,J)*QC
222 ! ELSEIF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere
223 ! RH(I,J) = RHmin
224 ! Q1(I,J) = RH(I,J)*QC
225 ! ENDIF
226 
227  ENDIF
228  ELSE
229  rh(i,j) = spval
230  ENDIF
231  ENDDO
232  ENDDO
233 
234  END SUBROUTINE calrh_gfs
235 !
236 !-------------------------------------------------------------------------------------
237 !
238  SUBROUTINE calrh_gsd(P1,T1,Q1,RHB)
239 !
240 ! Algorithm use at GSD for RUC and Rapid Refresh
241 !------------------------------------------------------------------
242 !
243 
244  use ctlblk_mod, only: ista, iend, jsta, jend, spval
245 
246  implicit none
247 
248  integer :: j, i
249  real :: tx, pol, esx, es, e
250  real, dimension(ista:iend,jsta:jend) :: p1, t1, q1, rhb
251 
252 
253  DO j=jsta,jend
254  DO i=ista,iend
255  IF (t1(i,j) < spval .AND. p1(i,j) < spval .AND. q1(i,j) < spval) THEN
256 ! - compute relative humidity
257  tx=t1(i,j)-273.15
258  pol = 0.99999683 + tx*(-0.90826951e-02 + &
259  tx*(0.78736169e-04 + tx*(-0.61117958e-06 + &
260  tx*(0.43884187e-08 + tx*(-0.29883885e-10 + &
261  tx*(0.21874425e-12 + tx*(-0.17892321e-14 + &
262  tx*(0.11112018e-16 + tx*(-0.30994571e-19)))))))))
263  esx = 6.1078/pol**8
264 
265  es = esx
266  e = p1(i,j)/100.*q1(i,j)/(0.62197+q1(i,j)*0.37803)
267  rhb(i,j) = min(1.,e/es)
268  ELSE
269  rhb(i,j) = spval
270  ENDIF
271  ENDDO
272  ENDDO
273 
274  END SUBROUTINE calrh_gsd
275 !
276 !-------------------------------------------------------------------------------------
277 !
278  SUBROUTINE calrh_pw(RHPW)
279 !
280 ! Algorithm use at GSD for RUC and Rapid Refresh
281 !------------------------------------------------------------------
282 !
283 
284  use vrbls3d, only: q, pmid, t
285  use params_mod, only: g
286  use ctlblk_mod, only: lm, ista, iend, jsta, jend, spval
287 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288  implicit none
289 
290  real,PARAMETER :: svp1=6.1153,svp2=17.67,svp3=29.65
291 
292  REAL, dimension(ista:iend,jsta:jend):: pw, pw_sat, rhpw
293  REAL deltp,sh,qv,temp,es,qs,qv_sat
294  integer i,j,l,k,ka,kb
295 
296  pw = 0.
297  pw_sat = 0.
298  rhpw = 0.
299 
300  DO l=1,lm
301  k=lm-l+1
302  DO j=jsta,jend
303  DO i=ista,iend
304 ! -- use specific humidity for PW calculation
305  if(t(i,j,k)<spval.and.q(i,j,k)<spval) then
306  sh = q(i,j,k)
307  qv = sh/(1.-sh)
308  ka = max(1,k-1)
309  kb = min(lm,k+1)
310 
311 ! assumes that P is in mb at this point - be careful!
312  deltp = 0.5*(pmid(i,j,kb)-pmid(i,j,ka))
313  pw(i,j) = pw(i,j) + sh *deltp/g
314 
315 !Csgb -- Add more for RH w.r.t. PW-sat
316 
317  temp = t(i,j,k)
318 ! --- use saturation mixing ratio w.r.t. water here
319 ! for this check.
320  es = svp1*exp(svp2*(temp-273.15)/(temp-svp3))
321 ! -- get saturation specific humidity (w.r.t. total air)
322  qs = 0.62198*es/(pmid(i,j,k)*1.e-2-0.37802*es)
323 ! -- get saturation mixing ratio (w.r.t. dry air)
324  qv_sat = qs/(1.-qs)
325 
326  pw_sat(i,j) = pw_sat(i,j) + max(sh,qs)*deltp/g
327 
328  if (i==120 .and. j==120 ) &
329  write (6,*)'pw-sat', temp, sh, qs, pmid(i,j,kb) &
330  ,pmid(i,j,ka),pw(i,j),pw_sat(i,j)
331 
332 !sgb - This IS RH w.r.t. PW-sat.
333  rhpw(i,j) = min(1.,pw(i,j) / pw_sat(i,j)) * 100.
334  else
335  rhpw(i,j) = spval
336  endif
337  ENDDO
338  ENDDO
339  ENDDO
340 
341  END SUBROUTINE calrh_pw
342 !
343 !-------------------------------------------------------------------------------------
344 !
345  elemental function fpvsnew(t)
346 
368  implicit none
369  integer,parameter:: nxpvs=7501
370  real,parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt
371  real,parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt
372  real,parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K)
373  real,parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq
374  real,parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond
375  real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
376  real,parameter:: con_csol =2.1060e+3 ! spec heat H2O ice
377  real,parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion
378  real,parameter:: tliq=con_ttp
379  real,parameter:: tice=con_ttp-20.0
380  real,parameter:: dldtl=con_cvap-con_cliq
381  real,parameter:: heatl=con_hvap
382  real,parameter:: xponal=-dldtl/con_rv
383  real,parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
384  real,parameter:: dldti=con_cvap-con_csol
385  real,parameter:: heati=con_hvap+con_hfus
386  real,parameter:: xponai=-dldti/con_rv
387  real,parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
388  real tr,w,pvl,pvi
389  real fpvsnew
390  real,intent(in):: t
391  integer jx
392  real xj,x,tbpvs(nxpvs),xp1
393  real xmin,xmax,xinc,c2xpvs,c1xpvs
394 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395  xmin=180.0
396  xmax=330.0
397  xinc=(xmax-xmin)/(nxpvs-1)
398 ! c1xpvs=1.-xmin/xinc
399  c2xpvs=1./xinc
400  c1xpvs=1.-xmin*c2xpvs
401 ! xj=min(max(c1xpvs+c2xpvs*t,1.0),real(nxpvs,krealfp))
402  xj=min(max(c1xpvs+c2xpvs*t,1.0),float(nxpvs))
403  jx=min(xj,float(nxpvs)-1.0)
404  x=xmin+(jx-1)*xinc
405 
406  tr=con_ttp/x
407  if(x>=tliq) then
408  tbpvs(jx)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
409  elseif(x<tice) then
410  tbpvs(jx)=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
411  else
412  w=(t-tice)/(tliq-tice)
413  pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
414  pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
415  tbpvs(jx)=w*pvl+(1.-w)*pvi
416  endif
417 
418  xp1=xmin+(jx-1+1)*xinc
419 
420  tr=con_ttp/xp1
421  if(xp1>=tliq) then
422  tbpvs(jx+1)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
423  elseif(xp1<tice) then
424  tbpvs(jx+1)=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
425  else
426  w=(t-tice)/(tliq-tice)
427  pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
428  pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
429  tbpvs(jx+1)=w*pvl+(1.-w)*pvi
430  endif
431 
432  fpvsnew=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
433 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434  end function fpvsnew
435 !
436 !-------------------------------------------------------------------------------------
527  SUBROUTINE calcape(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, &
528  cins,pparc,zeql,thund)
529  use vrbls3d, only: pmid, t, q, zint
530  use vrbls2d, only: teql,ieql
531  use masks, only: lmh
532  use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, &
533  oneps, g
534  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
535  plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, &
536  itbq, jtbq, rdpq, the0q, stheq, rdtheq
537  use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, me, spval, &
538  ista_2l, iend_2u, ista, iend
539 !
540 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
541  implicit none
542 !
543 ! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980).
544  real,PARAMETER :: ismthp=2,ismtht=2,ismthq=2
545 !
546 ! DECLARE VARIABLES.
547 !
548  integer,intent(in) :: itype
549  real,intent(in) :: dpbnd
550  integer, dimension(ista:iend,Jsta:jend),intent(in) :: l1d
551  real, dimension(ista:iend,Jsta:jend),intent(in) :: p1d,t1d
552  real, dimension(ista:iend,jsta:jend),intent(inout) :: q1d,cape,cins,pparc,zeql
553 !
554  integer, dimension(ista:iend,jsta:jend) :: iptb, ithtb, parcel, klres, khres, lcl, idx
555 !
556  real, dimension(ista:iend,jsta:jend) :: thesp, psp, cape20, qq, pp, thund
557  REAL, ALLOCATABLE :: tpar(:,:,:)
558 
559  LOGICAL thunder(ista:iend,jsta:jend), needthun
560  real psfck,pkl,tbtk,qbtk,apebtk,tthbtk,tthk,apespk,tpspk, &
561  bqs00k,sqs00k,bqs10k,sqs10k,bqk,sqk,tqk,presk,gdzkl,thetap, &
562  thetaa,p00k,p10k,p01k,p11k,tthesk,esatp,qsatp,tvp,tv
563 ! real,external :: fpvsnew
564  integer i,j,l,knuml,knumh,lbeg,lend,iq, kb,ittbk
565 
566 ! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK
567 !
568 !**************************************************************
569 ! START CALCAPE HERE.
570 !
571  ALLOCATE(tpar(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
572 !
573 ! COMPUTE CAPE/CINS
574 !
575 ! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
576 ! G * (LN(THETAP) - LN(THETAA)) * DZ
577 !
578 ! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
579 !
580 ! WHERE:
581 ! THETAP IS THE PARCEL THETA
582 ! THETAA IS THE AMBIENT THETA
583 ! DZ IS THE THICKNESS OF THE LAYER
584 !
585 ! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
586 ! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
587 !
588 ! IEQL = EQ LEVEL
589 ! P_thetaemax - real pressure of theta-e max parcel (Pa)
590 !
591 ! INITIALIZE CAPE AND CINS ARRAYS
592 !
593 !$omp parallel do
594  DO j=jsta,jend
595  DO i=ista,iend
596  cape(i,j) = d00
597  cape20(i,j) = d00
598  cins(i,j) = d00
599  lcl(i,j) = 0
600  thesp(i,j) = d00
601  ieql(i,j) = lm
602  parcel(i,j) = lm
603  psp(i,j) = d00
604  pparc(i,j) = d00
605  thunder(i,j) = .true.
606  ENDDO
607  ENDDO
608 !
609 !$omp parallel do
610  DO l=1,lm
611  DO j=jsta,jend
612  DO i=ista,iend
613  tpar(i,j,l) = d00
614  ENDDO
615  ENDDO
616  ENDDO
617 !
618 ! TYPE 2 CAPE/CINS:
619 ! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
620 ! ARE DUMMY ARRAYS.
621 !
622  IF (itype == 2) THEN
623 !$omp parallel do private(i,j)
624  DO j=jsta,jend
625  DO i=ista,iend
626  q1d(i,j) = min(max(h1m12,q1d(i,j)),h99999)
627  ENDDO
628  ENDDO
629  ENDIF
630 !-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
631 !-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
632 !--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
633 
634  DO kb=1,lm
635 !hc IF (ITYPE==2.AND.KB>1) cycle
636  IF (itype == 1 .OR. (itype == 2 .AND. kb == 1)) THEN
637 
638 !$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, &
639 !$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, &
640 !$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk)
641  DO j=jsta,jend
642  DO i=ista,iend
643  psfck = pmid(i,j,nint(lmh(i,j)))
644  pkl = pmid(i,j,kb)
645  IF(psfck<spval.and.pkl<spval)THEN
646 
647 !hc IF (ITYPE==1.AND.(PKL<PSFCK-DPBND.OR.PKL>PSFCK)) cycle
648  IF (itype ==2 .OR. &
649  (itype == 1 .AND. (pkl >= psfck-dpbnd .AND. pkl <= psfck)))THEN
650  IF (itype == 1) THEN
651  tbtk = t(i,j,kb)
652  qbtk = max(0.0, q(i,j,kb))
653  apebtk = (h10e5/pkl)**capa
654  ELSE
655  pkl = p1d(i,j)
656  tbtk = t1d(i,j)
657  qbtk = max(0.0, q1d(i,j))
658  apebtk = (h10e5/pkl)**capa
659  ENDIF
660 
661 !----------Breogan Gomez - 2009-02-06
662 ! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK
663 ! and a floating invalid.
664 
665 ! if(QBTK < 0) QBTK = 0
666 
667 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
668  tthbtk = tbtk*apebtk
669  tthk = (tthbtk-thl)*rdth
670  qq(i,j) = tthk - aint(tthk)
671  ittbk = int(tthk) + 1
672 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
673  IF(ittbk < 1) THEN
674  ittbk = 1
675  qq(i,j) = d00
676  ENDIF
677  IF(ittbk >= jtb) THEN
678  ittbk = jtb-1
679  qq(i,j) = d00
680  ENDIF
681 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
682  bqs00k = qs0(ittbk)
683  sqs00k = sqs(ittbk)
684  bqs10k = qs0(ittbk+1)
685  sqs10k = sqs(ittbk+1)
686 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
687  bqk = (bqs10k-bqs00k)*qq(i,j) + bqs00k
688  sqk = (sqs10k-sqs00k)*qq(i,j) + sqs00k
689  tqk = (qbtk-bqk)/sqk*rdq
690  pp(i,j) = tqk-aint(tqk)
691  iq = int(tqk)+1
692 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
693  IF(iq < 1) THEN
694  iq = 1
695  pp(i,j) = d00
696  ENDIF
697  IF(iq >= itb) THEN
698  iq = itb-1
699  pp(i,j) = d00
700  ENDIF
701 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
702  p00k = ptbl(iq ,ittbk )
703  p10k = ptbl(iq+1,ittbk )
704  p01k = ptbl(iq ,ittbk+1)
705  p11k = ptbl(iq+1,ittbk+1)
706 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
707  tpspk = p00k + (p10k-p00k)*pp(i,j) + (p01k-p00k)*qq(i,j) &
708  + (p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
709 
710 !!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA
711  if (tpspk > 1.0e-3) then
712  apespk = (max(0.,h10e5/ tpspk))**capa
713  else
714  apespk = 0.0
715  endif
716 
717  tthesk = tthbtk * exp(elocp*qbtk*apespk/tthbtk)
718 !--------------CHECK FOR MAXIMUM THETA E--------------------------------
719  IF(tthesk > thesp(i,j)) THEN
720  psp(i,j) = tpspk
721  thesp(i,j) = tthesk
722  parcel(i,j) = kb
723  ENDIF
724  END IF
725  ENDIF !end PSFCK<spval.and.PKL<spval
726  ENDDO ! I loop
727  ENDDO ! J loop
728  END IF
729  ENDDO ! KB loop
730 
731 !----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED
732 !$omp parallel do private(i,j)
733  DO j=jsta,jend
734  DO i=ista,iend
735  pparc(i,j) = pmid(i,j,parcel(i,j))
736  ENDDO
737  ENDDO
738 !
739 !-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
740 !-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
741 !-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
742  DO l=1,lm
743 !$omp parallel do private(i,j)
744  DO j=jsta,jend
745  DO i=ista,iend
746  IF (pmid(i,j,l) < psp(i,j)) lcl(i,j) = l+1
747  ENDDO
748  ENDDO
749  ENDDO
750 !$omp parallel do private(i,j)
751  DO j=jsta,jend
752  DO i=ista,iend
753  IF (lcl(i,j) > nint(lmh(i,j))) lcl(i,j) = nint(lmh(i,j))
754  IF (itype > 2) THEN
755  IF (t(i,j,lcl(i,j)) < 263.15) THEN
756  thunder(i,j) = .false.
757  ENDIF
758  ENDIF
759  ENDDO
760  ENDDO
761 !-----------------------------------------------------------------------
762 !---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
763 !-----------------------------------------------------------------------
764 
765  DO l=lm,1,-1
766 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
767  knuml = 0
768  knumh = 0
769  DO j=jsta,jend
770  DO i=ista,iend
771  klres(i,j) = 0
772  khres(i,j) = 0
773  IF(l <= lcl(i,j)) THEN
774  IF(pmid(i,j,l) < plq)THEN
775  knuml = knuml + 1
776  klres(i,j) = 1
777  ELSE
778  knumh = knumh + 1
779  khres(i,j) = 1
780  ENDIF
781  ENDIF
782  ENDDO
783  ENDDO
784 !***
785 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
786 !**
787  IF(knuml > 0) THEN
788  CALL ttblex(tpar(ista_2l,jsta_2l,l),ttbl,itb,jtb,klres &
789  , pmid(ista_2l,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
790  , rdthe,thesp,iptb,ithtb)
791  ENDIF
792 !***
793 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
794 !**
795  IF(knumh > 0) THEN
796  CALL ttblex(tpar(ista_2l,jsta_2l,l),ttblq,itbq,jtbq,khres &
797  , pmid(ista_2l,jsta_2l,l),plq,qq,pp,rdpq &
798  ,the0q,stheq,rdtheq,thesp,iptb,ithtb)
799  ENDIF
800 
801 !------------SEARCH FOR EQ LEVEL----------------------------------------
802 !$omp parallel do private(i,j)
803  DO j=jsta,jend
804  DO i=ista,iend
805  IF(khres(i,j) > 0) THEN
806  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
807  ENDIF
808  ENDDO
809  ENDDO
810 !
811 !$omp parallel do private(i,j)
812  DO j=jsta,jend
813  DO i=ista,iend
814  IF(klres(i,j) > 0) THEN
815  IF(tpar(i,j,l) > t(i,j,l) .AND. &
816  pmid(i,j,l)>100.) ieql(i,j) = l
817  ENDIF
818  ENDDO
819  ENDDO
820 !-----------------------------------------------------------------------
821  ENDDO ! end of do l=lm,1,-1 loop
822 !------------COMPUTE CAPE AND CINS--------------------------------------
823  lbeg = 1000
824  lend = 0
825  DO j=jsta,jend
826  DO i=ista,iend
827  lbeg = min(ieql(i,j),lbeg)
828  lend = max(lcl(i,j),lend)
829  ENDDO
830  ENDDO
831 !
832 !$omp parallel do private(i,j)
833  DO j=jsta,jend
834  DO i=ista,iend
835  IF(t(i,j,ieql(i,j)) > 255.65) THEN
836  thunder(i,j) = .false.
837  ENDIF
838  ENDDO
839  ENDDO
840 !
841  DO l=lbeg,lend
842 
843 !$omp parallel do private(i,j)
844  DO j=jsta,jend
845  DO i=ista,iend
846  idx(i,j) = 0
847  IF(l >= ieql(i,j).AND.l <= lcl(i,j)) THEN
848  idx(i,j) = 1
849  ENDIF
850  ENDDO
851  ENDDO
852 !
853 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv)
854  DO j=jsta,jend
855  DO i=ista,iend
856  IF(idx(i,j) > 0) THEN
857  presk = pmid(i,j,l)
858  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
859  esatp = min(fpvsnew(tpar(i,j,l)),presk)
860  qsatp = eps*esatp/(presk-esatp*oneps)
861 ! TVP = TPAR(I,J,L)*(1+0.608*QSATP)
862  tvp = tvirtual(tpar(i,j,l),qsatp)
863  thetap = tvp*(h10e5/presk)**capa
864 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
865  tv = tvirtual(t(i,j,l),q(i,j,l))
866  thetaa = tv*(h10e5/presk)**capa
867  IF(thetap < thetaa) THEN
868  cins(i,j) = cins(i,j) + (log(thetap)-log(thetaa))*gdzkl
869  ELSEIF(thetap > thetaa) THEN
870  cape(i,j) = cape(i,j) + (log(thetap)-log(thetaa))*gdzkl
871  IF (thunder(i,j) .AND. t(i,j,l) < 273.15 &
872  .AND. t(i,j,l) > 253.15) THEN
873  cape20(i,j) = cape20(i,j) + (log(thetap)-log(thetaa))*gdzkl
874  ENDIF
875  ENDIF
876  ENDIF
877  ENDDO
878  ENDDO
879  ENDDO
880 !
881 ! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
882 ! LIMIT OF 0.0 ON CINS.
883 !
884 !$omp parallel do private(i,j)
885  DO j=jsta,jend
886  DO i=ista,iend
887  cape(i,j) = max(d00,cape(i,j))
888  cins(i,j) = min(cins(i,j),d00)
889 ! add equillibrium height
890  zeql(i,j) = zint(i,j,ieql(i,j))
891  teql(i,j) = t(i,j,ieql(i,j))
892  IF (cape20(i,j) < 75.) THEN
893  thunder(i,j) = .false.
894  ENDIF
895  IF (thunder(i,j)) THEN
896  thund(i,j) = 1.0
897  ELSE
898  thund(i,j) = 0.0
899  ENDIF
900  ENDDO
901  ENDDO
902 !
903  DEALLOCATE(tpar)
904 !
905  END SUBROUTINE calcape
906 !
907 !-------------------------------------------------------------------------------------
1003  SUBROUTINE calcape2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, &
1004  cape,cins,lfc,esrhl,esrhh, &
1005  dcape,dgld,esp)
1006  use vrbls3d, only: pmid, t, q, zint
1007  use vrbls2d, only: fis,ieql
1008  use gridspec_mod, only: gridtype
1009  use masks, only: lmh
1010  use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, &
1011  oneps, g, tfrz
1012  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
1013  plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, &
1014  itbq, jtbq, rdpq, the0q, stheq, rdtheq
1015  use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, jm, me, jsta_m, jend_m, spval,&
1016  ista_2l, iend_2u, ista, iend, ista_m, iend_m
1017 !
1018 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1019  implicit none
1020 !
1021 ! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980).
1022  real,PARAMETER :: ismthp=2,ismtht=2,ismthq=2
1023 !
1024 ! DECLARE VARIABLES.
1025 !
1026  integer,intent(in) :: itype
1027  real,intent(in) :: dpbnd
1028  integer, dimension(ista:iend,Jsta:jend),intent(in) :: l1d
1029  real, dimension(ista:iend,Jsta:jend),intent(in) :: p1d,t1d
1030 ! real, dimension(ista:iend,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL
1031  real, dimension(ista:iend,jsta:jend),intent(inout) :: q1d,cape,cins
1032  real, dimension(ista:iend,jsta:jend) :: pparc,zeql
1033  real, dimension(ista:iend,jsta:jend),intent(inout) :: lfc,esrhl,esrhh
1034  real, dimension(ista:iend,jsta:jend),intent(inout) :: dcape,dgld,esp
1035  integer, dimension(ista:iend,jsta:jend) ::l12,l17,l3km
1036 !
1037  integer, dimension(ista:iend,jsta:jend) :: iptb, ithtb, parcel, klres, khres, lcl, idx
1038 !
1039  real, dimension(ista:iend,jsta:jend) :: thesp, psp, cape20, qq, pp, thund
1040  integer, dimension(ista:iend,jsta:jend) :: parcel2
1041  real, dimension(ista:iend,jsta:jend) :: thesp2,psp2
1042  real, dimension(ista:iend,jsta:jend) :: cape4,cins4
1043  REAL, ALLOCATABLE :: tpar(:,:,:)
1044  REAL, ALLOCATABLE :: tpar2(:,:,:)
1045 
1046  LOGICAL thunder(ista:iend,jsta:jend), needthun
1047  real psfck,pkl,tbtk,qbtk,apebtk,tthbtk,tthk,apespk,tpspk, &
1048  bqs00k,sqs00k,bqs10k,sqs10k,bqk,sqk,tqk,presk,gdzkl,thetap, &
1049  thetaa,p00k,p10k,p01k,p11k,tthesk,esatp,qsatp,tvp,tv
1050  real presk2, esatp2, qsatp2, tvp2, thetap2, tv2, thetaa2
1051 ! real,external :: fpvsnew
1052  integer i,j,l,knuml,knumh,lbeg,lend,iq, kb,ittbk
1053  integer ie,iw,jn,js,ive(jm),ivw(jm),jvn,jvs
1054  integer istart,istop,jstart,jstop
1055  real, dimension(ista:iend,jsta:jend) :: htsfc
1056 
1057 ! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK
1058 !
1059 !**************************************************************
1060 ! START CALCAPE HERE.
1061 !
1062  ALLOCATE(tpar(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
1063  ALLOCATE(tpar2(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
1064 !
1065 ! COMPUTE CAPE/CINS
1066 !
1067 ! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
1068 ! G * (LN(THETAP) - LN(THETAA)) * DZ
1069 !
1070 ! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
1071 !
1072 ! WHERE:
1073 ! THETAP IS THE PARCEL THETA
1074 ! THETAA IS THE AMBIENT THETA
1075 ! DZ IS THE THICKNESS OF THE LAYER
1076 !
1077 ! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
1078 ! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
1079 !
1080 ! IEQL = EQ LEVEL
1081 ! P_thetaemax - real pressure of theta-e max parcel (Pa)
1082 !
1083 ! INITIALIZE CAPE AND CINS ARRAYS
1084 !
1085 !$omp parallel do
1086  DO j=jsta,jend
1087  DO i=ista,iend
1088  cape(i,j) = d00
1089  cape20(i,j) = d00
1090  cape4(i,j) = d00
1091  cins(i,j) = d00
1092  cins4(i,j) = d00
1093  lcl(i,j) = 0
1094  thesp(i,j) = d00
1095  ieql(i,j) = lm
1096  parcel(i,j) = lm
1097  psp(i,j) = d00
1098  pparc(i,j) = d00
1099  thunder(i,j) = .true.
1100  lfc(i,j) = d00
1101  esrhl(i,j) = d00
1102  esrhh(i,j) = d00
1103  dcape(i,j) = d00
1104  dgld(i,j) = d00
1105  esp(i,j) = d00
1106  thesp2(i,j) = 1000.
1107  psp2(i,j) = d00
1108  parcel2(i,j) = lm
1109  ENDDO
1110  ENDDO
1111 !
1112 !$omp parallel do
1113  DO l=1,lm
1114  DO j=jsta,jend
1115  DO i=ista,iend
1116  tpar(i,j,l) = d00
1117  tpar2(i,j,l) = d00
1118  ENDDO
1119  ENDDO
1120  ENDDO
1121 !
1122 ! FIND SURFACE HEIGHT
1123 !
1124  IF(gridtype == 'E')THEN
1125  jvn = 1
1126  jvs = -1
1127  do j=jsta,jend
1128  ive(j) = mod(j,2)
1129  ivw(j) = ive(j)-1
1130  enddo
1131  istart = ista_m
1132  istop = iend_m
1133  jstart = jsta_m
1134  jstop = jend_m
1135  ELSE IF(gridtype == 'B')THEN
1136  jvn = 1
1137  jvs = 0
1138  do j=jsta,jend
1139  ive(j)=1
1140  ivw(j)=0
1141  enddo
1142  istart = ista_m
1143  istop = iend_m
1144  jstart = jsta_m
1145  jstop = jend_m
1146  ELSE
1147  jvn = 0
1148  jvs = 0
1149  do j=jsta,jend
1150  ive(j) = 0
1151  ivw(j) = 0
1152  enddo
1153  istart = ista
1154  istop = iend
1155  jstart = jsta
1156  jstop = jend
1157  END IF
1158 !!$omp parallel do private(htsfc,ie,iw)
1159  IF(gridtype /= 'A') CALL exch(fis(ista:iend,jsta:jend))
1160  DO j=jstart,jstop
1161  DO i=istart,istop
1162  ie = i+ive(j)
1163  iw = i+ivw(j)
1164  jn = j+jvn
1165  js = j+jvs
1166 !mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+
1167 !mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25
1168 !mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT
1169  IF (gridtype=='B')THEN
1170  htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(ie,jn))
1171  ELSE
1172  htsfc(i,j) = (0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(i,js))
1173  ENDIF
1174  ENDDO
1175  ENDDO
1176 !
1177 ! TYPE 2 CAPE/CINS:
1178 ! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
1179 ! ARE DUMMY ARRAYS.
1180 !
1181  IF (itype == 2) THEN
1182 !$omp parallel do private(i,j)
1183  DO j=jsta,jend
1184  DO i=ista,iend
1185  q1d(i,j) = min(max(h1m12,q1d(i,j)),h99999)
1186  ENDDO
1187  ENDDO
1188  ENDIF
1189 !-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
1190 !-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
1191 !--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
1192 
1193  DO kb=1,lm
1194 !hc IF (ITYPE==2.AND.KB>1) cycle
1195  IF (itype == 1 .OR. (itype == 2 .AND. kb == 1)) THEN
1196 
1197 !$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, &
1198 !$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, &
1199 !$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk)
1200  DO j=jsta,jend
1201  DO i=ista,iend
1202  psfck = pmid(i,j,nint(lmh(i,j)))
1203  pkl = pmid(i,j,kb)
1204 
1205 !hc IF (ITYPE==1.AND.(PKL<PSFCK-DPBND.OR.PKL>PSFCK)) cycle
1206  IF (itype ==2 .OR. &
1207  (itype == 1 .AND. (pkl >= psfck-dpbnd .AND. pkl <= psfck)))THEN
1208  IF (itype == 1) THEN
1209  tbtk = t(i,j,kb)
1210  qbtk = max(0.0, q(i,j,kb))
1211  apebtk = (h10e5/pkl)**capa
1212  ELSE
1213  pkl = p1d(i,j)
1214  tbtk = t1d(i,j)
1215  qbtk = max(0.0, q1d(i,j))
1216  apebtk = (h10e5/pkl)**capa
1217  ENDIF
1218 
1219 !----------Breogan Gomez - 2009-02-06
1220 ! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK
1221 ! and a floating invalid.
1222 
1223 ! if(QBTK < 0) QBTK = 0
1224 
1225 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
1226  tthbtk = tbtk*apebtk
1227  tthk = (tthbtk-thl)*rdth
1228  qq(i,j) = tthk - aint(tthk)
1229  ittbk = int(tthk) + 1
1230 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
1231  IF(ittbk < 1) THEN
1232  ittbk = 1
1233  qq(i,j) = d00
1234  ENDIF
1235  IF(ittbk >= jtb) THEN
1236  ittbk = jtb-1
1237  qq(i,j) = d00
1238  ENDIF
1239 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
1240  bqs00k = qs0(ittbk)
1241  sqs00k = sqs(ittbk)
1242  bqs10k = qs0(ittbk+1)
1243  sqs10k = sqs(ittbk+1)
1244 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
1245  bqk = (bqs10k-bqs00k)*qq(i,j) + bqs00k
1246  sqk = (sqs10k-sqs00k)*qq(i,j) + sqs00k
1247  tqk = (qbtk-bqk)/sqk*rdq
1248  pp(i,j) = tqk-aint(tqk)
1249  iq = int(tqk)+1
1250 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
1251  IF(iq < 1) THEN
1252  iq = 1
1253  pp(i,j) = d00
1254  ENDIF
1255  IF(iq >= itb) THEN
1256  iq = itb-1
1257  pp(i,j) = d00
1258  ENDIF
1259 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
1260  p00k = ptbl(iq ,ittbk )
1261  p10k = ptbl(iq+1,ittbk )
1262  p01k = ptbl(iq ,ittbk+1)
1263  p11k = ptbl(iq+1,ittbk+1)
1264 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
1265  tpspk = p00k + (p10k-p00k)*pp(i,j) + (p01k-p00k)*qq(i,j) &
1266  + (p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
1267 
1268 !!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA
1269  if (tpspk > 1.0e-3) then
1270  apespk = (max(0.,h10e5/ tpspk))**capa
1271  else
1272  apespk = 0.0
1273  endif
1274 
1275  tthesk = tthbtk * exp(elocp*qbtk*apespk/tthbtk)
1276 !--------------CHECK FOR MAXIMUM THETA E--------------------------------
1277  IF(tthesk > thesp(i,j)) THEN
1278  psp(i,j) = tpspk
1279  thesp(i,j) = tthesk
1280  parcel(i,j) = kb
1281  ENDIF
1282 !--------------CHECK FOR MINIMUM THETA E--------------------------------
1283  IF(tthesk < thesp2(i,j)) THEN
1284  psp2(i,j) = tpspk
1285  thesp2(i,j) = tthesk
1286  parcel2(i,j) = kb
1287  ENDIF
1288  END IF
1289  ENDDO ! I loop
1290  ENDDO ! J loop
1291  END IF
1292  ENDDO ! KB loop
1293 
1294 !----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED
1295 !$omp parallel do private(i,j)
1296  DO j=jsta,jend
1297  DO i=ista,iend
1298  pparc(i,j) = pmid(i,j,parcel(i,j))
1299  ENDDO
1300  ENDDO
1301 !
1302 !-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
1303 !-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
1304 !-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
1305  DO l=1,lm
1306 !$omp parallel do private(i,j)
1307  DO j=jsta,jend
1308  DO i=ista,iend
1309  IF (pmid(i,j,l) < psp(i,j)) lcl(i,j) = l+1
1310  ENDDO
1311  ENDDO
1312  ENDDO
1313 !$omp parallel do private(i,j)
1314  DO j=jsta,jend
1315  DO i=ista,iend
1316  IF (lcl(i,j) > nint(lmh(i,j))) lcl(i,j) = nint(lmh(i,j))
1317  IF (itype > 2) THEN
1318  IF (t(i,j,lcl(i,j)) < 263.15) THEN
1319  thunder(i,j) = .false.
1320  ENDIF
1321  ENDIF
1322  ENDDO
1323  ENDDO
1324 !-----------------------------------------------------------------------
1325 !---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
1326 !-----------------------------------------------------------------------
1327  DO l=lm,1,-1
1328 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
1329  knuml = 0
1330  knumh = 0
1331  DO j=jsta,jend
1332  DO i=ista,iend
1333  klres(i,j) = 0
1334  khres(i,j) = 0
1335  IF(l <= lcl(i,j)) THEN
1336  IF(pmid(i,j,l) < plq)THEN
1337  knuml = knuml + 1
1338  klres(i,j) = 1
1339  ELSE
1340  knumh = knumh + 1
1341  khres(i,j) = 1
1342  ENDIF
1343  ENDIF
1344  ENDDO
1345  ENDDO
1346 !***
1347 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
1348 !**
1349  IF(knuml > 0) THEN
1350  CALL ttblex(tpar(ista_2l,jsta_2l,l),ttbl,itb,jtb,klres &
1351  , pmid(ista_2l,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
1352  , rdthe,thesp,iptb,ithtb)
1353  ENDIF
1354 !***
1355 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
1356 !**
1357  IF(knumh > 0) THEN
1358  CALL ttblex(tpar(ista_2l,jsta_2l,l),ttblq,itbq,jtbq,khres &
1359  , pmid(ista_2l,jsta_2l,l),plq,qq,pp,rdpq &
1360  ,the0q,stheq,rdtheq,thesp,iptb,ithtb)
1361  ENDIF
1362 
1363 !------------SEARCH FOR EQ LEVEL----------------------------------------
1364 !$omp parallel do private(i,j)
1365  DO j=jsta,jend
1366  DO i=ista,iend
1367  IF(khres(i,j) > 0) THEN
1368  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
1369  ENDIF
1370  ENDDO
1371  ENDDO
1372 !
1373 !$omp parallel do private(i,j)
1374  DO j=jsta,jend
1375  DO i=ista,iend
1376  IF(klres(i,j) > 0) THEN
1377  IF(tpar(i,j,l) > t(i,j,l)) ieql(i,j) = l
1378  ENDIF
1379  ENDDO
1380  ENDDO
1381 !-----------------------------------------------------------------------
1382  ENDDO ! end of do l=lm,1,-1 loop
1383 !------------COMPUTE CAPE AND CINS--------------------------------------
1384  lbeg = 1000
1385  lend = 0
1386  DO j=jsta,jend
1387  DO i=ista,iend
1388  lbeg = min(ieql(i,j),lbeg)
1389  lend = max(lcl(i,j),lend)
1390  ENDDO
1391  ENDDO
1392 !
1393 !$omp parallel do private(i,j)
1394  DO j=jsta,jend
1395  DO i=ista,iend
1396  IF(t(i,j,ieql(i,j)) > 255.65) THEN
1397  thunder(i,j) = .false.
1398  ENDIF
1399  ENDDO
1400  ENDDO
1401 !
1402 !reverse L order from bottom up for ESRH calculation
1403 !
1404  esrhh = lcl
1405  esrhl = lcl
1406 ! DO L=LBEG,LEND
1407  DO l=lend,lbeg,-1
1408 
1409 !$omp parallel do private(i,j)
1410  DO j=jsta,jend
1411  DO i=ista,iend
1412  idx(i,j) = 0
1413  IF(l >= ieql(i,j).AND.l <= lcl(i,j)) THEN
1414  idx(i,j) = 1
1415  ENDIF
1416  ENDDO
1417  ENDDO
1418 !
1419 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv,&
1420 !$omp & presk2,esatp2,qsatp2,tvp2,thetap2,tv2,thetaa2)
1421  DO j=jsta,jend
1422  DO i=ista,iend
1423  IF(idx(i,j) > 0) THEN
1424  presk = pmid(i,j,l)
1425  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
1426  esatp = min(fpvsnew(tpar(i,j,l)),presk)
1427  qsatp = eps*esatp/(presk-esatp*oneps)
1428 ! TVP = TPAR(I,J,L)*(1+0.608*QSATP)
1429  tvp = tvirtual(tpar(i,j,l),qsatp)
1430  thetap = tvp*(h10e5/presk)**capa
1431 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
1432  tv = tvirtual(t(i,j,l),q(i,j,l))
1433  thetaa = tv*(h10e5/presk)**capa
1434  IF(thetap < thetaa) THEN
1435  cins4(i,j) = cins4(i,j) + (log(thetap)-log(thetaa))*gdzkl
1436  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1437  cins(i,j) = cins(i,j) + (log(thetap)-log(thetaa))*gdzkl
1438  ENDIF
1439  ELSEIF(thetap > thetaa) THEN
1440  cape4(i,j) = cape4(i,j) + (log(thetap)-log(thetaa))*gdzkl
1441  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1442  cape(i,j) = cape(i,j) + (log(thetap)-log(thetaa))*gdzkl
1443  ENDIF
1444  IF (thunder(i,j) .AND. t(i,j,l) < 273.15 &
1445  .AND. t(i,j,l) > 253.15) THEN
1446  cape20(i,j) = cape20(i,j) + (log(thetap)-log(thetaa))*gdzkl
1447  ENDIF
1448  ENDIF
1449 
1450 ! LFC
1451  IF (itype /= 1) THEN
1452  presk2 = pmid(i,j,l+1)
1453  esatp2 = min(fpvsnew(tpar(i,j,l+1)),presk2)
1454  qsatp2 = eps*esatp2/(presk2-esatp2*oneps)
1455 ! TVP2 = TPAR(I,J,L+1)*(1+0.608*QSATP2)
1456  tvp2 = tvirtual(tpar(i,j,l+1),qsatp2)
1457  thetap2 = tvp2*(h10e5/presk2)**capa
1458 ! TV2 = T(I,J,L+1)*(1+0.608*Q(I,J,L+1))
1459  tv2 = tvirtual(t(i,j,l+1),q(i,j,l+1))
1460  thetaa2 = tv2*(h10e5/presk2)**capa
1461  IF(thetap >= thetaa .AND. thetap2 <= thetaa2) THEN
1462  IF(lfc(i,j) == d00)THEN
1463  lfc(i,j) = zint(i,j,l)
1464  ENDIF
1465  ENDIF
1466  ENDIF
1467 !
1468 ! ESRH/CAPE threshold check
1469  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) THEN
1470  IF(cape4(i,j) >= 100. .AND. cins4(i,j) >= -250.) THEN
1471  IF(esrhl(i,j) == lcl(i,j)) esrhl(i,j)=l
1472  ENDIF
1473  esrhh(i,j)=l
1474  ENDIF
1475 
1476  ENDIF !(IDX(I,J) > 0)
1477  ENDDO
1478  ENDDO
1479  ENDDO
1480 
1481 !$omp parallel do private(i,j)
1482  DO j=jsta,jend
1483  DO i=ista,iend
1484  IF(esrhh(i,j) > esrhl(i,j)) esrhh(i,j)=ieql(i,j)
1485  ENDDO
1486  ENDDO
1487 !
1488 ! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
1489 ! LIMIT OF 0.0 ON CINS.
1490 ! ENFORCE LFC ABOVE LCL AND BELOW EL
1491 !
1492 !$omp parallel do private(i,j)
1493  DO j=jsta,jend
1494  DO i=ista,iend
1495  cape(i,j) = max(d00,cape(i,j))
1496  cins(i,j) = min(cins(i,j),d00)
1497 ! equillibrium height
1498  zeql(i,j) = zint(i,j,ieql(i,j))
1499  lfc(i,j) = min(lfc(i,j),zint(i,j,ieql(i,j)))
1500  lfc(i,j) = max(zint(i,j, lcl(i,j)),lfc(i,j))
1501  IF (cape20(i,j) < 75.) THEN
1502  thunder(i,j) = .false.
1503  ENDIF
1504  IF (thunder(i,j)) THEN
1505  thund(i,j) = 1.0
1506  ELSE
1507  thund(i,j) = 0.0
1508  ENDIF
1509  ENDDO
1510  ENDDO
1511 !------------COMPUTE DCAPE--------------------------------------
1512 !-----------------------------------------------------------------------
1513 !---------FIND TEMP OF PARCEL DESCENDED ALONG MOIST ADIABAT (TPAR)---------
1514 !-----------------------------------------------------------------------
1515  IF (itype == 1) THEN
1516 
1517  DO l=lm,1,-1
1518 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
1519  knuml = 0
1520  knumh = 0
1521  DO j=jsta,jend
1522  DO i=ista,iend
1523  klres(i,j) = 0
1524  khres(i,j) = 0
1525  psfck = pmid(i,j,nint(lmh(i,j)))
1526  pkl = pmid(i,j,l)
1527  IF(pkl >= psfck-dpbnd) THEN
1528  IF(pmid(i,j,l) < plq)THEN
1529  knuml = knuml + 1
1530  klres(i,j) = 1
1531  ELSE
1532  knumh = knumh + 1
1533  khres(i,j) = 1
1534  ENDIF
1535  ENDIF
1536  ENDDO
1537  ENDDO
1538 !***
1539 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
1540 !**
1541  IF(knuml > 0) THEN
1542  CALL ttblex(tpar2(ista_2l,jsta_2l,l),ttbl,itb,jtb,klres &
1543  , pmid(ista_2l,jsta_2l,l),pl,qq,pp,rdp,the0,sthe &
1544  , rdthe,thesp2,iptb,ithtb)
1545  ENDIF
1546 !***
1547 !*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
1548 !**
1549  IF(knumh > 0) THEN
1550  CALL ttblex(tpar2(ista_2l,jsta_2l,l),ttblq,itbq,jtbq,khres &
1551  , pmid(ista_2l,jsta_2l,l),plq,qq,pp,rdpq &
1552  , the0q,stheq,rdtheq,thesp2,iptb,ithtb)
1553  ENDIF
1554  ENDDO ! end of do l=lm,1,-1 loop
1555 
1556  lbeg = 1
1557  lend = lm
1558 
1559  DO l=lbeg,lend
1560 !$omp parallel do private(i,j)
1561  DO j=jsta,jend
1562  DO i=ista,iend
1563  idx(i,j) = 0
1564  IF(l >= parcel2(i,j).AND.l < nint(lmh(i,j))) THEN
1565  idx(i,j) = 1
1566  ENDIF
1567  ENDDO
1568  ENDDO
1569 !
1570 !$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv)
1571  DO j=jsta,jend
1572  DO i=ista,iend
1573  IF(idx(i,j) > 0) THEN
1574  presk = pmid(i,j,l)
1575  gdzkl = (zint(i,j,l)-zint(i,j,l+1)) * g
1576  esatp = min(fpvsnew(tpar2(i,j,l)),presk)
1577  qsatp = eps*esatp/(presk-esatp*oneps)
1578 ! TVP = TPAR2(I,J,L)*(1+0.608*QSATP)
1579  tvp = tvirtual(tpar2(i,j,l),qsatp)
1580  thetap = tvp*(h10e5/presk)**capa
1581 ! TV = T(I,J,L)*(1+0.608*Q(I,J,L))
1582  tv = tvirtual(t(i,j,l),q(i,j,l))
1583  thetaa = tv*(h10e5/presk)**capa
1584  !IF(THETAP < THETAA) THEN
1585  dcape(i,j) = dcape(i,j) + (log(thetap)-log(thetaa))*gdzkl
1586  !ENDIF
1587  ENDIF
1588  ENDDO
1589  ENDDO
1590  ENDDO
1591 
1592 !$omp parallel do private(i,j)
1593  DO j=jsta,jend
1594  DO i=ista,iend
1595  dcape(i,j) = min(d00,dcape(i,j))
1596  ENDDO
1597  ENDDO
1598 
1599  ENDIF !ITYPE=1 FOR DCAPE
1600 
1601 !
1602 ! Dendritic Growth Layer depth
1603 ! the layer with temperatures from -12 to -17 C in meters
1604 !
1605  l12=lm
1606  l17=lm
1607  DO l=lm,1,-1
1608 !$omp parallel do private(i,j)
1609  DO j=jsta,jend
1610  DO i=ista,iend
1611  IF(t(i,j,l) <= tfrz-12. .AND. l12(i,j)==lm) l12(i,j)=l
1612  IF(t(i,j,l) <= tfrz-17. .AND. l17(i,j)==lm) l17(i,j)=l
1613  ENDDO
1614  ENDDO
1615  ENDDO
1616 !$omp parallel do private(i,j)
1617  DO j=jsta,jend
1618  DO i=ista,iend
1619  IF(l12(i,j)/=lm .AND. l17(i,j)/=lm) THEN
1620  dgld(i,j)=zint(i,j,l17(i,j))-zint(i,j,l12(i,j))
1621  dgld(i,j)=max(dgld(i,j),0.)
1622  ENDIF
1623  ENDDO
1624  ENDDO
1625 !
1626 ! Enhanced Stretching Potential
1627 ! ESP = (0-3 km MLCAPE / 50 J kg-1) * ((0-3 km lapse rate - 7.0) / 1.0 C (km-1)
1628 ! https://www.spc.noaa.gov/exper/soundings/help/params4.html
1629 !
1630  l3km=lm
1631  DO l=lm,1,-1
1632 !$omp parallel do private(i,j)
1633  DO j=jsta,jend
1634  DO i=ista,iend
1635  IF(zint(i,j,l)-htsfc(i,j) <= 3000.) l3km(i,j)=l
1636  ENDDO
1637  ENDDO
1638  ENDDO
1639 !$omp parallel do private(i,j)
1640  DO j=jsta,jend
1641  DO i=ista,iend
1642  esp(i,j) = (cape(i,j) / 50.) * (t(i,j,lm) - t(i,j,l3km(i,j)) - 7.0)
1643  IF((t(i,j,lm) - t(i,j,l3km(i,j))) < 7.0) esp(i,j) = 0.
1644 ! IF(CAPE(I,J) < 250.) ESP(I,J) = 0.
1645  ENDDO
1646  ENDDO
1647 !
1648  DEALLOCATE(tpar)
1649  DEALLOCATE(tpar2)
1650 !
1651  END SUBROUTINE calcape2
1652 !
1653 !-------------------------------------------------------------------------------------
1654 !
1655 !
1656 !-------------------------------------------------------------------------------------
1657 !
1658  elemental function tvirtual(T,Q)
1659 !
1660 ! COMPUTE VIRTUAL TEMPERATURE
1661 !
1662  IMPLICIT NONE
1663  REAL tvirtual
1664  REAL, INTENT(IN) :: t, q
1665 
1666  tvirtual = t*(1+0.608*q)
1667 
1668  end function tvirtual
1669 !
1670 !-------------------------------------------------------------------------------------
1671 !
1697  SUBROUTINE calvor(UWND,VWND,ABSV)
1698 
1699 !
1700 !
1701  use vrbls2d, only: f
1702  use masks, only: gdlat, gdlon, dx, dy
1703  use params_mod, only: d00, dtr, small, erad
1704  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
1705  jsta, jend, im, jm, jsta_m, jend_m, gdsdegr,&
1706  ista, iend, ista_m, iend_m, ista_2l, iend_2u, me, num_procs
1707  use gridspec_mod, only: gridtype, dyval
1708  use upp_math, only: dvdxdudy, ddvdx, ddudy, uuavg
1709 
1710  implicit none
1711 !
1712 ! DECLARE VARIABLES.
1713 !
1714  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in) :: uwnd, vwnd
1715  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(inout) :: absv
1716  REAL, dimension(IM,2) :: glatpoles, coslpoles, upoles, avpoles
1717  REAL, dimension(IM,JSTA:JEND) :: cosltemp, avtemp
1718 !
1719  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
1720  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
1721 !
1722  integer, parameter :: npass2=2, npass3=3
1723  integer i,j,ip1,im1,ii,iir,iil,jj,jmt2,imb2, npass, nn, jtem
1724  real r2dx,r2dy,dvdx,dudy,uavg,tph1,tphi, tx1(im+2), tx2(im+2)
1725 !
1726 !***************************************************************************
1727 ! START CALVOR HERE.
1728 !
1729 ! LOOP TO COMPUTE ABSOLUTE VORTICITY FROM WINDS.
1730 !
1731  IF(modelname == 'RAPR') then
1732 !$omp parallel do private(i,j)
1733  DO j=jsta_2l,jend_2u
1734  DO i=ista_2l,iend_2u
1735  absv(i,j) = d00
1736  ENDDO
1737  ENDDO
1738  else
1739 !$omp parallel do private(i,j)
1740  DO j=jsta_2l,jend_2u
1741  DO i=ista_2l,iend_2u
1742  absv(i,j) = spval
1743  ENDDO
1744  ENDDO
1745  endif
1746 
1747 ! print*,'dyval in CALVOR= ',DYVAL
1748 
1749  CALL exch(uwnd)
1750  CALL exch(vwnd)
1751 !
1752  IF (modelname == 'GFS' .or. global) THEN
1753  CALL exch(gdlat(ista_2l,jsta_2l))
1754  CALL exch(gdlon(ista_2l,jsta_2l))
1755 
1756  allocate (wrk1(ista:iend,jsta:jend), wrk2(ista:iend,jsta:jend), &
1757  & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u))
1758  allocate(iw(im),ie(im))
1759 
1760  imb2 = im/2
1761 !$omp parallel do private(i)
1762  do i=ista,iend
1763  ie(i) = i+1
1764  iw(i) = i-1
1765  enddo
1766 ! iw(1) = im
1767 ! ie(im) = 1
1768 
1769 ! if(1>=jsta .and. 1<=jend)then
1770 ! if(cos(gdlat(1,1)*dtr)<small)poleflag=.T.
1771 ! end if
1772 ! call mpi_bcast(poleflag,1,MPI_LOGICAL,0,mpi_comm_comp,iret)
1773 
1774 !$omp parallel do private(i,j,ip1,im1)
1775  DO j=jsta,jend
1776  do i=ista,iend
1777  ip1 = ie(i)
1778  im1 = iw(i)
1779  cosl(i,j) = cos(gdlat(i,j)*dtr)
1780  IF(cosl(i,j) >= small) then
1781  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
1782  else
1783  wrk1(i,j) = 0.
1784  end if
1785  if(i == im .or. i == 1) then
1786  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
1787  else
1788  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
1789  end if
1790  enddo
1791  enddo
1792 ! CALL EXCH(cosl(1,JSTA_2L))
1793  CALL exch(cosl)
1794 
1795  call fullpole( cosl(ista_2l:iend_2u,jsta_2l:jend_2u),coslpoles)
1796  call fullpole(gdlat(ista_2l:iend_2u,jsta_2l:jend_2u),glatpoles)
1797 
1798  if(me==0 ) print*,'CALVOR ',me,glatpoles(ista,1),glatpoles(ista,2)
1799  if(me==num_procs-1) print*,'CALVOR ',me,glatpoles(ista,1),glatpoles(ista,2)
1800 
1801 !$omp parallel do private(i,j,ii)
1802  DO j=jsta,jend
1803  if (j == 1) then
1804  if(gdlat(ista,j) > 0.) then ! count from north to south
1805  do i=ista,iend
1806  ii = i + imb2
1807  if (ii > im) ii = ii - im
1808  ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J+1)-GDLAT(II,J))*DTR) !1/dphi
1809  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-glatpoles(ii,1))*dtr) !1/dphi
1810  enddo
1811  else ! count from south to north
1812  do i=ista,iend
1813  ii = i + imb2
1814  if (ii > im) ii = ii - im
1815  ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J+1)+GDLAT(II,J))*DTR) !1/dphi
1816  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+glatpoles(ii,1))*dtr) !1/dphi
1817 !
1818  enddo
1819  end if
1820  elseif (j == jm) then
1821  if(gdlat(ista,j) < 0.) then ! count from north to south
1822  do i=ista,iend
1823  ii = i + imb2
1824  if (ii > im) ii = ii - im
1825  ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J-1)+GDLAT(II,J))*DTR)
1826  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+glatpoles(ii,2))*dtr)
1827  enddo
1828  else ! count from south to north
1829  do i=ista,iend
1830  ii = i + imb2
1831  if (ii > im) ii = ii - im
1832  ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J-1)-GDLAT(II,J))*DTR)
1833  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-glatpoles(ii,2))*dtr)
1834  enddo
1835  end if
1836  else
1837  do i=ista,iend
1838  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
1839  enddo
1840  endif
1841  enddo
1842 
1843  npass = 0
1844 
1845  jtem = jm / 18 + 1
1846 
1847  call fullpole(uwnd(ista_2l:iend_2u,jsta_2l:jend_2u),upoles)
1848 
1849 !$omp parallel do private(i,j,ip1,im1,ii,jj,tx1,tx2)
1850  DO j=jsta,jend
1851 ! npass = npass2
1852 ! if (j > jm-jtem+1 .or. j < jtem) npass = npass3
1853  IF(j == 1) then ! Near North or South pole
1854  if(gdlat(ista,j) > 0.) then ! count from north to south
1855  IF(cosl(ista,j) >= small) THEN !not a pole point
1856  DO i=ista,iend
1857  ip1 = ie(i)
1858  im1 = iw(i)
1859  ii = i + imb2
1860  if (ii > im) ii = ii - im
1861  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
1862 ! UWND(II,J)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
1863  upoles(ii,1)==spval .or. uwnd(i,j+1)==spval) cycle
1864  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
1865 ! & + (UWND(II,J)*COSL(II,J) &
1866  & + (upoles(ii,1)*coslpoles(ii,1) &
1867  & + uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
1868  & + f(i,j)
1869  enddo
1870  ELSE !pole point, compute at j=2
1871  jj = 2
1872  DO i=ista,iend
1873  ip1 = ie(i)
1874  im1 = iw(i)
1875  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
1876  uwnd(i,j)==spval .or. uwnd(i,jj+1)==spval) cycle
1877  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
1878  & - (uwnd(i,j)*cosl(i,j) &
1879  - uwnd(i,jj+1)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
1880  & + f(i,jj)
1881  enddo
1882  ENDIF
1883  else
1884  IF(cosl(ista,j) >= small) THEN !not a pole point
1885  DO i=ista,iend
1886  ip1 = ie(i)
1887  im1 = iw(i)
1888  ii = i + imb2
1889  if (ii > im) ii = ii - im
1890  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
1891 ! UWND(II,J)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
1892  upoles(ii,1)==spval .or. uwnd(i,j+1)==spval) cycle
1893  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
1894 ! & - (UWND(II,J)*COSL(II,J) &
1895  & - (upoles(ii,1)*coslpoles(ii,1) &
1896  & + uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
1897  & + f(i,j)
1898  enddo
1899  ELSE !pole point, compute at j=2
1900  jj = 2
1901  DO i=ista,iend
1902  ip1 = ie(i)
1903  im1 = iw(i)
1904  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
1905  uwnd(i,j)==spval .or. uwnd(i,jj+1)==spval) cycle
1906  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
1907  & + (uwnd(i,j)*cosl(i,j) &
1908  - uwnd(i,jj+1)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
1909  & + f(i,jj)
1910  enddo
1911  ENDIF
1912  endif
1913  ELSE IF(j == jm) THEN ! Near North or South Pole
1914  if(gdlat(ista,j) < 0.) then ! count from north to south
1915  IF(cosl(ista,j) >= small) THEN !not a pole point
1916  DO i=ista,iend
1917  ip1 = ie(i)
1918  im1 = iw(i)
1919  ii = i + imb2
1920  if (ii > im) ii = ii - im
1921  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
1922 ! UWND(I,J-1)==SPVAL .or. UWND(II,J)==SPVAL) cycle
1923  uwnd(i,j-1)==spval .or. upoles(ii,2)==spval) cycle
1924  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
1925  & - (uwnd(i,j-1)*cosl(i,j-1) &
1926 ! & + UWND(II,J)*COSL(II,J))*wrk3(i,j)) * wrk1(i,j) &
1927  & + upoles(ii,2)*coslpoles(ii,2))*wrk3(i,j)) * wrk1(i,j) &
1928  & + f(i,j)
1929  enddo
1930  ELSE !pole point,compute at jm-1
1931  jj = jm-1
1932  DO i=ista,iend
1933  ip1 = ie(i)
1934  im1 = iw(i)
1935  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
1936  uwnd(i,jj-1)==spval .or. uwnd(i,j)==spval) cycle
1937  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
1938  & - (uwnd(i,jj-1)*cosl(i,jj-1) &
1939  & - uwnd(i,j)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj) &
1940  & + f(i,jj)
1941  enddo
1942  ENDIF
1943  else
1944  IF(cosl(ista,j) >= small) THEN !not a pole point
1945  DO i=ista,iend
1946  ip1 = ie(i)
1947  im1 = iw(i)
1948  ii = i + imb2
1949  if (ii > im) ii = ii - im
1950  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
1951 ! UWND(I,J-1)==SPVAL .or. UWND(II,J)==SPVAL) cycle
1952  uwnd(i,j-1)==spval .or. upoles(ii,2)==spval) cycle
1953  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
1954  & + (uwnd(i,j-1)*cosl(i,j-1) &
1955 ! & + UWND(II,J)*COSL(II,J))*wrk3(i,j)) * wrk1(i,j) &
1956  & + upoles(ii,2)*coslpoles(ii,2))*wrk3(i,j)) * wrk1(i,j) &
1957  & + f(i,j)
1958  enddo
1959  ELSE !pole point,compute at jm-1
1960  jj = jm-1
1961  DO i=ista,iend
1962  ip1 = ie(i)
1963  im1 = iw(i)
1964  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
1965  uwnd(i,jj-1)==spval .or. uwnd(i,j)==spval) cycle
1966  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
1967  & + (uwnd(i,jj-1)*cosl(i,jj-1) &
1968  & - uwnd(i,j)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj) &
1969  & + f(i,jj)
1970  enddo
1971  ENDIF
1972  endif
1973  ELSE
1974  DO i=ista,iend
1975  ip1 = ie(i)
1976  im1 = iw(i)
1977  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
1978  uwnd(i,j-1)==spval .or. uwnd(i,j+1)==spval) cycle
1979  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
1980  & - (uwnd(i,j-1)*cosl(i,j-1) &
1981  - uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
1982  + f(i,j)
1983  ENDDO
1984  END IF
1985 ! if(ABSV(I,J)>1.0)print*,'Debug CALVOR',i,j,VWND(ip1,J),VWND(im1,J), &
1986 ! wrk2(i,j),UWND(I,J-1),COSL(I,J-1),UWND(I,J+1),COSL(I,J+1),wrk3(i,j),cosl(i,j),F(I,J),ABSV(I,J)
1987  if (npass > 0) then
1988  do i=ista,iend
1989  tx1(i) = absv(i,j)
1990  enddo
1991  do nn=1,npass
1992  do i=ista,iend
1993  tx2(i+1) = tx1(i)
1994  enddo
1995  tx2(1) = tx2(im+1)
1996  tx2(im+2) = tx2(2)
1997  do i=2,im+1
1998  tx1(i-1) = 0.25 * (tx2(i-1) + tx2(i+1)) + 0.5*tx2(i)
1999  enddo
2000  enddo
2001  do i=ista,iend
2002  absv(i,j) = tx1(i)
2003  enddo
2004  endif
2005  END DO ! end of J loop
2006 
2007 ! deallocate (wrk1, wrk2, wrk3, cosl)
2008 ! GFS use lon avg as one scaler value for pole point
2009 
2010  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1,jsta),SPVAL,ABSV(1,jsta))
2011 
2012  call exch(absv(ista_2l:iend_2u,jsta_2l:jend_2u))
2013  call fullpole(absv(ista_2l:iend_2u,jsta_2l:jend_2u),avpoles)
2014 
2015  cosltemp=spval
2016  if(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
2017  if(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
2018  avtemp=spval
2019  if(jsta== 1) avtemp(1:im, 1)=avpoles(1:im,1)
2020  if(jend==jm) avtemp(1:im,jm)=avpoles(1:im,2)
2021 
2022  call poleavg(im,jm,jsta,jend,small,cosltemp(1,jsta),spval,avtemp(1,jsta))
2023 
2024  if(jsta== 1) absv(ista:iend, 1)=avtemp(ista:iend, 1)
2025  if(jend==jm) absv(ista:iend,jm)=avtemp(ista:iend,jm)
2026 
2027  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
2028 
2029  ELSE !(MODELNAME == 'GFS' .or. global)
2030 
2031  IF (gridtype == 'B')THEN
2032  CALL exch(vwnd)
2033  CALL exch(uwnd)
2034  ENDIF
2035 
2036  CALL dvdxdudy(uwnd,vwnd)
2037 
2038  IF(gridtype == 'A')THEN
2039 !$omp parallel do private(i,j,jmt2,tphi,r2dx,r2dy,dvdx,dudy,uavg)
2040  DO j=jsta_m,jend_m
2041  jmt2 = jm/2+1
2042  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
2043  DO i=ista_m,iend_m
2044  IF(vwnd(i+1,j)<spval.AND.vwnd(i-1,j)<spval.AND. &
2045  & uwnd(i,j+1)<spval.AND.uwnd(i,j-1)<spval) THEN
2046  dvdx = ddvdx(i,j)
2047  dudy = ddudy(i,j)
2048  uavg = uuavg(i,j)
2049 ! is there a (f+tan(phi)/erad)*u term?
2050  IF(modelname == 'RAPR') then
2051  absv(i,j) = dvdx - dudy + f(i,j) ! for run RAP over north pole
2052  else
2053  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(gdlat(i,j)*dtr)/erad ! not sure about this???
2054  endif
2055  END IF
2056  END DO
2057  END DO
2058 
2059  ELSE IF (gridtype == 'E')THEN
2060  allocate(ihw(jsta_2l:jend_2u), ihe(jsta_2l:jend_2u))
2061 !$omp parallel do private(j)
2062  DO j=jsta_2l,jend_2u
2063  ihw(j) = -mod(j,2)
2064  ihe(j) = ihw(j)+1
2065  ENDDO
2066 !$omp parallel do private(i,j,jmt2,tphi,r2dx,r2dy,dvdx,dudy,uavg)
2067  DO j=jsta_m,jend_m
2068  jmt2 = jm/2+1
2069  tphi = (j-jmt2)*(dyval/1000.)*dtr
2070  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
2071  DO i=ista_m,iend_m
2072  IF(vwnd(i+ihe(j),j) < spval.AND.vwnd(i+ihw(j),j) < spval .AND. &
2073  & uwnd(i,j+1) < spval .AND.uwnd(i,j-1) < spval) THEN
2074  dvdx = ddvdx(i,j)
2075  dudy = ddudy(i,j)
2076  uavg = uuavg(i,j)
2077 ! is there a (f+tan(phi)/erad)*u term?
2078  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(tphi)/erad
2079  END IF
2080  END DO
2081  END DO
2082  deallocate(ihw, ihe)
2083  ELSE IF (gridtype == 'B')THEN
2084 ! CALL EXCH(VWND) !done before dvdxdudy() Jesse 20200520
2085  DO j=jsta_m,jend_m
2086  jmt2 = jm/2+1
2087  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
2088  DO i=ista_m,iend_m
2089  if(vwnd(i, j)==spval .or. vwnd(i, j-1)==spval .or. &
2090  vwnd(i-1,j)==spval .or. vwnd(i-1,j-1)==spval .or. &
2091  uwnd(i, j)==spval .or. uwnd(i-1,j)==spval .or. &
2092  uwnd(i,j-1)==spval .or. uwnd(i-1,j-1)==spval) cycle
2093  dvdx = ddvdx(i,j)
2094  dudy = ddudy(i,j)
2095  uavg = uuavg(i,j)
2096 ! is there a (f+tan(phi)/erad)*u term?
2097  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(tphi)/erad
2098  END DO
2099  END DO
2100  END IF
2101  END IF
2102 !
2103 ! END OF ROUTINE.
2104 !
2105  RETURN
2106  END
2107 
2124  SUBROUTINE caldiv(UWND,VWND,DIV)
2125  use masks, only: gdlat, gdlon
2126  use params_mod, only: d00, dtr, small, erad
2127  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
2128  jsta, jend, im, jm, jsta_m, jend_m, lm, &
2129  ista, iend, ista_m, iend_m, ista_2l, iend_2u
2130  use gridspec_mod, only: gridtype
2131 
2132  implicit none
2133 !
2134 ! DECLARE VARIABLES.
2135 !
2136  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u,lm), intent(in) :: uwnd,vwnd
2137  REAL, dimension(ista:iend,jsta:jend,lm), intent(inout) :: div
2138  REAL, dimension(IM,2) :: glatpoles, coslpoles, upoles, vpoles, divpoles
2139  REAL, dimension(IM,JSTA:JEND) :: cosltemp, divtemp
2140 !
2141  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
2142  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
2143 !
2144  real :: dnpole, dspole, tem
2145  integer i,j,ip1,im1,ii,iir,iil,jj,imb2, l
2146 !
2147 !***************************************************************************
2148 ! START CALDIV HERE.
2149 !
2150 ! LOOP TO COMPUTE DIVERGENCE FROM WINDS.
2151 !
2152  CALL exch(gdlat(ista_2l,jsta_2l))
2153  CALL exch(gdlon(ista_2l,jsta_2l))
2154 
2155  allocate (wrk1(ista:iend,jsta:jend), wrk2(ista:iend,jsta:jend), &
2156  & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u))
2157  allocate(iw(im),ie(im))
2158 
2159  imb2 = im/2
2160 !$omp parallel do private(i)
2161  do i=ista,iend
2162  ie(i) = i+1
2163  iw(i) = i-1
2164  enddo
2165 ! iw(1) = im
2166 ! ie(im) = 1
2167 
2168 
2169 !$omp parallel do private(i,j,ip1,im1)
2170  DO j=jsta,jend
2171  do i=ista,iend
2172  ip1 = ie(i)
2173  im1 = iw(i)
2174  cosl(i,j) = cos(gdlat(i,j)*dtr)
2175  IF(cosl(i,j) >= small) then
2176  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
2177  else
2178  wrk1(i,j) = 0.
2179  end if
2180  if(i == im .or. i == 1) then
2181  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
2182  else
2183  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
2184  end if
2185  enddo
2186  ENDDO
2187 
2188  CALL exch(cosl)
2189  CALL fullpole(cosl,coslpoles)
2190  CALL fullpole(gdlat(ista_2l:iend_2u,jsta_2l:jend_2u),glatpoles)
2191 
2192 !$omp parallel do private(i,j,ii)
2193  DO j=jsta,jend
2194  if (j == 1) then
2195  if(gdlat(ista,j) > 0.) then ! count from north to south
2196  do i=ista,iend
2197  ii = i + imb2
2198  if (ii > im) ii = ii - im
2199  ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J+1)-GDLAT(II,J))*DTR) !1/dphi
2200  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-glatpoles(ii,1))*dtr) !1/dphi
2201  enddo
2202  else ! count from south to north
2203  do i=ista,iend
2204  ii = i + imb2
2205  if (ii > im) ii = ii - im
2206  ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J+1)+GDLAT(II,J))*DTR) !1/dphi
2207  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+glatpoles(ii,1))*dtr) !1/dphi
2208  enddo
2209  end if
2210  elseif (j == jm) then
2211  if(gdlat(ista,j) < 0.) then ! count from north to south
2212  do i=ista,iend
2213  ii = i + imb2
2214  if (ii > im) ii = ii - im
2215  ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J-1)+GDLAT(II,J))*DTR)
2216  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+glatpoles(ii,2))*dtr)
2217  enddo
2218  else ! count from south to north
2219  do i=ista,iend
2220  ii = i + imb2
2221  if (ii > im) ii = ii - im
2222  ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J-1)-GDLAT(II,J))*DTR)
2223  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-glatpoles(ii,2))*dtr)
2224  enddo
2225  end if
2226  else
2227  do i=ista,iend
2228  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
2229  enddo
2230  endif
2231  enddo
2232 
2233  do l=1,lm
2234 !$omp parallel do private(i,j)
2235  DO j=jsta,jend
2236  DO i=ista,iend
2237  div(i,j,l) = spval
2238  ENDDO
2239  ENDDO
2240 
2241  CALL exch(vwnd(ista_2l,jsta_2l,l))
2242  CALL exch(uwnd(ista_2l,jsta_2l,l))
2243 
2244  CALL fullpole(vwnd(ista_2l:iend_2u,jsta_2l:jend_2u,l),vpoles)
2245  CALL fullpole(uwnd(ista_2l:iend_2u,jsta_2l:jend_2u,l),upoles)
2246 
2247 !$omp parallel do private(i,j,ip1,im1,ii,jj)
2248  DO j=jsta,jend
2249  IF(j == 1) then ! Near North pole
2250  if(gdlat(ista,j) > 0.) then ! count from north to south
2251  IF(cosl(ista,j) >= small) THEN !not a pole point
2252  DO i=ista,iend
2253  ip1 = ie(i)
2254  im1 = iw(i)
2255  ii = i + imb2
2256  if (ii > im) ii = ii - im
2257  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
2258  !& ! - (VWND(II,J,l)*COSL(II,J) &
2259  & - (vpoles(ii,1)*coslpoles(ii,1) &
2260  & + vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
2261  enddo
2262 !--
2263  ELSE !North pole point, compute at j=2
2264  jj = 2
2265  do i=ista,iend
2266  ip1 = ie(i)
2267  im1 = iw(i)
2268  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
2269  & + (vwnd(i,j,l)*cosl(i,j) &
2270  - vwnd(i,jj+1,l)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj)
2271  enddo
2272 !--
2273  ENDIF
2274  else
2275  IF(cosl(ista,j) >= small) THEN !not a pole point
2276  DO i=ista,iend
2277  ip1 = ie(i)
2278  im1 = iw(i)
2279  ii = i + imb2
2280  if (ii > im) ii = ii - im
2281  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
2282  !& ! + (VWND(II,J,l)*COSL(II,J) &
2283  & + (vpoles(ii,1)*coslpoles(ii,1) &
2284  & + vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
2285  enddo
2286 !--
2287  ELSE !North pole point, compute at j=2
2288  jj = 2
2289  do i=ista,iend
2290  ip1 = ie(i)
2291  im1 = iw(i)
2292  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
2293  & - (vwnd(i,j,l)*cosl(i,j) &
2294  - vwnd(i,jj+1,l)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj)
2295  enddo
2296  ENDIF
2297  endif
2298  ELSE IF(j == jm) THEN ! Near South pole
2299  if(gdlat(ista,j) < 0.) then ! count from north to south
2300  IF(cosl(ista,j) >= small) THEN !not a pole point
2301  DO i=ista,iend
2302  ip1 = ie(i)
2303  im1 = iw(i)
2304  ii = i + imb2
2305  if (ii > im) ii = ii - im
2306  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
2307  & + (vwnd(i,j-1,l)*cosl(i,j-1) &
2308  !& ! + VWND(II,J,l)*COSL(II,J))*wrk3(i,j)) * wrk1(i,j)
2309  & + vpoles(ii,2)*coslpoles(ii,2))*wrk3(i,j)) * wrk1(i,j)
2310  enddo
2311 !--
2312  ELSE !South pole point,compute at jm-1
2313  jj = jm-1
2314  do i=ista,iend
2315  ip1 = ie(i)
2316  im1 = iw(i)
2317  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
2318  & + (vwnd(i,jj-1,l)*cosl(i,jj-1) &
2319  & - vwnd(i,j,l)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj)
2320 
2321  enddo
2322  ENDIF
2323  else
2324  IF(cosl(ista,j) >= small) THEN !not a pole point
2325  DO i=ista,iend
2326  ip1 = ie(i)
2327  im1 = iw(i)
2328  ii = i + imb2
2329  if (ii > im) ii = ii - im
2330  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
2331  & - (vwnd(i,j-1,l)*cosl(i,j-1) &
2332  !& ! + VWND(II,J,l)*COSL(II,J))*wrk3(i,j)) * wrk1(i,j)
2333  & + vpoles(ii,2)*coslpoles(ii,2))*wrk3(i,j)) * wrk1(i,j)
2334  enddo
2335 !--
2336  ELSE !South pole point,compute at jm-1
2337  jj = jm-1
2338  do i=ista,iend
2339  ip1 = ie(i)
2340  im1 = iw(i)
2341  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
2342  & - (vwnd(i,jj-1,l)*cosl(i,jj-1) &
2343  & - vwnd(i,j,l)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj)
2344 
2345  enddo
2346  ENDIF
2347  endif
2348  ELSE
2349  DO i=ista,iend
2350  ip1 = ie(i)
2351  im1 = iw(i)
2352  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
2353  & + (vwnd(i,j-1,l)*cosl(i,j-1) &
2354  - vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
2355 !sk06132016
2356  if(div(i,j,l)>1.0)print*,'Debug in CALDIV',i,j,uwnd(ip1,j,l),uwnd(im1,j,l), &
2357  & wrk2(i,j),vwnd(i,j-1,l),cosl(i,j-1),vwnd(i,j+1,l),cosl(i,j+1), &
2358  & wrk3(i,j),wrk1(i,j),div(i,j,l)
2359 !--
2360  ENDDO
2361  ENDIF
2362  ENDDO ! end of J loop
2363 
2364 ! GFS use lon avg as one scaler value for pole point
2365 ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1,jsta),SPVAL,DIV(1,jsta,l))
2366 
2367  call exch(div(ista_2l:iend_2u,jsta_2l:jend_2u,l))
2368  call fullpole(div(ista_2l:iend_2u,jsta_2l:jend_2u,l),divpoles)
2369 
2370  cosltemp=spval
2371  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
2372  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
2373  divtemp=spval
2374  IF(jsta== 1) divtemp(1:im, 1)=divpoles(1:im,1)
2375  IF(jend==jm) divtemp(1:im,jm)=divpoles(1:im,2)
2376 
2377  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
2378  ,spval,divtemp(1:im,jsta:jend))
2379 
2380  IF(jsta== 1) div(ista:iend, 1,l)=divtemp(ista:iend, 1)
2381  IF(jend==jm) div(ista:iend,jm,l)=divtemp(ista:iend,jm)
2382 
2383 !sk06142016e
2384  if(div(ista,jsta,l)>1.0)print*,'Debug in CALDIV',jsta,div(ista,jsta,l)
2385 ! print*,'Debug in CALDIV',' jsta= ',jsta,DIV(1,jsta,l)
2386 
2387  enddo ! end of l looop
2388 !--
2389  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
2390 
2391 
2392  END SUBROUTINE caldiv
2393 
2394  SUBROUTINE calgradps(PS,PSX,PSY)
2395 
2410  use masks, only: gdlat, gdlon
2411  use params_mod, only: dtr, d00, small, erad
2412  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
2413  jsta, jend, im, jm, jsta_m, jend_m, &
2414  ista, iend, ista_m, iend_m, ista_2l, iend_2u
2415 
2416  use gridspec_mod, only: gridtype
2417 
2418  implicit none
2419 !
2420 ! DECLARE VARIABLES.
2421 !
2422  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in) :: ps
2423  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(inout) :: psx,psy
2424 !
2425  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
2426  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
2427 !
2428  integer i,j,ip1,im1,ii,iir,iil,jj,imb2
2429 !
2430 !***************************************************************************
2431 ! START CALGRADPS HERE.
2432 !
2433 ! LOOP TO COMPUTE ZONAL AND MERIDIONAL GRADIENTS OF PS OR LNPS
2434 !
2435 !sk06162016 DO J=JSTA_2L,JEND_2U
2436 !$omp parallel do private(i,j)
2437  DO j=jsta,jend
2438  DO i=ista,iend
2439  psx(i,j) = spval
2440  psy(i,j) = spval
2441 !sk PSX(I,J) = D00
2442 !sk PSY(I,J) = D00
2443  ENDDO
2444  ENDDO
2445 
2446  CALL exch(ps)
2447 
2448 ! IF (MODELNAME == 'GFS' .or. global) THEN
2449  CALL exch(gdlat(ista_2l,jsta_2l))
2450  CALL exch(gdlon(ista_2l,jsta_2l))
2451 
2452  allocate (wrk1(ista:iend,jsta:jend), wrk2(ista:iend,jsta:jend), &
2453  & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u))
2454  allocate(iw(im),ie(im))
2455 
2456  imb2 = im/2
2457 !$omp parallel do private(i)
2458  do i=ista,iend
2459  ie(i) = i+1
2460  iw(i) = i-1
2461  enddo
2462 ! iw(1) = im
2463 ! ie(im) = 1
2464 
2465 
2466 !$omp parallel do private(i,j,ip1,im1)
2467  DO j=jsta,jend
2468  do i=ista,iend
2469  ip1 = ie(i)
2470  im1 = iw(i)
2471  cosl(i,j) = cos(gdlat(i,j)*dtr)
2472  if(cosl(i,j) >= small) then
2473  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
2474  else
2475  wrk1(i,j) = 0.
2476  end if
2477  if(i == im .or. i == 1) then
2478  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
2479  else
2480  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
2481  end if
2482  enddo
2483  ENDDO
2484 
2485  CALL exch(cosl)
2486 
2487 !$omp parallel do private(i,j,ii)
2488  DO j=jsta,jend
2489  if (j == 1) then
2490  if(gdlat(ista,j) > 0.) then ! count from north to south
2491  do i=ista,iend
2492  ii = i + imb2
2493  if (ii > im) ii = ii - im
2494  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-gdlat(ii,j))*dtr) !1/dphi
2495  enddo
2496  else ! count from south to north
2497  do i=ista,iend
2498  ii = i + imb2
2499  if (ii > im) ii = ii - im
2500  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+gdlat(ii,j))*dtr) !1/dphi
2501  enddo
2502  end if
2503  elseif (j == jm) then
2504  if(gdlat(ista,j) < 0.) then ! count from north to south
2505  do i=ista,iend
2506  ii = i + imb2
2507  if (ii > im) ii = ii - im
2508  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+gdlat(ii,j))*dtr)
2509  enddo
2510  else ! count from south to north
2511  do i=ista,iend
2512  ii = i + imb2
2513  if (ii > im) ii = ii - im
2514  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-gdlat(ii,j))*dtr)
2515  enddo
2516  end if
2517  else
2518  do i=ista,iend
2519  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
2520  enddo
2521  endif
2522  ENDDO
2523 
2524 !$omp parallel do private(i,j,ip1,im1,ii,jj)
2525  DO j=jsta,jend
2526  IF(j == 1) then ! Near North pole
2527  if(gdlat(ista,j) > 0.) then ! count from north to south
2528  IF(cosl(ista,j) >= small) THEN !not a pole point
2529  DO i=ista,iend
2530  ip1 = ie(i)
2531  im1 = iw(i)
2532  ii = i + imb2
2533  if (ii > im) ii = ii - im
2534  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
2535  psy(i,j) = (ps(ii,j)-ps(i,j+1))*wrk3(i,j)/erad
2536  enddo
2537  ELSE !North pole point, compute at j=2
2538  jj = 2
2539  DO i=ista,iend
2540  ip1 = ie(i)
2541  im1 = iw(i)
2542  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
2543  psy(i,j) = (ps(i,j)-ps(i,jj+1))*wrk3(i,jj)/erad
2544  enddo
2545  ENDIF
2546  else
2547  IF(cosl(ista,j) >= small) THEN !not a pole point
2548  DO i=ista,iend
2549  ip1 = ie(i)
2550  im1 = iw(i)
2551  ii = i + imb2
2552  if (ii > im) ii = ii - im
2553  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
2554  psy(i,j) = - (ps(ii,j)-ps(i,j+1))*wrk3(i,j)/erad
2555  enddo
2556  ELSE !North pole point, compute at j=2
2557  jj = 2
2558  DO i=ista,iend
2559  ip1 = ie(i)
2560  im1 = iw(i)
2561  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
2562  psy(i,j) = - (ps(i,j)-ps(i,jj+1))*wrk3(i,jj)/erad
2563  enddo
2564  ENDIF
2565  endif
2566  ELSE IF(j == jm) THEN ! Near South pole
2567  if(gdlat(ista,j) < 0.) then ! count from north to south
2568  IF(cosl(ista,j) >= small) THEN !not a pole point
2569  DO i=ista,iend
2570  ip1 = ie(i)
2571  im1 = iw(i)
2572  ii = i + imb2
2573  if (ii > im) ii = ii - im
2574  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
2575  psy(i,j) = (ps(i,j-1)-ps(ii,j))*wrk3(i,j)/erad
2576  enddo
2577  ELSE !South pole point,compute at jm-1
2578  jj = jm-1
2579  DO i=ista,iend
2580  ip1 = ie(i)
2581  im1 = iw(i)
2582  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
2583  psy(i,j) = (ps(i,jj-1)-ps(i,j))*wrk3(i,jj)/erad
2584  enddo
2585  ENDIF
2586  else
2587  IF(cosl(ista,j) >= small) THEN !not a pole point
2588  DO i=ista,iend
2589  ip1 = ie(i)
2590  im1 = iw(i)
2591  ii = i + imb2
2592  if (ii > im) ii = ii - im
2593  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
2594  psy(i,j) = - (ps(i,j-1)-ps(ii,j))*wrk3(i,j)/erad
2595  enddo
2596  ELSE !South pole point,compute at jm-1
2597  jj = jm-1
2598  DO i=ista,iend
2599  ip1 = ie(i)
2600  im1 = iw(i)
2601  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
2602  psy(i,j) = - (ps(i,jj-1)-ps(i,j))*wrk3(i,jj)/erad
2603  enddo
2604  ENDIF
2605  endif
2606  ELSE
2607  DO i=ista,iend
2608  ip1 = ie(i)
2609  im1 = iw(i)
2610  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
2611  psy(i,j) = (ps(i,j-1)-ps(i,j+1))*wrk3(i,j)/erad
2612 !sk06142016A
2613  if(psx(i,j)>100.0)print*,'Debug in CALGRADPS: PSX',i,j,ps(ip1,j),ps(im1,j), &
2614 ! print*,'Debug in CALGRADPS',i,j,PS(ip1,J),PS(im1,J), &
2615  & wrk2(i,j),wrk1(i,j),psx(i,j)
2616  if(psy(i,j)>100.0)print*,'Debug in CALGRADPS: PSY',i,j,ps(i,j-1),ps(i,j+1), &
2617 ! print*,'Debug in CALGRADPS',i,j,PS(i,J-1),PS(i,J+1), &
2618  & wrk3(i,j),erad,psy(i,j)
2619 !--
2620  ENDDO
2621  END IF
2622 !
2623  ENDDO ! end of J loop
2624 
2625  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
2626 
2627 ! END IF
2628 
2629  END SUBROUTINE calgradps
2630 !
2631 !-------------------------------------------------------------------------------------
2632 !
2633  end module upp_physics
2634 
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
upp_math
dvdxdudy() computes dudy, dvdx, uwnd
Definition: UPP_MATH.f:17
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:345
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27
lookup_mod
Definition: LOOKUP.f:1