UPP v11.0.0
Loading...
Searching...
No Matches
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
upp_math
dvdxdudy() computes dudy, dvdx, uwnd
Definition UPP_MATH.f:17
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition UPP_PHYSICS.f:27
upp_physics::calcape2
subroutine, public calcape2(itype, dpbnd, p1d, t1d, q1d, l1d, cape, cins, lfc, esrhl, esrhh, dcape, dgld, esp)
calcape2() computes CAPE and CINS.
Definition UPP_PHYSICS.f:1006
upp_physics::calcape
subroutine, public calcape(itype, dpbnd, p1d, t1d, q1d, l1d, cape, cins, pparc, zeql, thund)
calcape() computes CAPE and CINS.
Definition UPP_PHYSICS.f:529
upp_physics::calgradps
subroutine, public calgradps(ps, psx, psy)
Definition UPP_PHYSICS.f:2395
upp_physics::calrh_nam
subroutine, public calrh_nam(p1, t1, q1, rh)
calrh_nam() computes relative humidity.
Definition UPP_PHYSICS.f:96
upp_physics::caldiv
subroutine, public caldiv(uwnd, vwnd, div)
CALDIV computes divergence.
Definition UPP_PHYSICS.f:2125
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition UPP_PHYSICS.f:346
upp_physics::calrh_gfs
subroutine, public calrh_gfs(p1, t1, q1, rh)
calrh_gfs() computes relative humidity.
Definition UPP_PHYSICS.f:176