UPP (upp-srw-2.2.0)
Loading...
Searching...
No Matches
MISCLN.f
Go to the documentation of this file.
1
2!
89 SUBROUTINE miscln
90
91!
92!
93 use vrbls3d, only: pmid, uh, vh, t, zmid, zint, pint, alpint, q, omga
94 use vrbls3d, only: catedr,mwt,gtg
95 use vrbls2d, only: pblh, cprate, fis, t500, t700, z500, z700,&
96 teql,ieql, cape,cin
97 use masks, only: lmh
98 use params_mod, only: d00, d50, h99999, h100, h1, h1m12, pq0, a2, a3, a4, &
99 rhmin, rgamog, tfrz, small, g
100 use ctlblk_mod, only: grib, cfld, fld_info, datapd, im, jsta, jend, jm, jsta_m, jend_m, &
101 nbnd, nbin_du, lm, htfd, spval, pthresh, nfd, petabnd, me,&
102 jsta_2l, jend_2u, modelname, submodelname, &
103 ista, iend, ista_m, iend_m, ista_2l, iend_2u, &
104 ifi_flight_levels
105 use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
106 use grib2_module, only: pset
107 use upp_physics, only: fpvsnew,calrh_pw,calcape,calcape2,tvirtual
108 use gridspec_mod, only: gridtype
109!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110 implicit none
111!
112! SET LOCAL PARAMETERS. MAKE SURE NFD AND NBND AGREE
113! WITH THE VALUES SET IN SUBROUTINES FDLVL AND BNDLYR,
114! RESPECTIVELY.
115 real,PARAMETER :: C2K=273.15
116 real,parameter :: con_rd =2.8705e+2 ! gas constant air (J/kg/K)
117 real,parameter :: con_rv =4.6150e+2 ! gas constant H2O
118 real,parameter :: con_eps =con_rd/con_rv
119 real,parameter :: con_epsm1 =con_rd/con_rv-1
120 real,parameter :: cpthresh =0.000004
121 real,PARAMETER :: D1000=1000
122 real,PARAMETER :: D1500=1500
123 real,PARAMETER :: D2000=2000
124 real,PARAMETER :: HCONST=42000000.
125 real,PARAMETER :: K2C=273.16
126 REAL,PARAMETER :: DM9999=-9999.0
127
128!
129! DECLARE VARIABLES.
130!
131 LOGICAL NORTH, FIELD1,FIELD2, NEED_IFI
132 LOGICAL, dimension(ISTA:IEND,JSTA:JEND) :: DONE, DONE1
133
134 INTEGER, allocatable :: LVLBND(:,:,:),LB2(:,:)
135! INTEGER LVLBND(IM,JM,NBND),LB2(IM,JM),LPBL(IM,JM)
136
137 real,dimension(im,jm) :: GRID1, GRID2
138 real,dimension(ista:iend,jsta:jend) :: P1D, T1D, Q1D, U1D, V1D, SHR1D, Z1D, &
139 rh1d, egrid1, egrid2, egrid3, egrid4, &
140 egrid5, egrid6, egrid7, egrid8, &
141 mlcape,mlcin,mllcl,mucape,mucin, &
142 freezelvl,muq1d,slcl,the,maxthe
143 integer,dimension(ista:iend,jsta:jend) :: MAXTHEPOS
144 real, dimension(:,:,:),allocatable :: OMGBND, PWTBND, QCNVBND, &
145 pbnd, tbnd, qbnd, &
146 ubnd, vbnd, rhbnd, &
147 wbnd, t7d, q7d, &
148 u7d, v6d, p7d, &
149 icingfd,gtgfd,catfd,mwtfd,midcal
150
151 real, dimension(:,:),allocatable :: QM8510, RH4710, RH8498, &
152 rh4796, rh1847, ust, vst, &
153 rh3310, rh6610, rh3366, &
154 pw3310, rh4410, rh7294, &
155 rh4472, &
156 t78483, t89671, p78483, p89671
157
158 REAL, dimension(:,:,:),allocatable :: HELI
159 real, dimension(:,:), allocatable :: USHR1, VSHR1, USHR6, VSHR6, &
160 maxwp, maxwz, maxwu, maxwv, &
161 maxwt
162 INTEGER,dimension(:,:),allocatable :: LLOW,LUPP,LLOW_ZINT,IEQL_ZINT, &
163 z_midcal
164 REAL, dimension(:,:),allocatable :: CANGLE,ESHR,UVECT,VVECT,&
165 effust,effvst,fshr,htsfc,&
166 esrh,z_temp
167!
168 integer I,J,ii,jj,L,ITYPE,ISVALUE,LBND,ILVL,IFD,ITYPEFDLVL(NFD), &
169 iget1, iget2, iget3, llmh,imax,jmax,lmax
170 real DPBND,PKL1,PKU1,FAC1,FAC2,PL,TL,QL,QSAT,RHL,TVRL,TVRBLO, &
171 es1,es2,qs1,qs2,rh1,rh2,zsf,depth(2),work1,work2,work3, &
172 scintmp,mucapetmp,mucintmp,mllcltmp,eshrtmp,mlcapetmp,stp,&
173 fshrtmp,mlcintmp,slcltmp,lapse,ship
174
175! Variables introduced to allow FD levels from control file - Y Mao
176 integer :: N,NFDCTL
177 REAL, allocatable :: HTFDCTL(:)
178 integer, allocatable :: ITYPEFDLVLCTL(:)
179 integer IE,IW,JN,JS,IVE(JM),IVW(JM),JVN,JVS
180 integer ISTART,ISTOP,JSTART,JSTOP
181 real dummy(ista:iend,jsta:jend)
182 integer idummy(ista:iend,jsta:jend)
183! NEW VARIABLES USED FOR EFFECTIVE LAYER
184 INTEGER,dimension(:,:),allocatable :: EL_BASE, EL_TOPS
185 LOGICAL,dimension(:,:),allocatable :: FOUND_BASE, FOUND_TOPS
186 INTEGER,dimension(:,:),allocatable :: L_THETAE_MAX
187 INTEGER,dimension(:,:),allocatable :: CAPE9, CINS9
188 CHARACTER(LEN=5) :: IM_CH, JSTA_CH, JEND_CH, ME_CH
189 CHARACTER(LEN=60) :: EFFL_FNAME
190 CHARACTER(LEN=60) :: EFFL_FNAME2
191 INTEGER :: IREC, IUNIT
192 INTEGER :: IREC2, IUNIT2
193 LOGICAL :: debugprint
194 INTEGER :: LLL
195 INTEGER :: LLCL_PAR, LEQL_PAR
196 REAL :: LMASK, PSFC, CAPE_PAR, CINS_PAR, LPAR0
197 REAL, DIMENSION(4) :: PARCEL0
198 REAL, DIMENSION(:), ALLOCATABLE :: TPAR_B, TPAR_T
199 REAL, DIMENSION(:), ALLOCATABLE :: TPAR_TMP
200 REAL, DIMENSION(:), ALLOCATABLE :: P_AMB, T_AMB, Q_AMB, ZINT_AMB
201 REAL, DIMENSION(:,:,:), ALLOCATABLE :: TPAR_BASE, TPAR_TOPS
202
203!
204!****************************************************************************
205! START MISCLN HERE.
206!
207 debugprint = .false.
208
209 need_ifi = iget(1007)>0 .or. iget(1008)>0 .or. iget(1009)>0 .or. iget(1010)>0
210
211 allocate(ushr1(ista_2l:iend_2u,jsta_2l:jend_2u),vshr1(ista_2l:iend_2u,jsta_2l:jend_2u), &
212 ushr6(ista_2l:iend_2u,jsta_2l:jend_2u),vshr6(ista_2l:iend_2u,jsta_2l:jend_2u))
213 allocate(ust(ista_2l:iend_2u,jsta_2l:jend_2u),vst(ista_2l:iend_2u,jsta_2l:jend_2u), &
214 heli(ista_2l:iend_2u,jsta_2l:jend_2u,2),fshr(ista_2l:iend_2u,jsta_2l:jend_2u))
215!
216! HELICITY AND STORM MOTION.
217 iget1 = iget(162)
218 iget2 = -1
219 iget3 = -1
220 if (iget1 > 0) then
221 iget2 = lvls(1,iget1)
222 iget3 = lvls(2,iget1)
223 endif
224 IF (iget1 > 0 .OR. iget(163) > 0 .OR. iget(164) > 0) THEN
225 depth(1) = 3000.0
226 depth(2) = 1000.0
227 CALL calhel(depth,ust,vst,heli,ushr1,vshr1,ushr6,vshr6)
228 IF (iget2 > 0) then
229!$omp parallel do private(i,j)
230 DO j=jsta,jend
231 DO i=ista,iend
232 grid1(i,j) = heli(i,j,1)
233 ENDDO
234 ENDDO
235 if(grib=='grib2') then
236 cfld=cfld+1
237 fld_info(cfld)%ifld=iavblfld(iget1)
238 fld_info(cfld)%lvl=lvlsxml(1,iget1)
239!$omp parallel do private(i,j,ii,jj)
240 do j=1,jend-jsta+1
241 jj = jsta+j-1
242 do i=1,iend-ista+1
243 ii = ista+i-1
244 datapd(i,j,cfld) = grid1(ii,jj)
245 enddo
246 enddo
247 endif
248 ENDIF
249
250 IF (iget3 > 0) then
251!$omp parallel do private(i,j)
252 DO j=jsta,jend
253 DO i=ista,iend
254 grid1(i,j) = heli(i,j,2)
255 ENDDO
256 ENDDO
257 if(grib=='grib2') then
258 cfld=cfld+1
259 fld_info(cfld)%ifld=iavblfld(iget1)
260 fld_info(cfld)%lvl=lvlsxml(2,iget1)
261!$omp parallel do private(i,j,ii,jj)
262 do j=1,jend-jsta+1
263 jj = jsta+j-1
264 do i=1,iend-ista+1
265 ii = ista+i-1
266 datapd(i,j,cfld) = grid1(ii,jj)
267 enddo
268 enddo
269 endif
270 ENDIF
271
272 IF (iget(163) > 0) THEN
273!$omp parallel do private(i,j)
274 DO j=jsta,jend
275 DO i=ista,iend
276 grid1(i,j) = ust(i,j)
277 ENDDO
278 ENDDO
279 if(grib=='grib2') then
280 cfld=cfld+1
281 fld_info(cfld)%ifld=iavblfld(iget(163))
282!$omp parallel do private(i,j,ii,jj)
283 do j=1,jend-jsta+1
284 jj = jsta+j-1
285 do i=1,iend-ista+1
286 ii = ista+i-1
287 datapd(i,j,cfld) = grid1(ii,jj)
288 enddo
289 enddo
290 endif
291 ENDIF
292 IF (iget(164) > 0) THEN
293!$omp parallel do private(i,j)
294 DO j=jsta,jend
295 DO i=ista,iend
296 grid1(i,j) = vst(i,j)
297 ENDDO
298 ENDDO
299 if(grib=='grib2') then
300 cfld=cfld+1
301 fld_info(cfld)%ifld=iavblfld(iget(164))
302!$omp parallel do private(i,j,ii,jj)
303 do j=1,jend-jsta+1
304 jj = jsta+j-1
305 do i=1,iend-ista+1
306 ii = ista+i-1
307 datapd(i,j,cfld) = grid1(ii,jj)
308 enddo
309 enddo
310 endif
311 ENDIF
312 ENDIF
313
314! UPDRAFT HELICITY
315
316 if (iget(427) > 0) THEN
317 CALL calupdhel(grid1(ista_2l:iend_2u,jsta_2l:jend_2u))
318 if(grib=='grib2') then
319 cfld=cfld+1
320 fld_info(cfld)%ifld=iavblfld(iget(427))
321!$omp parallel do private(i,j,ii,jj)
322 do j=1,jend-jsta+1
323 jj = jsta+j-1
324 do i=1,iend-ista+1
325 ii = ista+i-1
326 datapd(i,j,cfld) = grid1(ii,jj)
327 enddo
328 enddo
329 endif
330
331 ENDIF
332
333! CRA 0-1 KM AND 0-6 KM SHEAR
334
335 IF(iget(430) > 0 .OR. iget(431) > 0 .OR. iget(432) > 0 &
336 .OR. iget(433) > 0) THEN
337
338 depth = 6000.0
339 CALL calhel(depth,ust,vst,heli,ushr1,vshr1,ushr6,vshr6)
340! 0-6 km shear magnitude
341!$omp parallel do private(i,j)
342 DO j=jsta,jend
343 DO i=ista,iend
344 fshr(i,j) = sqrt(ushr6(i,j)**2+vshr6(i,j)**2)
345 ENDDO
346 ENDDO
347 IF(iget(430) > 0) THEN
348!$omp parallel do private(i,j)
349 DO j=jsta,jend
350 DO i=ista,iend
351 grid1(i,j) = ushr1(i,j)
352 ENDDO
353 ENDDO
354 if(grib=='grib2') then
355 cfld=cfld+1
356 fld_info(cfld)%ifld=iavblfld(iget(430))
357!$omp parallel do private(i,j,ii,jj)
358 do j=1,jend-jsta+1
359 jj = jsta+j-1
360 do i=1,iend-ista+1
361 ii = ista+i-1
362 datapd(i,j,cfld) = grid1(ii,jj)
363 enddo
364 enddo
365 endif
366 ENDIF
367 IF(iget(431) > 0) THEN
368!$omp parallel do private(i,j)
369 DO j=jsta,jend
370 DO i=ista,iend
371 grid1(i,j) = vshr1(i,j)
372 ENDDO
373 ENDDO
374 if(grib=='grib2') then
375 cfld=cfld+1
376 fld_info(cfld)%ifld=iavblfld(iget(431))
377!$omp parallel do private(i,j,ii,jj)
378 do j=1,jend-jsta+1
379 jj = jsta+j-1
380 do i=1,iend-ista+1
381 ii = ista+i-1
382 datapd(i,j,cfld) = grid1(ii,jj)
383 enddo
384 enddo
385 endif
386 ENDIF
387 IF(iget(432) > 0) THEN
388!$omp parallel do private(i,j)
389 DO j=jsta,jend
390 DO i=ista,iend
391 grid1(i,j) = ushr6(i,j)
392 ENDDO
393 ENDDO
394 if(grib=='grib2') then
395 cfld=cfld+1
396 fld_info(cfld)%ifld=iavblfld(iget(432))
397!$omp parallel do private(i,j,ii,jj)
398 do j=1,jend-jsta+1
399 jj = jsta+j-1
400 do i=1,iend-ista+1
401 ii = ista+i-1
402 datapd(i,j,cfld) = grid1(ii,jj)
403 enddo
404 enddo
405 endif
406 ENDIF
407 IF(iget(433) > 0) THEN
408!$omp parallel do private(i,j)
409 DO j=jsta,jend
410 DO i=ista,iend
411 grid1(i,j) = vshr6(i,j)
412 ENDDO
413 ENDDO
414 if(grib=='grib2') then
415 cfld=cfld+1
416 fld_info(cfld)%ifld=iavblfld(iget(433))
417!$omp parallel do private(i,j,ii,jj)
418 do j=1,jend-jsta+1
419 jj = jsta+j-1
420 do i=1,iend-ista+1
421 ii = ista+i-1
422 datapd(i,j,cfld) = grid1(ii,jj)
423 enddo
424 enddo
425 endif
426 ENDIF
427 ENDIF
428!
429 if (allocated(ushr1)) deallocate(ushr1)
430 if (allocated(vshr1)) deallocate(vshr1)
431 if (allocated(ushr6)) deallocate(ushr6)
432 if (allocated(vshr6)) deallocate(vshr6)
433 if (allocated(ust)) deallocate(ust)
434 if (allocated(vst)) deallocate(vst)
435 if (allocated(heli)) deallocate(heli)
436! CRA
437!
438!
439! ***BLOCK 1: TROPOPAUSE P, Z, T, U, V, AND WIND SHEAR.
440!
441 IF ((iget(054)>0).OR.(iget(055)>0).OR. &
442 (iget(056)>0).OR.(iget(057)>0).OR. &
443 (iget(177)>0).OR. &
444 (iget(058)>0).OR.(iget(108)>0) ) THEN
445! Chuang: Use GFS algorithm per Iredell's and DiMego's decision on unification
446!$omp parallel do private(i,j)
447 DO j=jsta,jend
448 DO i=ista,iend
449
450 if(pmid(i,j,1)<spval) then
451! INPUT
452 CALL tpause(lm,pmid(i,j,1:lm),uh(i,j,1:lm) &
453! INPUT
454 ,vh(i,j,1:lm),t(i,j,1:lm),zmid(i,j,1:lm) &
455! OUTPUT
456 ,p1d(i,j),u1d(i,j),v1d(i,j),t1d(i,j) &
457! OUTPUT
458 ,z1d(i,j),shr1d(i,j)) ! OUTPUT
459 else
460 p1d(i,j) = spval
461 u1d(i,j) = spval
462 v1d(i,j) = spval
463 t1d(i,j) = spval
464 z1d(i,j) = spval
465 shr1d(i,j) = spval
466 endif
467
468 END DO
469 END DO
470!
471! TROPOPAUSE PRESSURE.
472 IF (iget(054) > 0) THEN
473!$omp parallel do private(i,j)
474 DO j=jsta,jend
475 DO i=ista,iend
476 grid1(i,j) = p1d(i,j)
477 ENDDO
478 ENDDO
479 if(grib=='grib2') then
480 cfld=cfld+1
481 fld_info(cfld)%ifld=iavblfld(iget(054))
482!$omp parallel do private(i,j,ii,jj)
483 do j=1,jend-jsta+1
484 jj = jsta+j-1
485 do i=1,iend-ista+1
486 ii = ista+i-1
487 datapd(i,j,cfld) = grid1(ii,jj)
488 enddo
489 enddo
490 endif
491 ENDIF
492
493! ICAO HEIGHT OF TROPOPAUSE
494 IF (iget(399)>0) THEN
495 CALL icaoheight(p1d, grid1(ista:iend,jsta:jend))
496! print*,'sample TROPOPAUSE ICAO HEIGHTS',GRID1(im/2,(jsta+jend)/2)
497 if(grib=='grib2') then
498 cfld=cfld+1
499 fld_info(cfld)%ifld=iavblfld(iget(399))
500!$omp parallel do private(i,j,ii,jj)
501 do j=1,jend-jsta+1
502 jj = jsta+j-1
503 do i=1,iend-ista+1
504 ii = ista+i-1
505 datapd(i,j,cfld) = grid1(ii,jj)
506 enddo
507 enddo
508 endif
509 ENDIF
510
511! TROPOPAUSE HEIGHT.
512 IF (iget(177) > 0) THEN
513!$omp parallel do private(i,j)
514 DO j=jsta,jend
515 DO i=ista,iend
516 grid1(i,j) = z1d(i,j)
517 ENDDO
518 ENDDO
519 if(grib=='grib2') then
520 cfld=cfld+1
521 fld_info(cfld)%ifld=iavblfld(iget(177))
522!$omp parallel do private(i,j,ii,jj)
523 do j=1,jend-jsta+1
524 jj = jsta+j-1
525 do i=1,iend-ista+1
526 ii = ista+i-1
527 datapd(i,j,cfld) = grid1(ii,jj)
528 enddo
529 enddo
530 endif
531 ENDIF
532!
533! TROPOPAUSE TEMPERATURE.
534 IF (iget(055) > 0) THEN
535!$omp parallel do private(i,j)
536 DO j=jsta,jend
537 DO i=ista,iend
538 grid1(i,j) = t1d(i,j)
539 ENDDO
540 ENDDO
541 if(grib=='grib2') then
542 cfld=cfld+1
543 fld_info(cfld)%ifld=iavblfld(iget(055))
544!$omp parallel do private(i,j,ii,jj)
545 do j=1,jend-jsta+1
546 jj = jsta+j-1
547 do i=1,iend-ista+1
548 ii = ista+i-1
549 datapd(i,j,cfld) = grid1(ii,jj)
550 enddo
551 enddo
552 endif
553 ENDIF
554!
555! TROPOPAUSE POTENTIAL TEMPERATURE.
556 IF (iget(108) > 0) THEN
557 CALL calpot(p1d,t1d,grid1(ista:iend,jsta:jend))
558 if(grib=='grib2') then
559 cfld=cfld+1
560 fld_info(cfld)%ifld=iavblfld(iget(108))
561!$omp parallel do private(i,j,ii,jj)
562 do j=1,jend-jsta+1
563 jj = jsta+j-1
564 do i=1,iend-ista+1
565 ii = ista+i-1
566 datapd(i,j,cfld) = grid1(ii,jj)
567 enddo
568 enddo
569 endif
570 ENDIF
571!
572! TROPOPAUSE U WIND AND/OR V WIND.
573 IF ((iget(056) > 0).OR.(iget(057) > 0)) THEN
574!$omp parallel do private(i,j)
575 DO j=jsta,jend
576 DO i=ista,iend
577 grid1(i,j)=u1d(i,j)
578 grid2(i,j)=v1d(i,j)
579 ENDDO
580 ENDDO
581 if(grib=='grib2') then
582 if(iget(056)>0) then
583 cfld=cfld+1
584 fld_info(cfld)%ifld=iavblfld(iget(056))
585!$omp parallel do private(i,j,ii,jj)
586 do j=1,jend-jsta+1
587 jj = jsta+j-1
588 do i=1,iend-ista+1
589 ii = ista+i-1
590 datapd(i,j,cfld) = grid1(ii,jj)
591 enddo
592 enddo
593 endif
594 if(iget(057)>0) then
595 cfld=cfld+1
596 fld_info(cfld)%ifld=iavblfld(iget(057))
597!$omp parallel do private(i,j,ii,jj)
598 do j=1,jend-jsta+1
599 jj = jsta+j-1
600 do i=1,iend-ista+1
601 ii = ista+i-1
602 datapd(i,j,cfld) = grid2(ii,jj)
603 enddo
604 enddo
605 endif
606 endif
607 ENDIF
608!
609! TROPOPAUSE WIND SHEAR.
610 IF (iget(058) > 0) THEN
611!$omp parallel do private(i,j)
612 DO j=jsta,jend
613 DO i=ista,iend
614 grid1(i,j) = shr1d(i,j)
615 ENDDO
616 ENDDO
617 if(grib=='grib2') then
618 cfld=cfld+1
619 fld_info(cfld)%ifld=iavblfld(iget(058))
620!$omp parallel do private(i,j,ii,jj)
621 do j=1,jend-jsta+1
622 jj = jsta+j-1
623 do i=1,iend-ista+1
624 ii = ista+i-1
625 datapd(i,j,cfld) = grid1(ii,jj)
626 enddo
627 enddo
628 endif
629 ENDIF
630 ENDIF
631!
632!
633! ***BLOCK 2: MAX WIND LEVEL P, Z, U, AND V
634!
635! MAX WIND LEVEL CALCULATIONS
636 IF ((iget(173)>0) .OR. (iget(174)>0) .OR. &
637 (iget(175)>0) .OR. (iget(176)>0)) THEN
638
639 allocate(maxwp(ista:iend,jsta:jend), maxwz(ista:iend,jsta:jend), &
640 maxwu(ista:iend,jsta:jend), maxwv(ista:iend,jsta:jend),maxwt(ista:iend,jsta:jend))
641!$omp parallel do private(i,j)
642 DO j=jsta,jend
643 DO i=ista,iend
644 maxwp(i,j)=spval
645 maxwz(i,j)=spval
646 maxwu(i,j)=spval
647 maxwv(i,j)=spval
648 ENDDO
649 ENDDO
650
651! CALL CALMXW(MAXWP,MAXWZ,MAXWU,MAXWV,MAXWT)
652! Chuang: Use GFS algorithm per Iredell's and DiMego's decision on unification
653!$omp parallel do private(i,j)
654 DO j=jsta,jend
655 loopi:DO i=ista,iend
656 DO l=1,lm
657 IF (abs(pmid(i,j,l)-spval)<=small .OR. &
658 abs(uh(i,j,l)-spval)<=small .OR. &
659 abs(uh(i,j,l)-spval)<=small .OR. &
660 abs(vh(i,j,l)-spval)<=small .OR. &
661 abs(t(i,j,l)-spval)<=small .OR. &
662 abs(zmid(i,j,l)-spval)<=small) cycle loopi
663 ENDDO
664! INPUT
665 CALL mxwind(lm,pmid(i,j,1:lm),uh(i,j,1:lm) &
666! INPUT
667 ,vh(i,j,1:lm),t(i,j,1:lm),zmid(i,j,1:lm) &
668! OUTPUT
669 ,maxwp(i,j),maxwu(i,j),maxwv(i,j) &
670! OUTPUT
671 ,maxwt(i,j),maxwz(i,j))
672 ENDDO loopi
673 END DO
674! PRESSURE OF MAX WIND LEVEL
675 IF (iget(173) > 0) THEN
676!$omp parallel do private(i,j)
677 DO j=jsta,jend
678 DO i=ista,iend
679 grid1(i,j) = maxwp(i,j)
680 ENDDO
681 ENDDO
682 if(grib=='grib2') then
683 cfld=cfld+1
684 fld_info(cfld)%ifld=iavblfld(iget(173))
685!$omp parallel do private(i,j,ii,jj)
686 do j=1,jend-jsta+1
687 jj = jsta+j-1
688 do i=1,iend-ista+1
689 ii = ista+i-1
690 datapd(i,j,cfld) = grid1(ii,jj)
691 enddo
692 enddo
693 endif
694 ENDIF
695! ICAO HEIGHT OF MAX WIND LEVEL
696 IF (iget(398)>0) THEN
697 CALL icaoheight(maxwp, grid1(ista:iend,jsta:jend))
698! print*,'sample MAX WIND ICAO HEIGHTS',GRID1(im/2,(jsta+jend)/2)
699 if(grib=='grib2') then
700 cfld=cfld+1
701 fld_info(cfld)%ifld=iavblfld(iget(398))
702!$omp parallel do private(i,j,ii,jj)
703 do j=1,jend-jsta+1
704 jj = jsta+j-1
705 do i=1,iend-ista+1
706 ii = ista+i-1
707 datapd(i,j,cfld) = grid1(ii,jj)
708 enddo
709 enddo
710 endif
711 ENDIF
712! HEIGHT OF MAX WIND LEVEL
713 IF (iget(174) > 0) THEN
714!$omp parallel do private(i,j)
715 DO j=jsta,jend
716 DO i=ista,iend
717 grid1(i,j) = maxwz(i,j)
718 ENDDO
719 ENDDO
720 if(grib=='grib2') then
721 cfld=cfld+1
722 fld_info(cfld)%ifld=iavblfld(iget(174))
723!$omp parallel do private(i,j,ii,jj)
724 do j=1,jend-jsta+1
725 jj = jsta+j-1
726 do i=1,iend-ista+1
727 ii = ista+i-1
728 datapd(i,j,cfld) = grid1(ii,jj)
729 enddo
730 enddo
731 endif
732 ENDIF
733
734! MAX WIND LEVEL U WIND AND/OR V WIND.
735 IF ((iget(175) > 0).OR.(iget(176) > 0)) THEN
736!$omp parallel do private(i,j)
737 DO j=jsta,jend
738 DO i=ista,iend
739 grid1(i,j) = maxwu(i,j)
740 grid2(i,j) = maxwv(i,j)
741 ENDDO
742 ENDDO
743 if(grib=='grib2') then
744 cfld=cfld+1
745 fld_info(cfld)%ifld=iavblfld(iget(175))
746!$omp parallel do private(i,j,ii,jj)
747 do j=1,jend-jsta+1
748 jj = jsta+j-1
749 do i=1,iend-ista+1
750 ii = ista+i-1
751 datapd(i,j,cfld) = grid1(ii,jj)
752 enddo
753 enddo
754 cfld=cfld+1
755 fld_info(cfld)%ifld=iavblfld(iget(176))
756!$omp parallel do private(i,j,ii,jj)
757 do j=1,jend-jsta+1
758 jj = jsta+j-1
759 do i=1,iend-ista+1
760 ii = ista+i-1
761 datapd(i,j,cfld) = grid2(ii,jj)
762 enddo
763 enddo
764 endif
765 ENDIF
766! TEMPERATURE OF MAX WIND LEVEL
767 IF (iget(314) > 0) THEN
768!$omp parallel do private(i,j)
769 DO j=jsta,jend
770 DO i=ista,iend
771 grid1(i,j)=maxwt(i,j)
772 ENDDO
773 ENDDO
774 if(grib=='grib2') then
775 cfld=cfld+1
776 fld_info(cfld)%ifld=iavblfld(iget(314))
777!$omp parallel do private(i,j,ii,jj)
778 do j=1,jend-jsta+1
779 jj = jsta+j-1
780 do i=1,iend-ista+1
781 ii = ista+i-1
782 datapd(i,j,cfld) = grid1(ii,jj)
783 enddo
784 enddo
785 endif
786 ENDIF
787 deallocate(maxwp,maxwz,maxwu,maxwv,maxwt)
788 ENDIF
789
790!
791! ***BLOCK 3-1: FD LEVEL (selected) T, Q, U, AND V.
792!
793 IF ( (iget(059)>0.or.iget(586)>0).OR.iget(911)>0.OR. &
794 (iget(060)>0.or.iget(576)>0).OR. &
795 (iget(061)>0.or.iget(577)>0).OR. &
796 (iget(451)>0.or.iget(578)>0).OR.iget(580)>0 ) THEN
797
798 ALLOCATE(t7d(ista:iend,jsta:jend,nfd), q7d(ista:iend,jsta:jend,nfd), &
799 u7d(ista:iend,jsta:jend,nfd), v6d(ista:iend,jsta:jend,nfd), &
800 p7d(ista:iend,jsta:jend,nfd), icingfd(ista:iend,jsta:jend,nfd))
801
802!
803! DETERMINE WHETHER TO DO MSL OR AGL FD LEVELS
804!
805 itypefdlvl=1
806 DO ifd = 1,nfd
807 IF (iget(059)>0) THEN
808 IF (lvls(ifd,iget(059))>1) itypefdlvl(ifd)=2
809 ENDIF
810 IF (iget(911)>0) THEN
811 IF (lvls(ifd,iget(911))>1) itypefdlvl(ifd)=2
812 ENDIF
813!for grib2, spec hgt only
814 IF (iget(586)>0) THEN
815 IF(lvls(ifd,iget(586))>0) itypefdlvl(ifd)=2
816 ENDIF
817 IF (iget(060)>0) THEN
818 IF (lvls(ifd,iget(060))>1) itypefdlvl(ifd)=2
819 ENDIF
820 IF (iget(576)>0) THEN
821 IF(lvls(ifd,iget(576))>0) itypefdlvl(ifd)=2
822 ENDIF
823 IF (iget(061)>0) THEN
824 IF (lvls(ifd,iget(061))>1) itypefdlvl(ifd)=2
825 ENDIF
826 IF (iget(577)>0) then
827 if(lvls(ifd,iget(577))>0) itypefdlvl(ifd)=2
828 ENDIF
829 IF (iget(451)>0) THEN
830 IF (lvls(ifd,iget(451))>1) itypefdlvl(ifd)=2
831 ENDIF
832 IF (iget(578)>0) then
833 if(lvls(ifd,iget(578))>0) itypefdlvl(ifd)=2
834 ENDIF
835
836 IF (iget(580)>0) then
837 if(lvls(ifd,iget(580))>1) itypefdlvl(ifd)=2
838 ENDIF
839 IF (iget(587)>0) then
840 if(lvls(ifd,iget(587))>0) itypefdlvl(ifd)=2
841 ENDIF
842
843 ENDDO
844! print *,'call FDLVL with ITYPEFDLVL: ', ITYPEFDLVL,'for tmp,lvls=',LVLS(1:15,iget(59)), &
845! 'grib2tmp lvs=',LVLS(1:15,iget(586))
846
847 CALL fdlvl(itypefdlvl,t7d,q7d,u7d,v6d,p7d,icingfd)
848!
849 loop_10: DO ifd = 1,nfd
850!
851! FD LEVEL TEMPERATURE.
852 iget1 = iget(059)
853 iget2 = iget(586)
854 if (iget1 > 0) then
855 work1 = lvls(ifd,iget1)
856 else
857 work1 = 0.0
858 endif
859 if (iget2 > 0) then
860 work2 = lvls(ifd,iget2)
861 else
862 work2 = 0.0
863 endif
864 IF (iget1 > 0 .or. iget2 > 0) THEN
865 IF (work1 > 0 .or. work2 > 0) THEN
866
867!$omp parallel do private(i,j)
868 DO j=jsta,jend
869 DO i=ista,iend
870 grid1(i,j) = t7d(i,j,ifd)
871 ENDDO
872 ENDDO
873 IF(work1 > 0) then
874 if(grib == 'grib2') then
875 cfld = cfld + 1
876 fld_info(cfld)%ifld = iavblfld(iget1)
877 fld_info(cfld)%lvl = lvlsxml(ifd,iget1)
878!$omp parallel do private(i,j,ii,jj)
879 do j=1,jend-jsta+1
880 jj = jsta+j-1
881 do i=1,iend-ista+1
882 ii = ista+i-1
883 datapd(i,j,cfld) = grid1(ii,jj)
884 enddo
885 enddo
886 endif
887 ENDIF
888 IF (work2 > 0) THEN
889 if(grib == 'grib2') then
890 cfld = cfld + 1
891 fld_info(cfld)%ifld = iavblfld(iget2)
892 fld_info(cfld)%lvl = lvlsxml(ifd,iget2)
893!$omp parallel do private(i,j,ii,jj)
894 do j=1,jend-jsta+1
895 jj = jsta+j-1
896 do i=1,iend-ista+1
897 ii = ista+i-1
898 datapd(i,j,cfld) = grid1(ii,jj)
899 enddo
900 enddo
901 endif
902 ENDIF
903 ENDIF
904 ENDIF
905
906! FD LEVEL VIRTUAL TEMPERATURE.
907 IF (iget(911)>0) THEN
908 IF (lvls(ifd,iget(911))>0) THEN
909 DO j=jsta,jend
910 DO i=ista,iend
911 if ( t7d(i,j,ifd) > 600 ) then
912 grid1(i,j)=spval
913 else
914 grid1(i,j)=t7d(i,j,ifd)*(1.+0.608*q7d(i,j,ifd))
915 endif
916 !print *, "grid value ",T7D(I,J,IFD),Q7D(I,J,IFD),T7D(I,J,IFD)*(1.+0.608*Q7D(I,J,IFD)),GRID1(I,J)
917 ENDDO
918 ENDDO
919 IF(lvls(ifd,iget(911))>0) then
920 if(grib=='grib2') then
921 cfld=cfld+1
922 fld_info(cfld)%ifld=iavblfld(iget(911))
923 fld_info(cfld)%lvl=lvlsxml(ifd,iget(911))
924 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
925 endif
926 ENDIF
927 ENDIF
928 ENDIF
929
930!
931! FD LEVEL SPEC HUMIDITY.
932 iget1 = iget(451)
933 iget2 = iget(578)
934 if (iget1 > 0) then
935 work1 = lvls(ifd,iget1)
936 else
937 work1 = 0.0
938 endif
939 if (iget2 > 0) then
940 work2 = lvls(ifd,iget2)
941 else
942 work2 = 0.0
943 endif
944 IF (iget1 > 0 .or. iget2 > 0) THEN
945 IF (work1 > 0 .or. work2 > 0)THEN
946!$omp parallel do private(i,j)
947 DO j=jsta,jend
948 DO i=ista,iend
949 grid1(i,j) = q7d(i,j,ifd)
950 ENDDO
951 ENDDO
952 if(work1 > 0) then
953 if(grib == 'grib2') then
954 cfld = cfld + 1
955 fld_info(cfld)%ifld = iavblfld(iget1)
956 fld_info(cfld)%lvl = lvlsxml(ifd,iget1)
957!$omp parallel do private(i,j,ii,jj)
958 do j=1,jend-jsta+1
959 jj = jsta+j-1
960 do i=1,iend-ista+1
961 ii = ista+i-1
962 datapd(i,j,cfld) = grid1(ii,jj)
963 enddo
964 enddo
965 endif
966 endif
967 if(work2 > 0) then
968 if(grib == 'grib2') then
969 cfld = cfld + 1
970 fld_info(cfld)%ifld = iavblfld(iget2)
971 fld_info(cfld)%lvl = lvlsxml(ifd,iget2)
972!$omp parallel do private(i,j,ii,jj)
973 do j=1,jend-jsta+1
974 jj = jsta+j-1
975 do i=1,iend-ista+1
976 ii = ista+i-1
977 datapd(i,j,cfld) = grid1(ii,jj)
978 enddo
979 enddo
980 endif
981 endif
982 ENDIF
983 ENDIF
984!
985! FD LEVEL PRESSURE
986 iget1 = iget(482)
987 iget2 = iget(579)
988 if (iget1 > 0) then
989 work1 = lvls(ifd,iget1)
990 else
991 work1 = 0.0
992 endif
993 if (iget2 > 0) then
994 work2 = lvls(ifd,iget2)
995 else
996 work2 = 0.0
997 endif
998 IF (iget1 > 0 .or. iget2 > 0) THEN
999 IF (work1 > 0 .or. work2 > 0) THEN
1000!$omp parallel do private(i,j)
1001 DO j=jsta,jend
1002 DO i=ista,iend
1003 grid1(i,j) = p7d(i,j,ifd)
1004 ENDDO
1005 ENDDO
1006 if(work1 > 0) then
1007 if(grib == 'grib2') then
1008 cfld = cfld + 1
1009 fld_info(cfld)%ifld = iavblfld(iget1)
1010 fld_info(cfld)%lvl = lvlsxml(ifd,iget1)
1011!$omp parallel do private(i,j,ii,jj)
1012 do j=1,jend-jsta+1
1013 jj = jsta+j-1
1014 do i=1,iend-ista+1
1015 ii = ista+i-1
1016 datapd(i,j,cfld) = grid1(ii,jj)
1017 enddo
1018 enddo
1019 endif
1020 endif
1021 if(work2 > 0) then
1022 if(grib == 'grib2') then
1023 cfld = cfld + 1
1024 fld_info(cfld)%ifld = iavblfld(iget2)
1025 fld_info(cfld)%lvl = lvlsxml(ifd,iget2)
1026!$omp parallel do private(i,j,ii,jj)
1027 do j=1,jend-jsta+1
1028 jj = jsta+j-1
1029 do i=1,iend-ista+1
1030 ii = ista+i-1
1031 datapd(i,j,cfld) = grid1(ii,jj)
1032 enddo
1033 enddo
1034 endif
1035 endif
1036 ENDIF
1037 ENDIF
1038!
1039! FD LEVEL ICING
1040 iget1 = iget(580)
1041 iget2 = iget(587)
1042 if (iget1 > 0) then
1043 work1 = lvls(ifd,iget1)
1044 else
1045 work1 = 0.0
1046 endif
1047 if (iget2 > 0) then
1048 work2 = lvls(ifd,iget2)
1049 else
1050 work2 = 0.0
1051 endif
1052 IF (iget1 > 0 .or. iget2 > 0) THEN
1053 IF (work1 > 0 .or. work2 > 0) THEN
1054!$omp parallel do private(i,j)
1055 DO j=jsta,jend
1056 DO i=ista,iend
1057 grid1(i,j) = icingfd(i,j,ifd)
1058 ENDDO
1059 ENDDO
1060 if(iget1 > 0) then
1061 if(grib == 'grib2') then
1062 cfld = cfld + 1
1063 fld_info(cfld)%ifld = iavblfld(iget1)
1064 fld_info(cfld)%lvl = lvlsxml(ifd,iget1)
1065!$omp parallel do private(i,j,ii,jj)
1066 do j=1,jend-jsta+1
1067 jj = jsta+j-1
1068 do i=1,iend-ista+1
1069 ii = ista+i-1
1070 datapd(i,j,cfld) = grid1(ii,jj)
1071 enddo
1072 enddo
1073 endif
1074 endif
1075 if(work2 > 0) then
1076 if(grib == 'grib2') then
1077 cfld = cfld + 1
1078 fld_info(cfld)%ifld = iavblfld(iget2)
1079 fld_info(cfld)%lvl = lvlsxml(ifd,iget2)
1080!$omp parallel do private(i,j,ii,jj)
1081 do j=1,jend-jsta+1
1082 jj = jsta+j-1
1083 do i=1,iend-ista+1
1084 ii = ista+i-1
1085 datapd(i,j,cfld) = grid1(ii,jj)
1086 enddo
1087 enddo
1088 endif
1089 endif
1090
1091 ENDIF
1092 ENDIF
1093!
1094
1095!
1096!
1097! FD LEVEL U WIND AND/OR V WIND.
1098 IF ((iget(060)>0).OR.(iget(061)>0)) THEN
1099!$omp parallel do private(i,j)
1100 DO j=jsta,jend
1101 DO i=ista,iend
1102 grid1(i,j)=u7d(i,j,ifd)
1103 grid2(i,j)=v6d(i,j,ifd)
1104 ENDDO
1105 ENDDO
1106 IF (iget(060)>0) THEN
1107 IF (lvls(ifd,iget(060))>0) then
1108 if(grib=='grib2') then
1109 cfld=cfld+1
1110 fld_info(cfld)%ifld=iavblfld(iget(060))
1111 fld_info(cfld)%lvl=lvlsxml(ifd,iget(060))
1112!$omp parallel do private(i,j,ii,jj)
1113 do j=1,jend-jsta+1
1114 jj = jsta+j-1
1115 do i=1,iend-ista+1
1116 ii = ista+i-1
1117 datapd(i,j,cfld) = grid1(ii,jj)
1118 enddo
1119 enddo
1120 endif
1121 ENDIF
1122 ENDIF
1123 IF (iget(061)>0) THEN
1124 IF (lvls(ifd,iget(061))>0) THEN
1125 if(grib=='grib2') then
1126 cfld=cfld+1
1127 fld_info(cfld)%ifld=iavblfld(iget(061))
1128 fld_info(cfld)%lvl=lvlsxml(ifd,iget(061))
1129!$omp parallel do private(i,j,ii,jj)
1130 do j=1,jend-jsta+1
1131 jj = jsta+j-1
1132 do i=1,iend-ista+1
1133 ii = ista+i-1
1134 datapd(i,j,cfld) = grid2(ii,jj)
1135 enddo
1136 enddo
1137 endif
1138 ENDIF
1139 ENDIF
1140 ENDIF
1141!
1142! FD LEVEL U WIND AND/OR V WIND.
1143 IF ((iget(576)>0).OR.(iget(577)>0)) THEN
1144!$omp parallel do private(i,j)
1145 DO j=jsta,jend
1146 DO i=ista,iend
1147 grid1(i,j) = u7d(i,j,ifd)
1148 grid2(i,j) = v6d(i,j,ifd)
1149 ENDDO
1150 ENDDO
1151 IF (iget(576)>0) THEN
1152 IF (lvls(ifd,iget(576))>0) then
1153 if(grib=='grib2') then
1154 cfld=cfld+1
1155 fld_info(cfld)%ifld=iavblfld(iget(576))
1156 fld_info(cfld)%lvl=lvlsxml(ifd,iget(576))
1157!$omp parallel do private(i,j,ii,jj)
1158 do j=1,jend-jsta+1
1159 jj = jsta+j-1
1160 do i=1,iend-ista+1
1161 ii = ista+i-1
1162 datapd(i,j,cfld) = grid1(ii,jj)
1163 enddo
1164 enddo
1165 endif
1166 ENDIF
1167 ENDIF
1168 IF (iget(577)>0) THEN
1169 IF (lvls(ifd,iget(577))>0) THEN
1170 if(grib=='grib2') then
1171 cfld=cfld+1
1172 fld_info(cfld)%ifld=iavblfld(iget(577))
1173 fld_info(cfld)%lvl=lvlsxml(ifd,iget(577))
1174!$omp parallel do private(i,j,ii,jj)
1175 do j=1,jend-jsta+1
1176 jj = jsta+j-1
1177 do i=1,iend-ista+1
1178 ii = ista+i-1
1179 datapd(i,j,cfld) = grid2(ii,jj)
1180 enddo
1181 enddo
1182 endif
1183 ENDIF
1184 ENDIF
1185 ENDIF
1186
1187 END DO loop_10
1188 DEALLOCATE(t7d,q7d,u7d,v6d,p7d,icingfd)
1189 ENDIF
1190
1191!
1192! ***BLOCK 3-2: FD LEVEL (from control file) GTG
1193!
1194 IF(iget(467)>0.or.iget(468)>0.or.iget(469)>0) THEN
1195 if(iget(467)>0) THEN ! GTG
1196 n=iavblfld(iget(467))
1197 nfdctl=size(pset%param(n)%level)
1198 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
1199 allocate(itypefdlvlctl(nfdctl))
1200 DO ifd = 1,nfdctl
1201 itypefdlvlctl(ifd)=lvls(ifd,iget(467))
1202 enddo
1203 if(allocated(htfdctl)) deallocate(htfdctl)
1204 allocate(htfdctl(nfdctl))
1205 htfdctl=pset%param(n)%level
1206! print *, "GTG 467 levels=",pset%param(N)%level
1207 allocate(gtgfd(ista:iend,jsta:jend,nfdctl))
1208 call fdlvl_mass(itypefdlvlctl,nfdctl,htfdctl,gtg,gtgfd)
1209! print *, "GTG 467 Done GTGFD=",me,GTGFD(IM/2,jend,1:NFDCTL)
1210
1211 ! Regional GTG has a legend of special defination
1212 ! 0 m holds the max value of the whole vertical column
1213 DO ifd = 1,nfdctl
1214 if(nint(htfdctl(ifd)) == 0) then
1215 n=ifd
1216 exit
1217 endif
1218 ENDDO
1219 DO ifd = 1,nfdctl
1220 DO j=jsta,jend
1221 DO i=ista,iend
1222 work1=gtgfd(i,j,ifd)
1223 if(gtgfd(i,j,n)>=spval) then
1224 gtgfd(i,j,n)=work1
1225 elseif(work1<spval) then
1226 if(gtgfd(i,j,n)<work1) gtgfd(i,j,n)=work1
1227 endif
1228 ENDDO
1229 ENDDO
1230 ENDDO
1231
1232 DO ifd = 1,nfdctl
1233 IF (lvls(ifd,iget(467))>0) THEN
1234!$omp parallel do private(i,j)
1235 DO j=jsta,jend
1236 DO i=ista,iend
1237 grid1(i,j)=gtgfd(i,j,ifd)
1238 ENDDO
1239 ENDDO
1240 if(grib=='grib2') then
1241 cfld=cfld+1
1242 fld_info(cfld)%ifld=iavblfld(iget(467))
1243 fld_info(cfld)%lvl=lvlsxml(ifd,iget(467))
1244!$omp parallel do private(i,j,ii,jj)
1245 do j=1,jend-jsta+1
1246 jj = jsta+j-1
1247 do i=1,iend-ista+1
1248 ii = ista+i-1
1249 datapd(i,j,cfld) = grid1(ii,jj)
1250 enddo
1251 enddo
1252 endif
1253 ENDIF
1254 ENDDO
1255 endif
1256
1257 if(iget(468)>0) THEN ! CAT
1258 n=iavblfld(iget(468))
1259 nfdctl=size(pset%param(n)%level)
1260 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
1261 allocate(itypefdlvlctl(nfdctl))
1262 DO ifd = 1,nfdctl
1263 itypefdlvlctl(ifd)=lvls(ifd,iget(468))
1264 enddo
1265 if(allocated(htfdctl)) deallocate(htfdctl)
1266 allocate(htfdctl(nfdctl))
1267 htfdctl=pset%param(n)%level
1268 allocate(catfd(ista:iend,jsta:jend,nfdctl))
1269 call fdlvl_mass(itypefdlvlctl,nfdctl,htfdctl,catedr,catfd)
1270 DO ifd = 1,nfdctl
1271 IF (lvls(ifd,iget(468))>0) THEN
1272!$omp parallel do private(i,j)
1273 DO j=jsta,jend
1274 DO i=ista,iend
1275 grid1(i,j)=catfd(i,j,ifd)
1276 ENDDO
1277 ENDDO
1278 if(grib=='grib2') then
1279 cfld=cfld+1
1280 fld_info(cfld)%ifld=iavblfld(iget(468))
1281 fld_info(cfld)%lvl=lvlsxml(ifd,iget(468))
1282!$omp parallel do private(i,j,ii,jj)
1283 do j=1,jend-jsta+1
1284 jj = jsta+j-1
1285 do i=1,iend-ista+1
1286 ii = ista+i-1
1287 datapd(i,j,cfld) = grid1(ii,jj)
1288 enddo
1289 enddo
1290 endif
1291 ENDIF
1292 ENDDO
1293 endif
1294
1295 if(iget(469)>0) THEN ! MWT
1296 n=iavblfld(iget(469))
1297 nfdctl=size(pset%param(n)%level)
1298 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
1299 allocate(itypefdlvlctl(nfdctl))
1300 DO ifd = 1,nfdctl
1301 itypefdlvlctl(ifd)=lvls(ifd,iget(469))
1302 enddo
1303 if(allocated(htfdctl)) deallocate(htfdctl)
1304 allocate(htfdctl(nfdctl))
1305 htfdctl=pset%param(n)%level
1306 allocate(mwtfd(ista:iend,jsta:jend,nfdctl))
1307 call fdlvl_mass(itypefdlvlctl,nfdctl,htfdctl,mwt,mwtfd)
1308 DO ifd = 1,nfdctl
1309 IF (lvls(ifd,iget(469))>0) THEN
1310!$omp parallel do private(i,j)
1311 DO j=jsta,jend
1312 DO i=ista,iend
1313 grid1(i,j)=mwtfd(i,j,ifd)
1314 ENDDO
1315 ENDDO
1316 if(grib=='grib2') then
1317 cfld=cfld+1
1318 fld_info(cfld)%ifld=iavblfld(iget(469))
1319 fld_info(cfld)%lvl=lvlsxml(ifd,iget(469))
1320!$omp parallel do private(i,j,ii,jj)
1321 do j=1,jend-jsta+1
1322 jj = jsta+j-1
1323 do i=1,iend-ista+1
1324 ii = ista+i-1
1325 datapd(i,j,cfld) = grid1(ii,jj)
1326 enddo
1327 enddo
1328 endif
1329 ENDIF
1330 ENDDO
1331 endif
1332
1333 if(allocated(gtgfd)) deallocate(gtgfd)
1334 if(allocated(catfd)) deallocate(catfd)
1335 if(allocated(mwtfd)) deallocate(mwtfd)
1336
1337 if(allocated(itypefdlvlctl)) deallocate(itypefdlvlctl)
1338 if(allocated(htfdctl)) deallocate(htfdctl)
1339
1340 ENDIF
1341!
1342!
1343! ***BLOCK 4: FREEZING LEVEL Z, RH AND P.
1344!
1345 IF ( (iget(062)>0).OR.(iget(063)>0) ) THEN
1346 CALL frzlvl(z1d,rh1d,p1d)
1347!
1348! FREEZING LEVEL HEIGHT.
1349 IF (iget(062)>0) THEN
1350!$omp parallel do private(i,j)
1351 DO j=jsta,jend
1352 DO i=ista,iend
1353 grid1(i,j)=z1d(i,j)
1354 IF (submodelname == 'RTMA') THEN
1355 freezelvl(i,j)=grid1(i,j)
1356 ENDIF
1357 ENDDO
1358 ENDDO
1359 CALL bound (grid1,d00,h99999)
1360 if(grib=='grib2') then
1361 cfld=cfld+1
1362 fld_info(cfld)%ifld=iavblfld(iget(062))
1363!$omp parallel do private(i,j,ii,jj)
1364 do j=1,jend-jsta+1
1365 jj = jsta+j-1
1366 do i=1,iend-ista+1
1367 ii = ista+i-1
1368 datapd(i,j,cfld) = grid1(ii,jj)
1369 enddo
1370 enddo
1371 endif
1372 ENDIF
1373!
1374! FREEZING LEVEL RELATIVE HUMIDITY.
1375 IF (iget(063)>0) THEN
1376!$omp parallel do private(i,j)
1377 DO j=jsta,jend
1378 DO i=ista,iend
1379 grid1(i,j) = rh1d(i,j)
1380 ENDDO
1381 ENDDO
1382 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
1383 CALL bound(grid1,h1,h100)
1384 if(grib=='grib2') then
1385 cfld=cfld+1
1386 fld_info(cfld)%ifld=iavblfld(iget(063))
1387!$omp parallel do private(i,j,ii,jj)
1388 do j=1,jend-jsta+1
1389 jj = jsta+j-1
1390 do i=1,iend-ista+1
1391 ii = ista+i-1
1392 datapd(i,j,cfld) = grid1(ii,jj)
1393 enddo
1394 enddo
1395 endif
1396 ENDIF
1397!
1398! FREEZING LEVEL PRESSURE
1399 IF (iget(753)>0) THEN
1400!$omp parallel do private(i,j)
1401 DO j=jsta,jend
1402 DO i=ista,iend
1403 grid1(i,j) = p1d(i,j)
1404 ENDDO
1405 ENDDO
1406 if(grib=='grib2') then
1407 cfld=cfld+1
1408 fld_info(cfld)%ifld=iavblfld(iget(753))
1409!$omp parallel do private(i,j,ii,jj)
1410 do j=1,jend-jsta+1
1411 jj = jsta+j-1
1412 do i=1,iend-ista+1
1413 ii = ista+i-1
1414 datapd(i,j,cfld) = grid1(ii,jj)
1415 enddo
1416 enddo
1417 endif
1418
1419 ENDIF
1420 ENDIF
1421
1422 IF (iget(165)>0 .OR. iget(350)>0.OR. iget(756)>0) THEN
1423 CALL frzlvl2(tfrz,z1d,rh1d,p1d)
1424!
1425! HIGHEST FREEZING LEVEL HEIGHT.
1426 IF (iget(165)>0)THEN
1427!$omp parallel do private(i,j)
1428 DO j=jsta,jend
1429 DO i=ista,iend
1430 grid1(i,j)=z1d(i,j)
1431 ENDDO
1432 ENDDO
1433 CALL bound (grid1,d00,h99999)
1434 if(grib=='grib2') then
1435 cfld=cfld+1
1436 fld_info(cfld)%ifld=iavblfld(iget(165))
1437!$omp parallel do private(i,j,ii,jj)
1438 do j=1,jend-jsta+1
1439 jj = jsta+j-1
1440 do i=1,iend-ista+1
1441 ii = ista+i-1
1442 datapd(i,j,cfld) = grid1(ii,jj)
1443 enddo
1444 enddo
1445 endif
1446 END IF
1447
1448! HIGHEST FREEZING LEVEL RELATIVE HUMIDITY
1449 IF (iget(350)>0)THEN
1450 grid1=spval
1451!$omp parallel do private(i,j)
1452 DO j=jsta,jend
1453 DO i=ista,iend
1454 IF(rh1d(i,j) < spval) grid1(i,j)=rh1d(i,j)*100.
1455 ENDDO
1456 ENDDO
1457 CALL bound (grid1,h1,h100)
1458 if(grib=='grib2') then
1459 cfld=cfld+1
1460 fld_info(cfld)%ifld=iavblfld(iget(350))
1461!$omp parallel do private(i,j,ii,jj)
1462 do j=1,jend-jsta+1
1463 jj = jsta+j-1
1464 do i=1,iend-ista+1
1465 ii = ista+i-1
1466 datapd(i,j,cfld) = grid1(ii,jj)
1467 enddo
1468 enddo
1469 endif
1470 END IF
1471
1472! HIGHEST FREEZING LEVEL PRESSURE
1473 IF (iget(756)>0) THEN
1474!$omp parallel do private(i,j)
1475 DO j=jsta,jend
1476 DO i=ista,iend
1477 grid1(i,j) = p1d(i,j)
1478 ENDDO
1479 ENDDO
1480 if(grib=='grib2') then
1481 cfld=cfld+1
1482 fld_info(cfld)%ifld=iavblfld(iget(756))
1483!$omp parallel do private(i,j,ii,jj)
1484 do j=1,jend-jsta+1
1485 jj = jsta+j-1
1486 do i=1,iend-ista+1
1487 ii = ista+i-1
1488 datapd(i,j,cfld) = grid1(ii,jj)
1489 enddo
1490 enddo
1491 endif
1492 ENDIF
1493
1494 ENDIF
1495!
1496! HIGHEST -10 C ISOTHERM VALUES
1497!
1498 IF (iget(776)>0 .OR. iget(777)>0.OR. iget(778)>0) THEN
1499 CALL frzlvl2(263.15,z1d,rh1d,p1d)
1500!
1501! HIGHEST -10C ISOTHERM LEVEL HEIGHT.
1502 IF (iget(776)>0)THEN
1503!$omp parallel do private(i,j)
1504 DO j=jsta,jend
1505 DO i=ista,iend
1506 grid1(i,j)=z1d(i,j)
1507 ENDDO
1508 ENDDO
1509 CALL bound (grid1,d00,h99999)
1510 if(grib=='grib2') then
1511 cfld=cfld+1
1512 fld_info(cfld)%ifld=iavblfld(iget(776))
1513!$omp parallel do private(i,j,ii,jj)
1514 do j=1,jend-jsta+1
1515 jj = jsta+j-1
1516 do i=1,iend-ista+1
1517 ii = ista+i-1
1518 datapd(i,j,cfld) = grid1(ii,jj)
1519 enddo
1520 enddo
1521 endif
1522 END IF
1523
1524! HIGHEST -10C ISOTHERM RELATIVE HUMIDITY
1525 IF (iget(777)>0)THEN
1526 grid1=spval
1527!$omp parallel do private(i,j)
1528 DO j=jsta,jend
1529 DO i=ista,iend
1530 IF(rh1d(i,j) < spval) grid1(i,j)=rh1d(i,j)*100.
1531 ENDDO
1532 ENDDO
1533 CALL bound (grid1,h1,h100)
1534 if(grib=='grib2') then
1535 cfld=cfld+1
1536 fld_info(cfld)%ifld=iavblfld(iget(777))
1537!$omp parallel do private(i,j,ii,jj)
1538 do j=1,jend-jsta+1
1539 jj = jsta+j-1
1540 do i=1,iend-ista+1
1541 ii = ista+i-1
1542 datapd(i,j,cfld) = grid1(ii,jj)
1543 enddo
1544 enddo
1545 endif
1546 END IF
1547
1548! HIGHEST -10C ISOTHERM LEVEL PRESSURE
1549 IF (iget(778)>0) THEN
1550!$omp parallel do private(i,j)
1551 DO j=jsta,jend
1552 DO i=ista,iend
1553 grid1(i,j)=p1d(i,j)
1554 ENDDO
1555 ENDDO
1556 if(grib=='grib2') then
1557 cfld=cfld+1
1558 fld_info(cfld)%ifld=iavblfld(iget(778))
1559!$omp parallel do private(i,j,ii,jj)
1560 do j=1,jend-jsta+1
1561 jj = jsta+j-1
1562 do i=1,iend-ista+1
1563 ii = ista+i-1
1564 datapd(i,j,cfld) = grid1(ii,jj)
1565 enddo
1566 enddo
1567 endif
1568 ENDIF
1569
1570 ENDIF
1571!
1572! HIGHEST -20 C ISOTHERM VALUES
1573!
1574 IF (iget(779)>0 .OR. iget(780)>0.OR. iget(781)>0) THEN
1575 CALL frzlvl2(253.15,z1d,rh1d,p1d)
1576!
1577! HIGHEST -20C ISOTHERM LEVEL HEIGHT.
1578 IF (iget(779)>0)THEN
1579!$omp parallel do private(i,j)
1580 DO j=jsta,jend
1581 DO i=ista,iend
1582 grid1(i,j)=z1d(i,j)
1583 ENDDO
1584 ENDDO
1585 CALL bound (grid1,d00,h99999)
1586 if(grib=='grib2') then
1587 cfld=cfld+1
1588 fld_info(cfld)%ifld=iavblfld(iget(779))
1589!$omp parallel do private(i,j,ii,jj)
1590 do j=1,jend-jsta+1
1591 jj = jsta+j-1
1592 do i=1,iend-ista+1
1593 ii = ista+i-1
1594 datapd(i,j,cfld) = grid1(ii,jj)
1595 enddo
1596 enddo
1597 endif
1598 END IF
1599
1600! HIGHEST -20C ISOTHERM RELATIVE HUMIDITY
1601 IF (iget(780)>0)THEN
1602 grid1=spval
1603!$omp parallel do private(i,j)
1604 DO j=jsta,jend
1605 DO i=ista,iend
1606 IF(rh1d(i,j) < spval) grid1(i,j)=rh1d(i,j)*100.
1607 ENDDO
1608 ENDDO
1609 CALL bound (grid1,h1,h100)
1610 if(grib=='grib2') then
1611 cfld=cfld+1
1612 fld_info(cfld)%ifld=iavblfld(iget(780))
1613!$omp parallel do private(i,j,ii,jj)
1614 do j=1,jend-jsta+1
1615 jj = jsta+j-1
1616 do i=1,iend-ista+1
1617 ii = ista+i-1
1618 datapd(i,j,cfld) = grid1(ii,jj)
1619 enddo
1620 enddo
1621 endif
1622 END IF
1623
1624! HIGHEST -20C ISOTHERM LEVEL PRESSURE
1625 IF (iget(781)>0) THEN
1626!$omp parallel do private(i,j)
1627 DO j=jsta,jend
1628 DO i=ista,iend
1629 grid1(i,j)=p1d(i,j)
1630 ENDDO
1631 ENDDO
1632 if(grib=='grib2') then
1633 cfld=cfld+1
1634 fld_info(cfld)%ifld=iavblfld(iget(781))
1635!$omp parallel do private(i,j,ii,jj)
1636 do j=1,jend-jsta+1
1637 jj = jsta+j-1
1638 do i=1,iend-ista+1
1639 ii = ista+i-1
1640 datapd(i,j,cfld) = grid1(ii,jj)
1641 enddo
1642 enddo
1643 endif
1644 ENDIF
1645
1646 ENDIF
1647!
1648 allocate(pbnd(ista:iend,jsta:jend,nbnd), tbnd(ista:iend,jsta:jend,nbnd), &
1649 qbnd(ista:iend,jsta:jend,nbnd), ubnd(ista:iend,jsta:jend,nbnd), &
1650 vbnd(ista:iend,jsta:jend,nbnd), rhbnd(ista:iend,jsta:jend,nbnd), &
1651 wbnd(ista:iend,jsta:jend,nbnd))
1652
1653!
1654! ***BLOCK 5: BOUNDARY LAYER FIELDS.
1655!
1656 IF ( (iget(067)>0).OR.(iget(068)>0).OR. &
1657 (iget(069)>0).OR.(iget(070)>0).OR. &
1658 (iget(071)>0).OR.(iget(072)>0).OR. &
1659 (iget(073)>0).OR.(iget(074)>0).OR. &
1660 (iget(088)>0).OR.(iget(089)>0).OR. &
1661 (iget(090)>0).OR.(iget(075)>0).OR. &
1662 (iget(109)>0).OR.(iget(110)>0).OR. &
1663 (iget(031)>0).OR.(iget(032)>0).OR. &
1664 (iget(573)>0).OR. need_ifi .OR. &
1665 (iget(107)>0).OR.(iget(091)>0).OR. &
1666 (iget(092)>0).OR.(iget(093)>0).OR. &
1667 (iget(094)>0).OR.(iget(095)>0).OR. &
1668 (iget(096)>0).OR.(iget(097)>0).OR. &
1669 (iget(098)>0).OR.(iget(221)>0) ) THEN
1670!
1671 call allocate_cape_arrays
1672! COMPUTE ETA BOUNDARY LAYER FIELDS.
1673 CALL bndlyr(pbnd,tbnd,qbnd,rhbnd,ubnd,vbnd, &
1674 wbnd,omgbnd,pwtbnd,qcnvbnd,lvlbnd)
1675
1676!$omp parallel do private(i,j)
1677 DO j=jsta,jend
1678 DO i=ista,iend
1679 egrid2(i,j) = spval
1680 ENDDO
1681 ENDDO
1682
1683!
1684! LOOP OVER NBND BOUNDARY LAYERS.
1685 boundary_layer_loop: DO lbnd = 1,nbnd
1686!
1687! BOUNDARY LAYER PRESSURE.
1688 IF (iget(067)>0) THEN
1689 IF (lvls(lbnd,iget(067))>0) THEN
1690!$omp parallel do private(i,j)
1691 DO j=jsta,jend
1692 DO i=ista,iend
1693 grid1(i,j) = pbnd(i,j,lbnd)
1694 ENDDO
1695 ENDDO
1696 if(grib=='grib2') then
1697 cfld=cfld+1
1698 fld_info(cfld)%ifld=iavblfld(iget(067))
1699 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(067))
1700!$omp parallel do private(i,j,ii,jj)
1701 do j=1,jend-jsta+1
1702 jj = jsta+j-1
1703 do i=1,iend-ista+1
1704 ii = ista+i-1
1705 datapd(i,j,cfld) = grid1(ii,jj)
1706 enddo
1707 enddo
1708 endif
1709 ENDIF
1710 ENDIF
1711!
1712! BOUNDARY LAYER TEMPERATURE.
1713 IF (iget(068)>0) THEN
1714 IF (lvls(lbnd,iget(068))>0) THEN
1715!$omp parallel do private(i,j)
1716 DO j=jsta,jend
1717 DO i=ista,iend
1718 grid1(i,j)=tbnd(i,j,lbnd)
1719 ENDDO
1720 ENDDO
1721 if(grib=='grib2') then
1722 cfld=cfld+1
1723 fld_info(cfld)%ifld=iavblfld(iget(068))
1724 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(068))
1725!$omp parallel do private(i,j,ii,jj)
1726 do j=1,jend-jsta+1
1727 jj = jsta+j-1
1728 do i=1,iend-ista+1
1729 ii = ista+i-1
1730 datapd(i,j,cfld) = grid1(ii,jj)
1731 enddo
1732 enddo
1733 endif
1734 ENDIF
1735 ENDIF
1736!
1737! BOUNDARY LAYER POTENTIAL TEMPERATURE.
1738 IF (iget(069)>0) THEN
1739 IF (lvls(lbnd,iget(069))>0) THEN
1740 CALL calpot(pbnd(ista,jsta,lbnd),tbnd(ista,jsta,lbnd),grid1(ista:iend,jsta:jend))
1741 if(grib=='grib2') then
1742 cfld=cfld+1
1743 fld_info(cfld)%ifld=iavblfld(iget(069))
1744 fld_info(cfld)%lvl=lvlsxml(ifd,iget(069))
1745!$omp parallel do private(i,j,ii,jj)
1746 do j=1,jend-jsta+1
1747 jj = jsta+j-1
1748 do i=1,iend-ista+1
1749 ii = ista+i-1
1750 datapd(i,j,cfld) = grid1(ii,jj)
1751 enddo
1752 enddo
1753 endif
1754 ENDIF
1755 ENDIF
1756!
1757! BOUNDARY LAYER RELATIVE HUMIDITY.
1758 IF (iget(072)>0) THEN
1759 IF (lvls(lbnd,iget(072))>0) THEN
1760!$omp parallel do private(i,j)
1761 DO j=jsta,jend
1762 DO i=ista,iend
1763 grid1(i,j)=rhbnd(i,j,lbnd)
1764 ENDDO
1765 ENDDO
1766 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
1767 CALL bound(grid1,h1,h100)
1768 if(grib=='grib2') then
1769 cfld=cfld+1
1770 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(072))
1771 fld_info(cfld)%ifld=iavblfld(iget(072))
1772!$omp parallel do private(i,j,ii,jj)
1773 do j=1,jend-jsta+1
1774 jj = jsta+j-1
1775 do i=1,iend-ista+1
1776 ii = ista+i-1
1777 datapd(i,j,cfld) = grid1(ii,jj)
1778 enddo
1779 enddo
1780 endif
1781 ENDIF
1782 ENDIF
1783!
1784! BOUNDARY LAYER DEWPOINT TEMPERATURE.
1785 IF (iget(070)>0) THEN
1786 IF (lvls(lbnd,iget(070))>0) THEN
1787 CALL caldwp(pbnd(ista:iend,jsta:jend,lbnd), qbnd(ista:iend,jsta:jend,lbnd), &
1788 grid1(ista:iend,jsta:jend), tbnd(ista:iend,jsta:jend,lbnd))
1789 if(grib=='grib2') then
1790 cfld=cfld+1
1791 fld_info(cfld)%ifld=iavblfld(iget(070))
1792 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(070))
1793!$omp parallel do private(i,j,ii,jj)
1794 do j=1,jend-jsta+1
1795 jj = jsta+j-1
1796 do i=1,iend-ista+1
1797 ii = ista+i-1
1798 datapd(i,j,cfld) = grid1(ii,jj)
1799 enddo
1800 enddo
1801 endif
1802 ENDIF
1803 ENDIF
1804!
1805! BOUNDARY LAYER SPECIFIC HUMIDITY.
1806 IF (iget(071)>0) THEN
1807 IF (lvls(lbnd,iget(071))>0) THEN
1808!$omp parallel do private(i,j)
1809 DO j=jsta,jend
1810 DO i=ista,iend
1811 grid1(i,j)=qbnd(i,j,lbnd)
1812 ENDDO
1813 ENDDO
1814 CALL bound(grid1,h1m12,h99999)
1815 if(grib=='grib2') then
1816 cfld=cfld+1
1817 fld_info(cfld)%ifld=iavblfld(iget(071))
1818 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(071))
1819!$omp parallel do private(i,j,ii,jj)
1820 do j=1,jend-jsta+1
1821 jj = jsta+j-1
1822 do i=1,iend-ista+1
1823 ii = ista+i-1
1824 datapd(i,j,cfld) = grid1(ii,jj)
1825 enddo
1826 enddo
1827 endif
1828 ENDIF
1829 ENDIF
1830!
1831! BOUNDARY LAYER MOISTURE CONVERGENCE.
1832 IF (iget(088)>0) THEN
1833 IF (lvls(lbnd,iget(088))>0) THEN
1834!$omp parallel do private(i,j)
1835 DO j=jsta,jend
1836 DO i=ista,iend
1837 grid1(i,j) = qcnvbnd(i,j,lbnd)
1838 ENDDO
1839 ENDDO
1840 if(grib=='grib2') then
1841 cfld=cfld+1
1842 fld_info(cfld)%ifld=iavblfld(iget(088))
1843 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(088))
1844!$omp parallel do private(i,j,ii,jj)
1845 do j=1,jend-jsta+1
1846 jj = jsta+j-1
1847 do i=1,iend-ista+1
1848 ii = ista+i-1
1849 datapd(i,j,cfld) = grid1(ii,jj)
1850 enddo
1851 enddo
1852 endif
1853 ENDIF
1854 ENDIF
1855!
1856! BOUNDARY LAYER U WIND AND/OR V WIND.
1857!
1858 field1=.false.
1859 field2=.false.
1860!
1861 IF(iget(073)>0)THEN
1862 IF(lvls(lbnd,iget(073))>0)field1=.true.
1863 ENDIF
1864 IF(iget(074)>0)THEN
1865 IF(lvls(lbnd,iget(074))>0)field2=.true.
1866 ENDIF
1867!
1868 IF(field1.OR.field2)THEN
1869!$omp parallel do private(i,j)
1870 DO j=jsta,jend
1871 DO i=ista,iend
1872 grid1(i,j) = ubnd(i,j,lbnd)
1873 grid2(i,j) = vbnd(i,j,lbnd)
1874 ENDDO
1875 ENDDO
1876!
1877 IF (iget(073)>0) THEN
1878 IF (lvls(lbnd,iget(073))>0) then
1879 if(grib=='grib2') then
1880 cfld=cfld+1
1881 fld_info(cfld)%ifld=iavblfld(iget(073))
1882 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(073))
1883!$omp parallel do private(i,j,ii,jj)
1884 do j=1,jend-jsta+1
1885 jj = jsta+j-1
1886 do i=1,iend-ista+1
1887 ii = ista+i-1
1888 datapd(i,j,cfld) = grid1(ii,jj)
1889 enddo
1890 enddo
1891 endif
1892 ENDIF
1893 ENDIF
1894 IF (iget(074)>0) THEN
1895 IF (lvls(lbnd,iget(074))>0) THEN
1896 if(grib=='grib2') then
1897 cfld=cfld+1
1898 fld_info(cfld)%ifld=iavblfld(iget(074))
1899 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(074))
1900!$omp parallel do private(i,j,ii,jj)
1901 do j=1,jend-jsta+1
1902 jj = jsta+j-1
1903 do i=1,iend-ista+1
1904 ii = ista+i-1
1905 datapd(i,j,cfld) = grid2(ii,jj)
1906 enddo
1907 enddo
1908 endif
1909 ENDIF
1910 ENDIF
1911 ENDIF
1912!
1913! BOUNDARY LAYER OMEGA.
1914 IF (iget(090)>0) THEN
1915 IF (lvls(lbnd,iget(090))>0) THEN
1916!$omp parallel do private(i,j)
1917 DO j=jsta,jend
1918 DO i=ista,iend
1919 grid1(i,j) = omgbnd(i,j,lbnd)
1920 ENDDO
1921 ENDDO
1922 if(grib=='grib2') then
1923 cfld=cfld+1
1924 fld_info(cfld)%ifld=iavblfld(iget(090))
1925 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(090))
1926!$omp parallel do private(i,j,ii,jj)
1927 do j=1,jend-jsta+1
1928 jj = jsta+j-1
1929 do i=1,iend-ista+1
1930 ii = ista+i-1
1931 datapd(i,j,cfld) = grid1(ii,jj)
1932 enddo
1933 enddo
1934 endiF
1935 ENDIF
1936 ENDIF
1937!
1938! BOUNDARY LAYER PRECIPITBLE WATER.
1939 IF (iget(089)>0) THEN
1940 IF (lvls(lbnd,iget(089))>0) THEN
1941!$omp parallel do private(i,j)
1942 DO j=jsta,jend
1943 DO i=ista,iend
1944 grid1(i,j) = pwtbnd(i,j,lbnd)
1945 ENDDO
1946 ENDDO
1947 CALL bound(grid1,d00,h99999)
1948 if(grib=='grib2') then
1949 cfld=cfld+1
1950 fld_info(cfld)%ifld=iavblfld(iget(089))
1951 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(089))
1952!$omp parallel do private(i,j,ii,jj)
1953 do j=1,jend-jsta+1
1954 jj = jsta+j-1
1955 do i=1,iend-ista+1
1956 ii = ista+i-1
1957 datapd(i,j,cfld) = grid1(ii,jj)
1958 enddo
1959 enddo
1960 endif
1961 ENDIF
1962 ENDIF
1963!
1964! BOUNDARY LAYER LIFTED INDEX.
1965 IF (iget(075)>0 .OR. iget(031)>0 .OR. iget(573)>0) THEN
1966 CALL otlft(pbnd(ista,jsta,lbnd),tbnd(ista,jsta,lbnd), &
1967 qbnd(ista,jsta,lbnd),grid1(ista:iend,jsta:jend))
1968 IF(iget(075)>0)THEN
1969 IF (lvls(lbnd,iget(075))>0) THEN
1970 if(grib=='grib2') then
1971 cfld=cfld+1
1972 fld_info(cfld)%ifld=iavblfld(iget(075))
1973 fld_info(cfld)%lvl=lvlsxml(lbnd,iget(075))
1974!$omp parallel do private(i,j,ii,jj)
1975 do j=1,jend-jsta+1
1976 jj = jsta+j-1
1977 do i=1,iend-ista+1
1978 ii = ista+i-1
1979 datapd(i,j,cfld) = grid1(ii,jj)
1980 enddo
1981 enddo
1982 endif
1983 ENDIF
1984 END IF
1985 IF(iget(031)>0 .or. iget(573)>0)THEN
1986!$omp parallel do private(i,j)
1987 DO j=jsta,jend
1988 DO i=ista,iend
1989 egrid2(i,j) = min(egrid2(i,j),grid1(i,j))
1990 END DO
1991 END DO
1992 END IF
1993 ENDIF
1994!
1995! END OF ETA BOUNDARY LAYER LOOP.
1996 END DO boundary_layer_loop
1997 deallocate(omgbnd,pwtbnd,qcnvbnd)
1998!
1999! BEST LIFTED INDEX FROM BOUNDARY LAYER FIELDS.
2000!
2001 IF (iget(031)>0 .OR. iget(573)>0 ) THEN
2002! DO J=JSTA,JEND
2003! DO I=ISTA,IEND
2004! EGRID1(I,J) = H99999
2005! EGRID2(I,J) = H99999
2006! ENDDO
2007! ENDDO
2008!
2009! DO 50 LBND = 1,NBND
2010! CALL OTLFT(PBND(1,1,LBND),TBND(1,1,LBND), &
2011! QBND(1,1,LBND),EGRID2)
2012! DO J=JSTA,JEND
2013! DO I=ISTA,IEND
2014! EGRID1(I,J)=AMIN1(EGRID1(I,J),EGRID2(I,J))
2015! ENDDO
2016! ENDDO
2017! 50 CONTINUE
2018!$omp parallel do private(i,j)
2019 DO j=jsta,jend
2020 DO i=ista,iend
2021 grid1(i,j)=egrid2(i,j)
2022 ENDDO
2023 ENDDO
2024! print*,'writting out best lifted index'
2025
2026 if (iget(031)>0) then
2027 if(grib=='grib2') then
2028 cfld=cfld+1
2029 fld_info(cfld)%ifld=iavblfld(iget(031))
2030 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2031 endif
2032 endif
2033
2034 if(iget(573)> 0 ) THEN
2035 if(grib=='grib2') then
2036 cfld=cfld+1
2037 fld_info(cfld)%ifld=iavblfld(iget(573))
2038!$omp parallel do private(i,j,ii,jj)
2039 do j=1,jend-jsta+1
2040 jj = jsta+j-1
2041 do i=1,iend-ista+1
2042 ii = ista+i-1
2043 datapd(i,j,cfld) = grid1(ii,jj)
2044 enddo
2045 enddo
2046 endif
2047 endif
2048
2049 END IF
2050!
2051! BEST BOUNDARY LAYER CAPE AND CINS.
2052!
2053 field1=.false.
2054 field2=.false.
2055!
2056 IF(iget(032)>0)THEN
2057 IF(lvls(2,iget(032))>0)field1=.true.
2058 ENDIF
2059 IF(iget(107)>0)THEN
2060 IF(lvls(2,iget(107))>0)field2=.true.
2061 ENDIF
2062!
2063 IF(iget(566)>0)THEN
2064 field1=.true.
2065 ENDIF
2066 IF(iget(567)>0)THEN
2067 field2=.true.
2068 ENDIF
2069!
2070 !if(grib=="grib2") print *,'in MISCLN.f,iget(566)=', &
2071 ! iget(566), 'iget(567)=',iget(567),'LVLSXML(1,IGET(566)=', &
2072 ! LVLSXML(1,IGET(566)),'LVLSXML(1,IGET(567)=', &
2073 ! LVLSXML(1,IGET(567)),'field1=',field1,'field2=',field2
2074!
2075 IF(field1.OR.field2.OR.need_ifi)THEN
2076 itype = 2
2077 call allocate_cape_arrays
2078!
2079!$omp parallel do private(i,j)
2080 DO j=jsta,jend
2081 DO i=ista,iend
2082 egrid1(i,j) = -h99999
2083 egrid2(i,j) = -h99999
2084 ENDDO
2085 ENDDO
2086!
2087 loop_80: DO lbnd = 1,nbnd
2088 CALL calthte(pbnd(ista,jsta,lbnd),tbnd(ista,jsta,lbnd), &
2089 qbnd(ista,jsta,lbnd),egrid1)
2090!$omp parallel do private(i,j)
2091 DO j=jsta,jend
2092 DO i=ista,iend
2093 IF (egrid1(i,j) > egrid2(i,j)) THEN
2094 egrid2(i,j) = egrid1(i,j)
2095 lb2(i,j) = lvlbnd(i,j,lbnd)
2096 p1d(i,j) = pbnd(i,j,lbnd)
2097 t1d(i,j) = tbnd(i,j,lbnd)
2098 q1d(i,j) = qbnd(i,j,lbnd)
2099 ENDIF
2100 ENDDO
2101 ENDDO
2102 ENDDO loop_80
2103!
2104 dpbnd = 0.
2105 CALL calcape(itype,dpbnd,p1d,t1d,q1d,lb2,egrid1, &
2106 egrid2,egrid3,egrid4,egrid5)
2107
2108!
2109 IF(iget(566)>0 .or. need_ifi) THEN
2110 grid1=spval
2111!$omp parallel do private(i,j)
2112 DO j=jsta,jend
2113 DO i=ista,iend
2114 IF(t1d(i,j) < spval) grid1(i,j) = egrid1(i,j)
2115 ENDDO
2116 ENDDO
2117 CALL bound(grid1,d00,h99999)
2118!$omp parallel do private(i,j)
2119 DO j=jsta,jend
2120 DO i=ista,iend
2121 cape(i,j) = grid1(i,j)
2122 ENDDO
2123 ENDDO
2124 ENDIF
2125
2126 IF (iget(566)>0) THEN
2127 if(grib=='grib2') then
2128 cfld=cfld+1
2129 fld_info(cfld)%ifld=iavblfld(iget(566))
2130 fld_info(cfld)%lvl=lvlsxml(1,iget(566))
2131!$omp parallel do private(i,j,ii,jj)
2132 do j=1,jend-jsta+1
2133 jj = jsta+j-1
2134 do i=1,iend-ista+1
2135 ii = ista+i-1
2136 datapd(i,j,cfld) = cape(ii,jj)
2137 enddo
2138 enddo
2139 endif
2140 ENDIF
2141!
2142 IF (iget(567) > 0 .or. need_ifi) THEN
2143! dong add missing value for CIN
2144 grid1=spval
2145!$omp parallel do private(i,j)
2146 DO j=jsta,jend
2147 DO i=ista,iend
2148 IF(t1d(i,j) < spval) grid1(i,j) = - egrid2(i,j)
2149 ENDDO
2150 ENDDO
2151!
2152 CALL bound(grid1,d00,h99999)
2153!
2154!$omp parallel do private(i,j)
2155 DO j=jsta,jend
2156 DO i=ista,iend
2157 IF(t1d(i,j) < spval) grid1(i,j) = - grid1(i,j)
2158 cin(i,j) = grid1(i,j)
2159 ENDDO
2160 ENDDO
2161 ENDIF
2162
2163 IF(iget(567) > 0) THEN
2164 if(grib=='grib2') then
2165 cfld=cfld+1
2166 fld_info(cfld)%ifld=iavblfld(iget(567))
2167 fld_info(cfld)%lvl=lvlsxml(1,iget(567))
2168!$omp parallel do private(i,j,ii,jj)
2169 do j=1,jend-jsta+1
2170 jj = jsta+j-1
2171 do i=1,iend-ista+1
2172 ii = ista+i-1
2173 datapd(i,j,cfld) = cin(ii,jj)
2174 enddo
2175 enddo
2176 endif
2177 ENDIF
2178 ENDIF
2179!
2180
2181! PBL HEIGHT
2182 IF(iget(221) > 0) THEN
2183!$omp parallel do private(i,j)
2184 DO j=jsta,jend
2185 DO i=ista,iend
2186 grid1(i,j) = pblh(i,j)
2187 ENDDO
2188 ENDDO
2189 if(grib=='grib2') then
2190 cfld=cfld+1
2191 fld_info(cfld)%ifld=iavblfld(iget(221))
2192!$omp parallel do private(i,j,ii,jj)
2193 do j=1,jend-jsta+1
2194 jj = jsta+j-1
2195 do i=1,iend-ista+1
2196 ii = ista+i-1
2197 datapd(i,j,cfld) = grid1(ii,jj)
2198 enddo
2199 enddo
2200 endif
2201 END IF
2202! BOUNDARY LAYER LIFTING CONDENSATION PRESSURE AND HEIGHT.
2203! EGRID1 IS LCL PRESSURE. EGRID2 IS LCL HEIGHT.
2204!
2205 IF ( (iget(109)>0).OR.(iget(110)>0) ) THEN
2206 CALL callcl(pbnd(ista,jsta,1),tbnd(ista,jsta,1), &
2207 qbnd(ista,jsta,1),egrid1,egrid2)
2208 IF (iget(109)>0) THEN
2209 grid1=spval
2210!$omp parallel do private(i,j)
2211 DO j=jsta,jend
2212 DO i=ista,iend
2213 IF(tbnd(i,j,1) < spval) grid1(i,j) = egrid2(i,j)
2214 ENDDO
2215 ENDDO
2216 if(grib=='grib2') then
2217 cfld=cfld+1
2218 fld_info(cfld)%ifld=iavblfld(iget(109))
2219!$omp parallel do private(i,j,ii,jj)
2220 do j=1,jend-jsta+1
2221 jj = jsta+j-1
2222 do i=1,iend-ista+1
2223 ii = ista+i-1
2224 datapd(i,j,cfld) = grid1(ii,jj)
2225 enddo
2226 enddo
2227 endif
2228 ENDIF
2229 IF (iget(110)>0) THEN
2230 grid1=spval
2231!$omp parallel do private(i,j)
2232 DO j=jsta,jend
2233 DO i=ista,iend
2234 IF(tbnd(i,j,1) < spval) grid1(i,j) = egrid1(i,j)
2235 ENDDO
2236 ENDDO
2237 if(grib=='grib2') then
2238 cfld=cfld+1
2239 fld_info(cfld)%ifld=iavblfld(iget(110))
2240!$omp parallel do private(i,j,ii,jj)
2241 do j=1,jend-jsta+1
2242 jj = jsta+j-1
2243 do i=1,iend-ista+1
2244 ii = ista+i-1
2245 datapd(i,j,cfld) = grid1(ii,jj)
2246 enddo
2247 enddo
2248 endif
2249 ENDIF
2250 ENDIF
2251!
2252! NGM BOUNDARY LAYER FIELDS.
2253!
2254 IF ( (iget(091)>0).OR.(iget(092)>0).OR. &
2255 (iget(093)>0).OR.(iget(094)>0).OR. &
2256 (iget(095)>0).OR.(iget(095)>0).OR. &
2257 (iget(096)>0).OR.(iget(097)>0).OR. &
2258 (iget(098)>0) ) THEN
2259
2260 allocate(t78483(ista:iend,jsta:jend), t89671(ista:iend,jsta:jend), &
2261 p78483(ista:iend,jsta:jend), p89671(ista:iend,jsta:jend))
2262!
2263! COMPUTE SIGMA 0.89671 AND 0.78483 TEMPERATURES
2264! INTERPOLATE LINEAR IN LOG P
2265 IF (iget(097)>0.OR.iget(098)>0) THEN
2266!$omp parallel do private(i,j)
2267 DO j=jsta,jend
2268 DO i=ista,iend
2269 p78483(i,j) = log(pint(i,j,nint(lmh(i,j)))*0.78483)
2270 p89671(i,j) = log(pint(i,j,nint(lmh(i,j)))*0.89671)
2271 ENDDO
2272 ENDDO
2273 done =.false.
2274 done1=.false.
2275!!$omp parallel do private(fac1,fac2,pkl1,pku1,t78483,t89671)
2276 DO l=2,lm
2277 DO j=jsta,jend
2278 DO i=ista,iend
2279 pkl1=0.5*(alpint(i,j,l)+alpint(i,j,l+1))
2280 pku1=0.5*(alpint(i,j,l)+alpint(i,j,l-1))
2281! IF(I==1 .AND. J==1)PRINT*,'L,P89671,PKL1,PKU1= ', &
2282! L,P89671(I,J), PKL1, PKU1
2283 IF(p78483(i,j) < pkl1.AND.p78483(i,j) > pku1) THEN
2284 fac1 = (pkl1-p78483(i,j))/(pkl1-pku1)
2285 fac2 = (p78483(i,j)-pku1)/(pkl1-pku1)
2286 t78483(i,j) = t(i,j,l)*fac2 + t(i,j,l-1)*fac1
2287 done1(i,j)=.true.
2288 ENDIF
2289 IF(p89671(i,j) < pkl1.AND.p89671(i,j) > pku1)THEN
2290 fac1 = (pkl1-p89671(i,j))/(pkl1-pku1)
2291 fac2 = (p89671(i,j)-pku1)/(pkl1-pku1)
2292 t89671(i,j) = t(i,j,l)*fac2 + t(i,j,l-1)*fac1
2293 done(i,j) = .true.
2294! IF(I==1 .AND. J==1)PRINT*,'done(1,1)= ',done(1,1)
2295 ENDIF
2296 ENDDO
2297 ENDDO
2298 ENDDO
2299! print*,'done(1,1)= ',done(1,1)
2300!$omp parallel do private(i,j,pl,tl,ql,qsat,rhl)
2301 DO j=jsta,jend
2302 DO i=ista,iend
2303 IF(.NOT. done(i,j)) THEN
2304 pl = pint(i,j,lm-1)
2305 tl = 0.5*(t(i,j,lm-2)+t(i,j,lm-1))
2306 ql = 0.5*(q(i,j,lm-2)+q(i,j,lm-1))
2307 qsat = pq0/pl *exp(a2*(tl-a3)/(tl-a4))
2308!
2309 rhl = ql/qsat
2310!
2311 IF(rhl > 1.)THEN
2312 rhl = 1.
2313 ql = rhl*qsat
2314 ENDIF
2315!
2316 IF(rhl < rhmin)THEN
2317 rhl = rhmin
2318 ql = rhl*qsat
2319 ENDIF
2320!
2321! TVRL = TL*(1.+0.608*QL)
2322! TVRBLO = TVRL*(P89671(I,J)/PL)**RGAMOG
2323! T89671(I,J) = TVRBLO/(1.+0.608*QL)
2324
2325 t89671(i,j) = tl * (p89671(i,j)/pl)**rgamog
2326!
2327
2328! PKL1=0.5*(ALPINT(I,J,LM)+ALPINT(I,J,LM+1))
2329! PKU1=0.5*(ALPINT(I,J,LM-1)+ALPINT(I,J,LM))
2330! T89671(I,J)=T(I,J,LM)+(T(I,J,LM)-T(I,J,LM-1))*
2331! + (P89671(I,J)-PKL1)/(PKL1-PKU1)
2332
2333! print*,'Debug T89671= ',i,j
2334! + ,P89671(I,J), PKL1, PKU1
2335! + ,T89671(I,J),T(I,J,LM-1),T(I,J,LM)
2336 END IF
2337
2338 IF(.NOT. done1(i,j))THEN
2339 pl = pint(i,j,lm-1)
2340 tl = 0.5*(t(i,j,lm-2)+t(i,j,lm-1))
2341 ql = 0.5*(q(i,j,lm-2)+q(i,j,lm-1))
2342 qsat = pq0/pl *exp(a2*(tl-a3)/(tl-a4))
2343!
2344 rhl = ql/qsat
2345!
2346 IF(rhl>1.)THEN
2347 rhl = 1.
2348 ql = rhl*qsat
2349 ENDIF
2350!
2351 IF(rhl<rhmin)THEN
2352 rhl = rhmin
2353 ql = rhl*qsat
2354 ENDIF
2355!
2356! TVRL = TL*(1.+0.608*QL)
2357! TVRBLO = TVRL*(P78483(I,J)/PL)**RGAMOG
2358! T78483(I,J) =TVRBLO/(1.+0.608*QL)
2359
2360 t78483(i,j) = tl * (p78483(i,j)/pl)**rgamog
2361!
2362 END IF
2363
2364 END DO
2365 END DO
2366!
2367! SIGMA 0.89671 TEMPERATURE
2368 IF (iget(097) > 0) THEN
2369 grid1=spval
2370!$omp parallel do private(i,j)
2371 DO j=jsta,jend
2372 DO i=ista,iend
2373 IF(t(i,j,lm) < spval) grid1(i,j) = t89671(i,j)
2374! IF(T89671(I,J)>350.)PRINT*,'LARGE T89671 ', &
2375! I,J,T89671(I,J)
2376 ENDDO
2377 ENDDO
2378 if(grib=='grib2') then
2379 cfld=cfld+1
2380 fld_info(cfld)%ifld=iavblfld(iget(097))
2381 fld_info(cfld)%lvl=lvlsxml(1,iget(097))
2382!$omp parallel do private(i,j,ii,jj)
2383 do j=1,jend-jsta+1
2384 jj = jsta+j-1
2385 do i=1,iend-ista+1
2386 ii = ista+i-1
2387 datapd(i,j,cfld) = grid1(ii,jj)
2388 enddo
2389 enddo
2390 endif
2391 ENDIF
2392!
2393! SIGMA 0.78483 TEMPERATURE
2394 IF (iget(098)>0) THEN
2395 grid1=spval
2396!$omp parallel do private(i,j)
2397 DO j=jsta,jend
2398 DO i=ista,iend
2399 IF(t(i,j,lm) < spval) grid1(i,j) = t78483(i,j)
2400 ENDDO
2401 ENDDO
2402 if(grib=='grib2') then
2403 cfld=cfld+1
2404 fld_info(cfld)%ifld=iavblfld(iget(098))
2405 fld_info(cfld)%lvl=lvlsxml(1,iget(098))
2406!$omp parallel do private(i,j,ii,jj)
2407 do j=1,jend-jsta+1
2408 jj = jsta+j-1
2409 do i=1,iend-ista+1
2410 ii = ista+i-1
2411 datapd(i,j,cfld) = grid1(ii,jj)
2412 enddo
2413 enddo
2414 endif
2415 ENDIF
2416 deallocate(t78483, t89671, p78483, p89671)
2417 ENDIF
2418!
2419! NGM SIGMA LAYER 0.98230 FIELDS. THESE FIELDS ARE
2420! THE FIRST ETA LAYER BOUNDARY LAYER FIELDS.
2421!
2422!
2423 IF ( (iget(091)>0).OR.(iget(092)>0).OR. &
2424 (iget(093)>0).OR.(iget(094)>0).OR. &
2425 (iget(095)>0).OR.(iget(095)>0).OR. &
2426 (iget(096)>0) ) THEN
2427!
2428!
2429! PRESSURE.
2430 IF (iget(091)>0) THEN
2431!$omp parallel do private(i,j)
2432 DO j=jsta,jend
2433 DO i=ista,iend
2434 grid1(i,j) = pbnd(i,j,1)
2435 ENDDO
2436 ENDDO
2437 if(grib=='grib2') then
2438 cfld=cfld+1
2439 fld_info(cfld)%ifld=iavblfld(iget(091))
2440!$omp parallel do private(i,j,ii,jj)
2441 do j=1,jend-jsta+1
2442 jj = jsta+j-1
2443 do i=1,iend-ista+1
2444 ii = ista+i-1
2445 datapd(i,j,cfld) = grid1(ii,jj)
2446 enddo
2447 enddo
2448 endif
2449 ENDIF
2450!
2451! TEMPERATURE.
2452 IF (iget(092)>0) THEN
2453!$omp parallel do private(i,j)
2454 DO j=jsta,jend
2455 DO i=ista,iend
2456 grid1(i,j) = tbnd(i,j,1)
2457 ENDDO
2458 ENDDO
2459 if(grib=='grib2') then
2460 cfld=cfld+1
2461 fld_info(cfld)%ifld=iavblfld(iget(092))
2462 fld_info(cfld)%lvl=lvlsxml(1,iget(092))
2463!$omp parallel do private(i,j,ii,jj)
2464 do j=1,jend-jsta+1
2465 jj = jsta+j-1
2466 do i=1,iend-ista+1
2467 ii = ista+i-1
2468 datapd(i,j,cfld) = grid1(ii,jj)
2469 enddo
2470 enddo
2471 endif
2472 ENDIF
2473!
2474! SPECIFIC HUMIDITY.
2475 IF (iget(093)>0) THEN
2476!$omp parallel do private(i,j)
2477 DO j=jsta,jend
2478 DO i=ista,iend
2479 grid1(i,j) = qbnd(i,j,1)
2480 ENDDO
2481 ENDDO
2482 CALL bound(grid1,h1m12,h99999)
2483 if(grib=='grib2') then
2484 cfld=cfld+1
2485 fld_info(cfld)%ifld=iavblfld(iget(093))
2486 fld_info(cfld)%lvl=lvlsxml(1,iget(093))
2487!$omp parallel do private(i,j,ii,jj)
2488 do j=1,jend-jsta+1
2489 jj = jsta+j-1
2490 do i=1,iend-ista+1
2491 ii = ista+i-1
2492 datapd(i,j,cfld) = grid1(ii,jj)
2493 enddo
2494 enddo
2495 endif
2496 ENDIF
2497!
2498! RELATIVE HUMIDITY.
2499 IF (iget(094)>0) THEN
2500!$omp parallel do private(i,j)
2501 DO j=jsta,jend
2502 DO i=ista,iend
2503 grid1(i,j) = rhbnd(i,j,1)
2504 ENDDO
2505 ENDDO
2506 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2507 CALL bound(grid1,h1,h100)
2508 if(grib=='grib2') then
2509 cfld=cfld+1
2510 fld_info(cfld)%ifld=iavblfld(iget(094))
2511 fld_info(cfld)%lvl=lvlsxml(1,iget(094))
2512!$omp parallel do private(i,j,ii,jj)
2513 do j=1,jend-jsta+1
2514 jj = jsta+j-1
2515 do i=1,iend-ista+1
2516 ii = ista+i-1
2517 datapd(i,j,cfld) = grid1(ii,jj)
2518 enddo
2519 enddo
2520 endif
2521 ENDIF
2522!
2523! U AND/OR V WIND.
2524 IF ((iget(095)>0).OR.(iget(096)>0)) THEN
2525!$omp parallel do private(i,j)
2526 DO j=jsta,jend
2527 DO i=ista,iend
2528 grid1(i,j) = ubnd(i,j,1)
2529 grid2(i,j) = vbnd(i,j,1)
2530 ENDDO
2531 ENDDO
2532 IF (iget(095)>0) then
2533 if(grib=='grib2') then
2534 cfld=cfld+1
2535 fld_info(cfld)%ifld=iavblfld(iget(095))
2536 fld_info(cfld)%lvl=lvlsxml(1,iget(095))
2537!$omp parallel do private(i,j,ii,jj)
2538 do j=1,jend-jsta+1
2539 jj = jsta+j-1
2540 do i=1,iend-ista+1
2541 ii = ista+i-1
2542 datapd(i,j,cfld) = grid1(ii,jj)
2543 enddo
2544 enddo
2545 endif
2546 ENDIF
2547 IF (iget(096)>0) then
2548 if(grib=='grib2') then
2549 cfld=cfld+1
2550 fld_info(cfld)%ifld=iavblfld(iget(096))
2551 fld_info(cfld)%lvl=lvlsxml(1,iget(096))
2552!$omp parallel do private(i,j,ii,jj)
2553 do j=1,jend-jsta+1
2554 jj = jsta+j-1
2555 do i=1,iend-ista+1
2556 ii = ista+i-1
2557 datapd(i,j,cfld) = grid2(ii,jj)
2558 enddo
2559 enddo
2560 endif
2561 ENDIF
2562 ENDIF
2563 ENDIF
2564 ENDIF
2565!
2566! ENDIF FOR BOUNDARY LAYER BLOCK.
2567!
2568 ENDIF
2569!
2570!
2571! ***BLOCK 6: MISCELLANEOUS LAYER MEAN LFM AND NGM FIELDS.
2572!
2573 IF ( (iget(066)>0).OR.(iget(081)>0).OR. &
2574 (iget(082)>0).OR.(iget(104)>0).OR. &
2575 (iget(099)>0).OR.(iget(100)>0).OR. &
2576 (iget(101)>0).OR.(iget(102)>0).OR. &
2577 (iget(103)>0) ) THEN
2578!
2579! LFM "MEAN" RELATIVE HUMIDITIES AND PRECIPITABLE WATER.
2580!
2581 IF ( (iget(066)>0).OR.(iget(081)>0).OR. &
2582 (iget(082)>0).OR.(iget(104)>0) ) THEN
2583 allocate(rh3310(ista:iend,jsta:jend),rh6610(ista:iend,jsta:jend), &
2584 rh3366(ista:iend,jsta:jend),pw3310(ista:iend,jsta:jend))
2585 CALL lfmfld(rh3310,rh6610,rh3366,pw3310)
2586!
2587! SIGMA 0.33-1.00 MEAN RELATIVE HUMIIDITY.
2588 IF (iget(066)>0) THEN
2589!$omp parallel do private(i,j)
2590 DO j=jsta,jend
2591 DO i=ista,iend
2592 grid1(i,j) = rh3310(i,j)
2593 ENDDO
2594 ENDDO
2595 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2596 CALL bound(grid1,h1,h100)
2597 if(grib=='grib2') then
2598 cfld=cfld+1
2599 fld_info(cfld)%ifld=iavblfld(iget(066))
2600 fld_info(cfld)%lvl=lvlsxml(1,iget(066))
2601!$omp parallel do private(i,j,ii,jj)
2602 do j=1,jend-jsta+1
2603 jj = jsta+j-1
2604 do i=1,iend-ista+1
2605 ii = ista+i-1
2606 datapd(i,j,cfld) = grid1(ii,jj)
2607 enddo
2608 enddo
2609! print *,'in miscln,RH0.33-1.0,cfld=',cfld,'fld=', &
2610! IAVBLFLD(IGET(066)),'lvl=',fld_info(cfld)%lvl
2611 endif
2612 ENDIF
2613!
2614! SIGMA 0.66-1.00 MEAN RELATIVE HUMIIDITY.
2615 IF (iget(081)>0) THEN
2616!$omp parallel do private(i,j)
2617 DO j=jsta,jend
2618 DO i=ista,iend
2619 grid1(i,j) = rh6610(i,j)
2620 ENDDO
2621 ENDDO
2622 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2623 CALL bound(grid1,h1,h100)
2624 if(grib=='grib2') then
2625 cfld=cfld+1
2626 fld_info(cfld)%ifld=iavblfld(iget(081))
2627 fld_info(cfld)%lvl=lvlsxml(1,iget(081))
2628!$omp parallel do private(i,j,ii,jj)
2629 do j=1,jend-jsta+1
2630 jj = jsta+j-1
2631 do i=1,iend-ista+1
2632 ii = ista+i-1
2633 datapd(i,j,cfld) = grid1(ii,jj)
2634 enddo
2635 enddo
2636 endif
2637 ENDIF
2638!
2639! SIGMA 0.33-0.66 MEAN RELATIVE HUMIIDITY.
2640 IF (iget(082)>0) THEN
2641!$omp parallel do private(i,j)
2642 DO j=jsta,jend
2643 DO i=ista,iend
2644 grid1(i,j) = rh3366(i,j)
2645 ENDDO
2646 ENDDO
2647 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2648 CALL bound(grid1,h1,h100)
2649 if(grib=='grib2') then
2650 cfld=cfld+1
2651 fld_info(cfld)%ifld=iavblfld(iget(082))
2652 fld_info(cfld)%lvl=lvlsxml(1,iget(082))
2653!$omp parallel do private(i,j,ii,jj)
2654 do j=1,jend-jsta+1
2655 jj = jsta+j-1
2656 do i=1,iend-ista+1
2657 ii = ista+i-1
2658 datapd(i,j,cfld) = grid1(ii,jj)
2659 enddo
2660 enddo
2661 endif
2662 ENDIF
2663!
2664! SIGMA 0.33-1.00 PRECIPITABLE WATER.
2665 IF (iget(104)>0) THEN
2666!$omp parallel do private(i,j)
2667 DO j=jsta,jend
2668 DO i=ista,iend
2669 grid1(i,j) = pw3310(i,j)
2670 ENDDO
2671 ENDDO
2672 CALL bound(grid1,d00,h99999)
2673 if(grib=='grib2') then
2674 cfld=cfld+1
2675 fld_info(cfld)%ifld=iavblfld(iget(104))
2676 fld_info(cfld)%lvl=lvlsxml(1,iget(104))
2677!$omp parallel do private(i,j,ii,jj)
2678 do j=1,jend-jsta+1
2679 jj = jsta+j-1
2680 do i=1,iend-ista+1
2681 ii = ista+i-1
2682 datapd(i,j,cfld) = grid1(ii,jj)
2683 enddo
2684 enddo
2685 endif
2686 ENDIF
2687 deallocate(rh3310,rh6610,rh3366,pw3310)
2688 ENDIF
2689!
2690! VARIOUS LAYER MEAN NGM SIGMA FIELDS.
2691!
2692 IF ( (iget(099)>0).OR.(iget(100)>0).OR. &
2693 (iget(101)>0).OR.(iget(102)>0).OR. &
2694 (iget(103)>0) ) THEN
2695 allocate(rh4710(ista_2l:iend_2u,jsta_2l:jend_2u),rh4796(ista_2l:iend_2u,jsta_2l:jend_2u), &
2696 rh1847(ista_2l:iend_2u,jsta_2l:jend_2u))
2697 allocate(rh8498(ista_2l:iend_2u,jsta_2l:jend_2u),qm8510(ista_2l:iend_2u,jsta_2l:jend_2u))
2698
2699 CALL ngmfld(rh4710,rh4796,rh1847,rh8498,qm8510)
2700!
2701! SIGMA 0.47191-1.00000 RELATIVE HUMIDITY.
2702 IF (iget(099)>0) THEN
2703!$omp parallel do private(i,j)
2704 DO j=jsta,jend
2705 DO i=ista,iend
2706 grid1(i,j) = rh4710(i,j)
2707 ENDDO
2708 ENDDO
2709 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2710 CALL bound(grid1,h1,h100)
2711 if(grib=='grib2') then
2712 cfld=cfld+1
2713 fld_info(cfld)%ifld=iavblfld(iget(099))
2714 fld_info(cfld)%lvl=lvlsxml(1,iget(099))
2715!$omp parallel do private(i,j,ii,jj)
2716 do j=1,jend-jsta+1
2717 jj = jsta+j-1
2718 do i=1,iend-ista+1
2719 ii = ista+i-1
2720 datapd(i,j,cfld) = grid1(ii,jj)
2721 enddo
2722 enddo
2723 endif
2724 ENDIF
2725!
2726! SIGMA 0.47191-0.96470 RELATIVE HUMIDITY.
2727 IF (iget(100)>0) THEN
2728!$omp parallel do private(i,j)
2729 DO j=jsta,jend
2730 DO i=ista,iend
2731 grid1(i,j) = rh4796(i,j)
2732 ENDDO
2733 ENDDO
2734 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2735 CALL bound(grid1,h1,h100)
2736 if(grib=='grib2') then
2737 cfld=cfld+1
2738 fld_info(cfld)%ifld=iavblfld(iget(100))
2739 fld_info(cfld)%lvl=lvlsxml(1,iget(100))
2740!$omp parallel do private(i,j,ii,jj)
2741 do j=1,jend-jsta+1
2742 jj = jsta+j-1
2743 do i=1,iend-ista+1
2744 ii = ista+i-1
2745 datapd(i,j,cfld) = grid1(ii,jj)
2746 enddo
2747 enddo
2748 endif
2749 ENDIF
2750!
2751! SIGMA 0.18019-0.47191 RELATIVE HUMIDITY.
2752 IF (iget(101)>0) THEN
2753!$omp parallel do private(i,j)
2754 DO j=jsta,jend
2755 DO i=ista,iend
2756 grid1(i,j) = rh1847(i,j)
2757 ENDDO
2758 ENDDO
2759 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2760 CALL bound(grid1,h1,h100)
2761 if(grib=='grib2') then
2762 cfld=cfld+1
2763 fld_info(cfld)%ifld=iavblfld(iget(101))
2764 fld_info(cfld)%lvl=lvlsxml(1,iget(101))
2765!$omp parallel do private(i,j,ii,jj)
2766 do j=1,jend-jsta+1
2767 jj = jsta+j-1
2768 do i=1,iend-ista+1
2769 ii = ista+i-1
2770 datapd(i,j,cfld) = grid1(ii,jj)
2771 enddo
2772 enddo
2773 endif
2774 ENDIF
2775!
2776! SIGMA 0.84368-0.98230 RELATIVE HUMIDITY.
2777 IF (iget(102)>0) THEN
2778!$omp parallel do private(i,j)
2779 DO j=jsta,jend
2780 DO i=ista,iend
2781 grid1(i,j) = rh8498(i,j)
2782 ENDDO
2783 ENDDO
2784 CALL sclfld(grid1(ista:iend,jsta:jend),h100,im,jm)
2785 CALL bound(grid1,h1,h100)
2786 if(grib=='grib2') then
2787 cfld=cfld+1
2788 fld_info(cfld)%ifld=iavblfld(iget(102))
2789 fld_info(cfld)%lvl=lvlsxml(1,iget(102))
2790!$omp parallel do private(i,j,ii,jj)
2791 do j=1,jend-jsta+1
2792 jj = jsta+j-1
2793 do i=1,iend-ista+1
2794 ii = ista+i-1
2795 datapd(i,j,cfld) = grid1(ii,jj)
2796 enddo
2797 enddo
2798 endif
2799 ENDIF
2800!
2801! SIGMA 0.85000-1.00000 MOISTURE CONVERGENCE.
2802 IF (iget(103)>0) THEN
2803 grid1=spval
2804! CONVERT TO DIVERGENCE FOR GRIB
2805!$omp parallel do private(i,j)
2806 DO j=jsta,jend
2807 DO i=ista,iend
2808 IF(qm8510(i,j) < spval) grid1(i,j) = -1.0*qm8510(i,j)
2809 ENDDO
2810 ENDDO
2811 if(grib=='grib2') then
2812 cfld=cfld+1
2813 fld_info(cfld)%ifld=iavblfld(iget(103))
2814 fld_info(cfld)%lvl=lvlsxml(1,iget(103))
2815!$omp parallel do private(i,j,ii,jj)
2816 do j=1,jend-jsta+1
2817 jj = jsta+j-1
2818 do i=1,iend-ista+1
2819 ii = ista+i-1
2820 datapd(i,j,cfld) = grid1(ii,jj)
2821 enddo
2822 enddo
2823 endif
2824 ENDIF
2825 deallocate(rh4710,rh4796,rh1847)
2826 deallocate(rh8498,qm8510)
2827 ENDIF
2828 ENDIF
2829
2830 IF ( (iget(318)>0).OR.(iget(319)>0).OR. &
2831 (iget(320)>0))THEN
2832 allocate(rh4410(ista:iend,jsta:jend),rh7294(ista:iend,jsta:jend), &
2833 rh4472(ista:iend,jsta:jend),rh3310(ista:iend,jsta:jend))
2834 CALL lfmfld_gfs(rh4410,rh7294,rh4472,rh3310)
2835!
2836! SIGMA 0.44-1.00 MEAN RELATIVE HUMIIDITY.
2837 IF (iget(318)>0) THEN
2838 grid1=spval
2839!$omp parallel do private(i,j)
2840 DO j=jsta,jend
2841 DO i=ista,iend
2842 IF(rh4410(i,j) < spval) grid1(i,j) = rh4410(i,j)*100.
2843 ENDDO
2844 ENDDO
2845 CALL bound(grid1,d00,h100)
2846 if(grib=='grib2') then
2847 cfld=cfld+1
2848 fld_info(cfld)%ifld=iavblfld(iget(318))
2849 fld_info(cfld)%lvl=lvlsxml(1,iget(318))
2850!$omp parallel do private(i,j,ii,jj)
2851 do j=1,jend-jsta+1
2852 jj = jsta+j-1
2853 do i=1,iend-ista+1
2854 ii = ista+i-1
2855 datapd(i,j,cfld) = grid1(ii,jj)
2856 enddo
2857 enddo
2858 endif
2859 ENDIF
2860!
2861! SIGMA 0.72-0.94 MEAN RELATIVE HUMIIDITY.
2862 IF (iget(319)>0) THEN
2863 grid1=spval
2864!$omp parallel do private(i,j)
2865 DO j=jsta,jend
2866 DO i=ista,iend
2867 IF(rh7294(i,j) < spval) grid1(i,j) = rh7294(i,j)*100.
2868 ENDDO
2869 ENDDO
2870 CALL bound(grid1,d00,h100)
2871 if(grib=='grib2') then
2872 cfld=cfld+1
2873 fld_info(cfld)%ifld=iavblfld(iget(319))
2874 fld_info(cfld)%lvl=lvlsxml(1,iget(319))
2875!$omp parallel do private(i,j,ii,jj)
2876 do j=1,jend-jsta+1
2877 jj = jsta+j-1
2878 do i=1,iend-ista+1
2879 ii = ista+i-1
2880 datapd(i,j,cfld) = grid1(ii,jj)
2881 enddo
2882 enddo
2883 endif
2884 ENDIF
2885!
2886! SIGMA 0.44-0.72 MEAN RELATIVE HUMIIDITY.
2887 IF (iget(320)>0) THEN
2888 grid1=spval
2889!$omp parallel do private(i,j)
2890 DO j=jsta,jend
2891 DO i=ista,iend
2892 IF(rh4472(i,j) < spval) grid1(i,j)=rh4472(i,j)*100.
2893 ENDDO
2894 ENDDO
2895 CALL bound(grid1,d00,h100)
2896 if(grib=='grib2') then
2897 cfld=cfld+1
2898 fld_info(cfld)%ifld=iavblfld(iget(320))
2899 fld_info(cfld)%lvl=lvlsxml(1,iget(320))
2900!$omp parallel do private(i,j,ii,jj)
2901 do j=1,jend-jsta+1
2902 jj = jsta+j-1
2903 do i=1,iend-ista+1
2904 ii = ista+i-1
2905 datapd(i,j,cfld) = grid1(ii,jj)
2906 enddo
2907 enddo
2908 endif
2909 ENDIF
2910 deallocate(rh4410,rh7294,rh4472,rh3310)
2911 END IF
2912
2913! GFS computes sigma=0.9950 T, THETA, U, V from lowest two model level fields
2914 IF ( (iget(321)>0).OR.(iget(322)>0).OR. &
2915 (iget(323)>0).OR.(iget(324)>0).OR. &
2916 (iget(325)>0).OR.(iget(326)>0)) THEN
2917!$omp parallel do private(i,j)
2918 DO j=jsta,jend
2919 DO i=ista,iend
2920 egrid2(i,j) = 0.995*pint(i,j,lm+1)
2921 egrid1(i,j) = log(pmid(i,j,lm)/egrid2(i,j)) &
2922 / log(pmid(i,j,lm)/pmid(i,j,lm-1))
2923
2924 IF (modelname == 'GFS' .OR. modelname == 'FV3R') THEN
2925 egrid1(i,j) = log(pmid(i,j,lm)/egrid2(i,j)) &
2926 / max(1.e-6,log(pmid(i,j,lm)/pmid(i,j,lm-1)))
2927 egrid1(i,j) =max(-10.0,min(egrid1(i,j), 10.0))
2928 IF ( abs(pmid(i,j,lm)-pmid(i,j,lm-1)) < 0.5 ) THEN
2929 egrid1(i,j) = -1.
2930 ENDIF
2931 ENDIF
2932
2933 END DO
2934 END DO
2935! Temperature
2936 IF (iget(321)>0) THEN
2937 grid1=spval
2938!$omp parallel do private(i,j)
2939 DO j=jsta,jend
2940 DO i=ista,iend
2941 IF(t(i,j,lm)<spval.and.t(i,j,lm-1)<spval.and.egrid1(i,j)<spval)&
2942 grid1(i,j) = t(i,j,lm)+(t(i,j,lm-1)-t(i,j,lm)) &
2943 * egrid1(i,j)
2944 ENDDO
2945 ENDDO
2946 if(grib=='grib2') then
2947 cfld=cfld+1
2948 fld_info(cfld)%ifld=iavblfld(iget(321))
2949 fld_info(cfld)%lvl=lvlsxml(1,iget(321))
2950!$omp parallel do private(i,j,ii,jj)
2951 do j=1,jend-jsta+1
2952 jj = jsta+j-1
2953 do i=1,iend-ista+1
2954 ii = ista+i-1
2955 datapd(i,j,cfld) = grid1(ii,jj)
2956 enddo
2957 enddo
2958 endif
2959! print *,'in miscln,iget(321,temp,sigmadata=',maxval(GRID1(1:im,jsta:jend)), &
2960! minval(GRID1(1:im,jsta:jend)),'grib=',grib
2961 ENDIF
2962! Potential Temperature
2963 IF (iget(322)>0) THEN
2964 grid2=spval
2965!$omp parallel do private(i,j)
2966 DO j=jsta,jend
2967 DO i=ista,iend
2968 IF(t(i,j,lm)<spval.and.t(i,j,lm-1)<spval.and.egrid1(i,j)<spval)&
2969 grid2(i,j) = t(i,j,lm)+(t(i,j,lm-1)-t(i,j,lm)) &
2970 * egrid1(i,j)
2971 ENDDO
2972 ENDDO
2973 CALL calpot(egrid2,grid2(ista:iend,jsta:jend),grid1(ista:iend,jsta:jend))
2974 if(grib=='grib2') then
2975 cfld=cfld+1
2976 fld_info(cfld)%ifld=iavblfld(iget(322))
2977 fld_info(cfld)%lvl=lvlsxml(1,iget(322))
2978!$omp parallel do private(i,j,ii,jj)
2979 do j=1,jend-jsta+1
2980 jj = jsta+j-1
2981 do i=1,iend-ista+1
2982 ii = ista+i-1
2983 datapd(i,j,cfld) = grid1(ii,jj)
2984 enddo
2985 enddo
2986 endif
2987 ENDIF
2988! RH
2989 IF (iget(323)>0) THEN
2990 grid1=spval
2991!$omp parallel do private(i,j,es1,qs1,rh1,es2,qs2,rh2)
2992 DO j=jsta,jend
2993 DO i=ista,iend
2994 IF(pmid(i,j,lm)<spval.and.pmid(i,j,lm-1)<spval.and.&
2995 q(i,j,lm)<spval.and.q(i,j,lm-1)<spval) THEN
2996 es1 = min(pmid(i,j,lm),fpvsnew(t(i,j,lm)))
2997 qs1 = con_eps*es1/(pmid(i,j,lm)+con_epsm1*es1)
2998 rh1 = q(i,j,lm)/qs1
2999 es2 = min(pmid(i,j,lm-1),fpvsnew(t(i,j,lm-1)))
3000 qs2 = con_eps*es2/(pmid(i,j,lm-1)+con_epsm1*es2)
3001 rh2 = q(i,j,lm-1)/qs2
3002 grid1(i,j) = (rh1+(rh2-rh1)*egrid1(i,j))*100.
3003 ENDIF
3004 ENDDO
3005 ENDDO
3006 CALL bound(grid1,d00,h100)
3007 if(grib=='grib2') then
3008 cfld=cfld+1
3009 fld_info(cfld)%ifld=iavblfld(iget(323))
3010 fld_info(cfld)%lvl=lvlsxml(1,iget(323))
3011!$omp parallel do private(i,j,ii,jj)
3012 do j=1,jend-jsta+1
3013 jj = jsta+j-1
3014 do i=1,iend-ista+1
3015 ii = ista+i-1
3016 datapd(i,j,cfld) = grid1(ii,jj)
3017 enddo
3018 enddo
3019 endif
3020 ENDIF
3021! U
3022 IF (iget(324)>0) THEN
3023 grid1=spval
3024!$omp parallel do private(i,j)
3025 DO j=jsta,jend
3026 DO i=ista,iend
3027 IF(uh(i,j,lm)<spval.and.uh(i,j,lm-1)<spval.and.egrid1(i,j)<spval)&
3028 grid1(i,j) = uh(i,j,lm)+(uh(i,j,lm-1)-uh(i,j,lm)) &
3029 * egrid1(i,j)
3030 ENDDO
3031 ENDDO
3032 if(grib=='grib2') then
3033 cfld=cfld+1
3034 fld_info(cfld)%ifld=iavblfld(iget(324))
3035 fld_info(cfld)%lvl=lvlsxml(1,iget(324))
3036!$omp parallel do private(i,j,ii,jj)
3037 do j=1,jend-jsta+1
3038 jj = jsta+j-1
3039 do i=1,iend-ista+1
3040 ii = ista+i-1
3041 datapd(i,j,cfld) = grid1(ii,jj)
3042 enddo
3043 enddo
3044 endif
3045 ENDIF
3046! V
3047 IF (iget(325)>0) THEN
3048 grid1=spval
3049!$omp parallel do private(i,j)
3050 DO j=jsta,jend
3051 DO i=ista,iend
3052 IF(vh(i,j,lm)<spval.and.vh(i,j,lm-1)<spval.and.egrid1(i,j)<spval)&
3053 grid1(i,j) = vh(i,j,lm)+(vh(i,j,lm-1)-vh(i,j,lm)) &
3054 * egrid1(i,j)
3055 ENDDO
3056 ENDDO
3057 if(grib=='grib2') then
3058 cfld=cfld+1
3059 fld_info(cfld)%ifld=iavblfld(iget(325))
3060 fld_info(cfld)%lvl=lvlsxml(1,iget(325))
3061!$omp parallel do private(i,j,ii,jj)
3062 do j=1,jend-jsta+1
3063 jj = jsta+j-1
3064 do i=1,iend-ista+1
3065 ii = ista+i-1
3066 datapd(i,j,cfld) = grid1(ii,jj)
3067 enddo
3068 enddo
3069 endif
3070 ENDIF
3071! OMGA
3072 IF (iget(326)>0) THEN
3073 grid1=spval
3074!$omp parallel do private(i,j)
3075 DO j=jsta,jend
3076 DO i=ista,iend
3077 IF(omga(i,j,lm)<spval.and.omga(i,j,lm-1)<spval.and.egrid1(i,j)<spval)&
3078 grid1(i,j) = omga(i,j,lm)+(omga(i,j,lm-1)-omga(i,j,lm))&
3079 * egrid1(i,j)
3080 ENDDO
3081 ENDDO
3082 if(grib=='grib2') then
3083 cfld=cfld+1
3084 fld_info(cfld)%ifld=iavblfld(iget(326))
3085 fld_info(cfld)%lvl=lvlsxml(1,iget(326))
3086!$omp parallel do private(i,j,ii,jj)
3087 do j=1,jend-jsta+1
3088 jj = jsta+j-1
3089 do i=1,iend-ista+1
3090 ii = ista+i-1
3091 datapd(i,j,cfld) = grid1(ii,jj)
3092 enddo
3093 enddo
3094 endif
3095 ENDIF
3096 END IF
3097
3098! MIXED LAYER CAPE AND CINS
3099!
3100 field1=.false.
3101 field2=.false.
3102!
3103 IF(iget(032)>0)THEN
3104 IF(lvls(3,iget(032))>0)field1=.true.
3105 ENDIF
3106 IF(iget(107)>0)THEN
3107 IF(lvls(3,iget(107))>0)field2=.true.
3108 ENDIF
3109!
3110 IF(iget(582)>0)THEN
3111 field1=.true.
3112 ENDIF
3113 IF(iget(583)>0)THEN
3114 field2=.true.
3115 ENDIF
3116
3117 IF(field1.OR.field2)THEN
3118 itype = 2
3119!
3120!$omp parallel do private(i,j)
3121 DO j=jsta,jend
3122 DO i=ista,iend
3123 egrid1(i,j) = -h99999
3124 egrid2(i,j) = -h99999
3125 lb2(i,j) = (lvlbnd(i,j,1) + lvlbnd(i,j,2) + &
3126 lvlbnd(i,j,3))/3
3127 p1d(i,j) = (pbnd(i,j,1) + pbnd(i,j,2) + pbnd(i,j,3))/3
3128 t1d(i,j) = (tbnd(i,j,1) + tbnd(i,j,2) + tbnd(i,j,3))/3
3129 q1d(i,j) = (qbnd(i,j,1) + qbnd(i,j,2) + qbnd(i,j,3))/3
3130 ENDDO
3131 ENDDO
3132!
3133 dpbnd = 0.
3134 CALL calcape(itype,dpbnd,p1d,t1d,q1d,lb2,egrid1, &
3135 egrid2,egrid3,egrid4,egrid5)
3136
3137 IF (iget(582)>0) THEN
3138! dong add missing value for cape
3139 grid1=spval
3140!$omp parallel do private(i,j)
3141 DO j=jsta,jend
3142 DO i=ista,iend
3143 IF(t1d(i,j) < spval) THEN
3144 grid1(i,j) = egrid1(i,j)
3145 IF (submodelname == 'RTMA')mlcape(i,j)=grid1(i,j)
3146 ENDIF
3147 ENDDO
3148 ENDDO
3149 CALL bound(grid1,d00,h99999)
3150 if(grib=='grib2') then
3151 cfld=cfld+1
3152 fld_info(cfld)%ifld=iavblfld(iget(582))
3153 fld_info(cfld)%lvl=lvlsxml(1,iget(582))
3154!$omp parallel do private(i,j,ii,jj)
3155 do j=1,jend-jsta+1
3156 jj = jsta+j-1
3157 do i=1,iend-ista+1
3158 ii = ista+i-1
3159 datapd(i,j,cfld) = grid1(ii,jj)
3160 enddo
3161 enddo
3162 endif
3163 ENDIF
3164 IF (iget(583)>0) THEN
3165! dong add missing value for cape
3166 grid1=spval
3167!$omp parallel do private(i,j)
3168 DO j=jsta,jend
3169 DO i=ista,iend
3170 IF(t1d(i,j) < spval) grid1(i,j) = - egrid2(i,j)
3171 ENDDO
3172 ENDDO
3173!
3174 CALL bound(grid1,d00,h99999)
3175!
3176!$omp parallel do private(i,j)
3177 DO j=jsta,jend
3178 DO i=ista,iend
3179 IF(t1d(i,j) < spval) THEN
3180 grid1(i,j) = - grid1(i,j)
3181 IF (submodelname == 'RTMA') mlcin(i,j)=grid1(i,j)
3182 ENDIF
3183 ENDDO
3184 ENDDO
3185!
3186 if(grib=='grib2') then
3187 cfld=cfld+1
3188 fld_info(cfld)%ifld=iavblfld(iget(583))
3189 fld_info(cfld)%lvl=lvlsxml(1,iget(583))
3190!$omp parallel do private(i,j,ii,jj)
3191 do j=1,jend-jsta+1
3192 jj = jsta+j-1
3193 do i=1,iend-ista+1
3194 ii = ista+i-1
3195 datapd(i,j,cfld) = grid1(ii,jj)
3196 enddo
3197 enddo
3198 endif
3199 ENDIF
3200 ENDIF
3201
3202! MIXED LAYER LIFTING CONDENSATION PRESSURE AND HEIGHT.
3203! EGRID1 IS LCL PRESSURE. EGRID2 IS LCL HEIGHT.
3204!
3205 IF ( (iget(109)>0).OR.(iget(110)>0) ) THEN
3206 CALL callcl(p1d,t1d,q1d,egrid1,egrid2)
3207 IF (iget(109)>0) THEN
3208 grid1=spval
3209 DO j=jsta,jend
3210 DO i=ista,iend
3211 IF(t1d(i,j) < spval) grid1(i,j)=egrid2(i,j)
3212 IF (submodelname == 'RTMA') mllcl(i,j) = grid1(i,j)
3213 ENDDO
3214 ENDDO
3215!
3216! ID(1:25) = 0
3217!
3218! CALL GRIBIT(IGET(109),1,
3219! X GRID1,IM,JM)
3220 ENDIF
3221!
3222! IF (IGET(110)>0) THEN
3223! DO J=JSTA,JEND
3224! DO I=ISTA,IEND
3225! GRID1(I,J)=EGRID1(I,J)
3226! ENDDO
3227! ENDDO
3228!
3229! ID(1:25) = 0
3230!
3231! CALL GRIBIT(IGET(110),1,
3232! X GRID1,IM,JM)
3233! ENDIF
3234 ENDIF
3235!
3236! MOST UNSTABLE CAPE-LOWEST 300 MB
3237!
3238 field1=.false.
3239 field2=.false.
3240!
3241 IF(iget(032)>0)THEN
3242 IF(lvls(4,iget(032))>0)field1=.true.
3243 ENDIF
3244
3245 IF(iget(107)>0)THEN
3246 IF(lvls(4,iget(107))>0)field2=.true.
3247 ENDIF
3248!
3249 IF(iget(584)>0)THEN
3250 field1=.true.
3251 ENDIF
3252 IF(iget(585)>0)THEN
3253 field2=.true.
3254 ENDIF
3255!
3256 IF(field1.OR.field2.OR.need_ifi)THEN
3257 itype = 1
3258 call allocate_cape_arrays
3259!
3260!$omp parallel do private(i,j)
3261 DO j=jsta,jend
3262 DO i=ista,iend
3263 egrid1(i,j) = -h99999
3264 egrid2(i,j) = -h99999
3265 ENDDO
3266 ENDDO
3267
3268 dpbnd = 300.e2
3269 CALL calcape(itype,dpbnd,p1d,t1d,q1d,lb2,egrid1, &
3270 egrid2,egrid3,egrid4,egrid5)
3271 IF (iget(584)>0 .or. need_ifi) THEN
3272! dong add missing value to cin
3273 grid1 = spval
3274!$omp parallel do private(i,j)
3275 DO j=jsta,jend
3276 DO i=ista,iend
3277 IF(t1d(i,j) < spval) THEN
3278 grid1(i,j) = egrid1(i,j)
3279 IF (submodelname == 'RTMA') mucape(i,j)=grid1(i,j)
3280 ENDIF
3281 ENDDO
3282 ENDDO
3283 CALL bound(grid1,d00,h99999)
3284! IF (SUBMODELNAME == 'RTMA') THEN
3285! CALL BOUND(MUCAPE,D00,H99999)
3286! ENDIF
3287!$omp parallel do private(i,j)
3288 DO j=jsta,jend
3289 DO i=ista,iend
3290 cape(i,j) = grid1(i,j)
3291 ENDDO
3292 ENDDO
3293 if(iget(584)>0 .and. grib=='grib2') then
3294 cfld=cfld+1
3295 fld_info(cfld)%ifld=iavblfld(iget(584))
3296 fld_info(cfld)%lvl=lvlsxml(1,iget(584))
3297!$omp parallel do private(i,j,ii,jj)
3298 do j=1,jend-jsta+1
3299 jj = jsta+j-1
3300 do i=1,iend-ista+1
3301 ii = ista+i-1
3302 datapd(i,j,cfld) = grid1(ii,jj)
3303 enddo
3304 enddo
3305 endif
3306
3307 ENDIF
3308
3309 IF (iget(585)>0 .or. need_ifi) THEN
3310! dong add missing value to cin
3311 grid1 = spval
3312!$omp parallel do private(i,j)
3313 DO j=jsta,jend
3314 DO i=ista,iend
3315 IF(t1d(i,j) < spval) grid1(i,j) = - egrid2(i,j)
3316 ENDDO
3317 ENDDO
3318 CALL bound(grid1,d00,h99999)
3319 DO j=jsta,jend
3320 DO i=ista,iend
3321 IF(t1d(i,j) < spval) THEN
3322 grid1(i,j) = - grid1(i,j)
3323 IF (submodelname == 'RTMA')THEN
3324 mucape(i,j) = grid1(i,j)
3325 muq1d(i,j) = q1d(i,j)
3326 ENDIF
3327 ENDIF
3328 ENDDO
3329 ENDDO
3330
3331!$omp parallel do private(i,j)
3332 DO j=jsta,jend
3333 DO i=ista,iend
3334 cin(i,j) = grid1(i,j)
3335 ENDDO
3336 ENDDO
3337
3338 if(iget(585)>0 .and. grib=='grib2') then
3339 cfld=cfld+1
3340 fld_info(cfld)%ifld=iavblfld(iget(585))
3341 fld_info(cfld)%lvl=lvlsxml(1,iget(585))
3342!$omp parallel do private(i,j,ii,jj)
3343 do j=1,jend-jsta+1
3344 jj = jsta+j-1
3345 do i=1,iend-ista+1
3346 ii = ista+i-1
3347 datapd(i,j,cfld) = grid1(ii,jj)
3348 enddo
3349 enddo
3350 endif
3351
3352 ENDIF
3353
3354! EQUILLIBRIUM HEIGHT
3355 IF (iget(443)>0) THEN
3356 grid1 = spval
3357!$omp parallel do private(i,j)
3358 DO j=jsta,jend
3359 DO i=ista,iend
3360 IF(t1d(i,j) < spval) grid1(i,j) = egrid4(i,j)
3361 ENDDO
3362 ENDDO
3363 if(grib=='grib2') then
3364 cfld=cfld+1
3365 fld_info(cfld)%ifld=iavblfld(iget(443))
3366 fld_info(cfld)%lvl=lvlsxml(1,iget(443))
3367!$omp parallel do private(i,j,ii,jj)
3368 do j=1,jend-jsta+1
3369 jj = jsta+j-1
3370 do i=1,iend-ista+1
3371 ii = ista+i-1
3372 datapd(i,j,cfld) = grid1(ii,jj)
3373 enddo
3374 enddo
3375 endif
3376 ENDIF
3377!Equilibrium Temperature
3378 IF (iget(982)>0) THEN
3379 DO j=jsta,jend
3380 DO i=ista,iend
3381 grid1(i,j) = teql(i,j)
3382 ENDDO
3383 ENDDO
3384 if(grib=='grib2') then
3385 cfld=cfld+1
3386 fld_info(cfld)%ifld=iavblfld(iget(982))
3387 fld_info(cfld)%lvl=lvlsxml(1,iget(982))
3388!$omp parallel do private(i,j,ii,jj)
3389 do j=1,jend-jsta+1
3390 jj = jsta+j-1
3391 do i=1,iend-ista+1
3392 ii = ista+i-1
3393 datapd(i,j,cfld) = grid1(ii,jj)
3394 enddo
3395 enddo
3396 endif
3397 ENDIF
3398
3399
3400! PRESSURE OF LEVEL FROM WHICH 300 MB MOST UNSTABLE CAPE
3401! PARCEL WAS LIFTED (eq. PRESSURE OF LEVEL OF HIGHEST THETA-E)
3402 IF (iget(246)>0) THEN
3403 grid1 = spval
3404!$omp parallel do private(i,j)
3405 DO j=jsta,jend
3406 DO i=ista,iend
3407 IF(t1d(i,j) < spval) grid1(i,j) = egrid3(i,j)
3408 ENDDO
3409 ENDDO
3410 CALL bound(grid1,d00,h99999)
3411! print *,'in miscln,PLPL=',maxval(grid1(1:im,jsta:jend)), &
3412! minval(grid1(1:im,jsta:jend))
3413 if(grib=='grib2') then
3414 cfld=cfld+1
3415 fld_info(cfld)%ifld=iavblfld(iget(246))
3416 fld_info(cfld)%lvl=lvlsxml(1,iget(246))
3417!$omp parallel do private(i,j,ii,jj)
3418 do j=1,jend-jsta+1
3419 jj = jsta+j-1
3420 do i=1,iend-ista+1
3421 ii = ista+i-1
3422 datapd(i,j,cfld) = grid1(ii,jj)
3423 enddo
3424 enddo
3425 endif
3426 ENDIF ! 246
3427
3428! GENERAL THUNDER PARAMETER ??? 458 ???
3429 IF (iget(444)>0) THEN
3430 grid1 = spval
3431!$omp parallel do private(i,j)
3432 DO j=jsta,jend
3433 DO i=ista,iend
3434 IF(cprate(i,j) < spval) THEN
3435 IF (cprate(i,j) > pthresh) THEN
3436 grid1(i,j) = egrid5(i,j)
3437 ELSE
3438 grid1(i,j) = 0
3439 ENDIF
3440 ENDIF
3441 ENDDO
3442 ENDDO
3443 CALL bound(grid1,d00,h99999)
3444 if(grib=='grib2') then
3445 cfld=cfld+1
3446 fld_info(cfld)%ifld=iavblfld(iget(444))
3447 fld_info(cfld)%lvl=lvlsxml(1,iget(444))
3448!$omp parallel do private(i,j,ii,jj)
3449 do j=1,jend-jsta+1
3450 jj = jsta+j-1
3451 do i=1,iend-ista+1
3452 ii = ista+i-1
3453 datapd(i,j,cfld) = grid1(ii,jj)
3454 enddo
3455 enddo
3456 endif
3457 ENDIF
3458 ENDIF
3459
3460
3461 IF (submodelname == 'RTMA')THEN
3462
3463!
3464! --- Effective (inflow) Layer (EL)
3465!
3466
3467 ALLOCATE(el_base(ista_2l:iend_2u,jsta_2l:jend_2u))
3468 ALLOCATE(el_tops(ista_2l:iend_2u,jsta_2l:jend_2u))
3469 ALLOCATE(found_base(ista_2l:iend_2u,jsta_2l:jend_2u))
3470 ALLOCATE(found_tops(ista_2l:iend_2u,jsta_2l:jend_2u))
3471!$omp parallel do private(i,j)
3472 DO j=jsta,jend
3473 DO i=ista,iend
3474 el_base(i,j) = lm
3475 el_tops(i,j) = lm
3476 found_base(i,j) = .false.
3477 found_tops(i,j) = .false.
3478 ENDDO
3479 ENDDO
3480
3481!
3482
3483 itype = 2
3484 dpbnd = 0.
3485
3486 DO l = lm, 1, -1
3487
3488! SET AIR PARCELS FOR LEVEL L
3489!$omp parallel do private(i,j)
3490 DO j=jsta,jend
3491 DO i=ista,iend
3492 egrid1(i,j) = -h99999
3493 egrid2(i,j) = -h99999
3494 idummy(i,j) = 0
3495 p1d(i,j) = pmid(i,j,l)
3496 t1d(i,j) = t(i,j,l)
3497 q1d(i,j) = q(i,j,l)
3498 ENDDO
3499 ENDDO
3500
3501!--- CALCULATE CAPE/CIN FOR ALL AIR PARCELS on LEVEL L
3502 IF (debugprint) WRITE(1000+me,'(A,I3)') &
3503 ' CALCULATING CAPE/CINS ON LEVEL:',l
3504 CALL calcape(itype,dpbnd,p1d,t1d,q1d,idummy,egrid1, &
3505 egrid2,egrid3,egrid4,egrid5)
3506
3507!--- CHECK CAPE/CIN OF EACH AIR PARCELS WITH EL CRITERIA
3508!$omp parallel do private(i,j)
3509 DO j=jsta,jend
3510 DO i=ista,iend
3511 IF ( .NOT. found_base(i,j) ) THEN
3512 IF ( egrid1(i,j) >= 100. .AND. egrid2(i,j) >= -250. ) THEN
3513 el_base(i,j) = l
3514 found_base(i,j) = .true.
3515 ELSE
3516 el_base(i,j) = lm
3517 found_base(i,j) = .false.
3518 END IF
3519 ELSE
3520 IF ( .NOT. found_tops(i,j) ) THEN
3521 IF ( egrid1(i,j) < 100. .OR. egrid2(i,j) < -250. ) THEN
3522 el_tops(i,j) = l + 1
3523 found_tops(i,j) = .true.
3524 ELSE
3525 el_tops(i,j) = lm
3526 found_tops(i,j) = .false.
3527 END IF
3528 END IF
3529 END IF
3530 ENDDO
3531 ENDDO
3532
3533 END DO ! L
3534
3535
3536 IF (ALLOCATED(found_base)) DEALLOCATE(found_base)
3537 IF (ALLOCATED(found_tops)) DEALLOCATE(found_tops)
3538
3539 IF (debugprint) THEN
3540 WRITE(im_ch,'(I5.5)') im
3541 WRITE(jsta_ch,'(I5.5)') jsta
3542 WRITE(jend_ch,'(I5.5)') jend
3543 effl_fname="EFFL_NEW_"//im_ch//"_"//jsta_ch//"_"//jend_ch &
3544 //".dat"
3545 effl_fname2="EFFL_NEW_LVLS_"//im_ch//"_"//jsta_ch//"_"//jend_ch &
3546 //".dat"
3547 iunit=10000+jsta
3548 iunit2=20000+jsta
3549 irec=0
3550 irec2=0
3551 OPEN(iunit,file=trim(adjustl(effl_fname)),form='FORMATTED')
3552
3553
3554 DO j=jsta,jend
3555 DO i=ista,iend
3556 irec = irec + 1
3557 irec2 = irec2 + 1
3558 WRITE(iunit,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') i, j, &
3559 el_base(i,j),pmid(i,j,el_base(i,j)), &
3560 el_tops(i,j),pmid(i,j,el_tops(i,j))
3561 END DO
3562 ENDDO
3563 CLOSE(iunit)
3564 ENDIF
3565
3566 IF(ALLOCATED(tpar_base)) DEALLOCATE(tpar_base)
3567 IF(ALLOCATED(tpar_tops)) DEALLOCATE(tpar_tops)
3568
3569 ENDIF
3570!
3571! EXPAND HRRR CAPE/CIN RELATED VARIABLES
3572!
3573! CAPE AND CINS 0-3KM, FOLLOW ML PROCEDURE WITH HEIGHT 0-3KM
3574!
3575 field1=.false.
3576 field2=.false.
3577!
3578 IF(iget(032)>0)THEN
3579 IF(lvls(3,iget(032))>0)field1=.true.
3580 ENDIF
3581 IF(iget(107)>0)THEN
3582 IF(lvls(3,iget(107))>0)field2=.true.
3583 ENDIF
3584!
3585 IF(iget(950)>0)THEN
3586 field1=.true.
3587 ENDIF
3588 IF(iget(951)>0)THEN
3589 field2=.true.
3590 ENDIF
3591 IF(modelname == "FV3R" .and. submodelname == "RTMA") THEN
3592 field1=.true.
3593 field2=.true.
3594 ENDIF
3595!
3596! IF(FIELD1)ITYPE=2
3597! IF(FIELD2)ITYPE=2
3598
3599 IF(field1.OR.field2)THEN
3600 itype = 2
3601
3602 call allocate_cape_arrays
3603
3604!
3605!$omp parallel do private(i,j)
3606 DO j=jsta,jend
3607 DO i=ista,iend
3608 egrid1(i,j) = -h99999
3609 egrid2(i,j) = -h99999
3610 egrid3(i,j) = -h99999
3611 egrid4(i,j) = -h99999
3612 egrid5(i,j) = -h99999
3613 egrid6(i,j) = -h99999
3614 egrid7(i,j) = -h99999
3615 egrid8(i,j) = -h99999
3616! ENDDO
3617! ENDDO
3618! DO J=JSTA,JEND
3619! DO I=ISTA,IEND
3620 lb2(i,j) = (lvlbnd(i,j,1) + lvlbnd(i,j,2) + &
3621 lvlbnd(i,j,3))/3
3622 p1d(i,j) = (pbnd(i,j,1) + pbnd(i,j,2) + pbnd(i,j,3))/3
3623 t1d(i,j) = (tbnd(i,j,1) + tbnd(i,j,2) + tbnd(i,j,3))/3
3624 q1d(i,j) = (qbnd(i,j,1) + qbnd(i,j,2) + qbnd(i,j,3))/3
3625 ENDDO
3626 ENDDO
3627
3628 dpbnd = 0.
3629 CALL calcape2(itype,dpbnd,p1d,t1d,q1d,lb2, &
3630 egrid1,egrid2,egrid3,egrid4,egrid5, &
3631 egrid6,egrid7,egrid8)
3632
3633! CAPE1, CINS2, LFC3, ESRHL4,ESRHH5,
3634! DCAPE6,DGLD7, ESP8)
3635!
3636 IF (iget(950)>0) THEN
3637! dong add missing value for cape
3638 grid1=spval
3639!$omp parallel do private(i,j)
3640 DO j=jsta,jend
3641 DO i=ista,iend
3642 IF(t1d(i,j) < spval) grid1(i,j) = egrid1(i,j)
3643 ENDDO
3644 ENDDO
3645 CALL bound(grid1,d00,h99999)
3646 if(grib=='grib2') then
3647 cfld=cfld+1
3648 fld_info(cfld)%ifld=iavblfld(iget(950))
3649 fld_info(cfld)%lvl=lvlsxml(1,iget(950))
3650!$omp parallel do private(i,j,ii,jj)
3651 do j=1,jend-jsta+1
3652 jj = jsta+j-1
3653 do i=1,iend-ista+1
3654 ii = ista+i-1
3655 datapd(i,j,cfld) = grid1(ii,jj)
3656 enddo
3657 enddo
3658 endif
3659 ENDIF !950
3660!
3661 IF (iget(951)>0) THEN
3662! dong add missing value for cape
3663 grid1=spval
3664!$omp parallel do private(i,j)
3665 DO j=jsta,jend
3666 DO i=ista,iend
3667 IF(t1d(i,j) < spval) grid1(i,j) = - egrid2(i,j)
3668 ENDDO
3669 ENDDO
3670!
3671 CALL bound(grid1,d00,h99999)
3672!
3673!$omp parallel do private(i,j)
3674 DO j=jsta,jend
3675 DO i=ista,iend
3676 IF(t1d(i,j) < spval) grid1(i,j) = - grid1(i,j)
3677 ENDDO
3678 ENDDO
3679!
3680 if(grib=='grib2') then
3681 cfld=cfld+1
3682 fld_info(cfld)%ifld=iavblfld(iget(951))
3683 fld_info(cfld)%lvl=lvlsxml(1,iget(951))
3684!$omp parallel do private(i,j,ii,jj)
3685 do j=1,jend-jsta+1
3686 jj = jsta+j-1
3687 do i=1,iend-ista+1
3688 ii = ista+i-1
3689 datapd(i,j,cfld) = grid1(ii,jj)
3690 enddo
3691 enddo
3692 endif
3693 ENDIF !951
3694
3695! LFC HEIGHT
3696
3697 IF (iget(952)>0) THEN
3698 grid1=spval
3699!$omp parallel do private(i,j)
3700 DO j=jsta,jend
3701 DO i=ista,iend
3702 IF(t1d(i,j) < spval) grid1(i,j) = egrid3(i,j)
3703 ENDDO
3704 ENDDO
3705 CALL bound(grid1,d00,h99999)
3706 if(grib=='grib2') then
3707 cfld=cfld+1
3708 fld_info(cfld)%ifld=iavblfld(iget(952))
3709 fld_info(cfld)%lvl=lvlsxml(1,iget(952))
3710!$omp parallel do private(i,j,ii,jj)
3711 do j=1,jend-jsta+1
3712 jj = jsta+j-1
3713 do i=1,iend-ista+1
3714 ii = ista+i-1
3715 datapd(i,j,cfld) = grid1(ii,jj)
3716 enddo
3717 enddo
3718 endif
3719 ENDIF !952
3720
3721
3722! EFFECTIVE STORM RELATIVE HELICITY AND STORM MOTION.
3723
3724 allocate(ust(ista_2l:iend_2u,jsta_2l:jend_2u),vst(ista_2l:iend_2u,jsta_2l:jend_2u), &
3725 heli(ista_2l:iend_2u,jsta_2l:jend_2u,2))
3726 allocate(llow(ista_2l:iend_2u,jsta_2l:jend_2u),lupp(ista_2l:iend_2u,jsta_2l:jend_2u), &
3727 cangle(ista_2l:iend_2u,jsta_2l:jend_2u))
3728 allocate(llow_zint(ista_2l:iend_2u,jsta_2l:jend_2u), &
3729 ieql_zint(ista_2l:iend_2u,jsta_2l:jend_2u),z_temp(ista_2l:iend_2u,jsta_2l:jend_2u))
3730 allocate(midcal(ista_2l:iend_2u,jsta_2l:jend_2u,1:lm))
3731 allocate(z_midcal(ista_2l:iend_2u,jsta_2l:jend_2u))
3732 iget1 = iget(953)
3733 iget2 = -1
3734 iget3 = -1
3735 if (iget1 > 0) then
3736 iget2 = lvls(1,iget1)
3737 iget3 = lvls(2,iget1)
3738 endif
3739 if(me==0) write(*,*) '953 ',iget1,iget2,iget3
3740 IF (iget1 > 0 .OR. iget(162) > 0 .OR. iget(953) > 0) THEN
3741 depth(1) = 3000.0
3742 depth(2) = 1000.0
3743 IF (submodelname == 'RTMA') THEN
3744!--- IF USSING EL BASE & TOP COMPUTED BY NEW SCHEME FOR THE
3745!RELATED VARIABLES
3746!$omp parallel do private(i,j)
3747 DO j=jsta,jend
3748 DO i=ista,iend
3749 llow(i,j) = el_base(i,j)
3750 lupp(i,j) = el_tops(i,j)
3751 ENDDO
3752 ENDDO
3753 ELSE
3754!$omp parallel do private(i,j)
3755 DO j=jsta,jend
3756 DO i=ista,iend
3757 llow(i,j) = int(egrid4(i,j))
3758 lupp(i,j) = int(egrid5(i,j))
3759 ENDDO
3760 ENDDO
3761 ENDIF
3762!--- OUTPUT EL BASE & TOP COMPUTED BY OLD SCHEME
3763 IF (debugprint) THEN
3764 WRITE(im_ch,'(I5.5)') im
3765 WRITE(jsta_ch,'(I5.5)') jsta
3766 WRITE(jend_ch,'(I5.5)') jend
3767 effl_fname="EFFL_OLD_"//im_ch//"_"//jsta_ch//"_"//jend_ch &
3768 //".dat"
3769 iunit=10000+jsta
3770 irec=0
3771 OPEN(iunit,file=trim(adjustl(effl_fname)),form='FORMATTED')
3772 DO j=jsta,jend
3773 DO i=ista,iend
3774 irec = irec + 1
3775 WRITE(iunit,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') i, j, &
3776 llow(i,j),pmid(i,j,llow(i,j)), &
3777 lupp(i,j),pmid(i,j,lupp(i,j))
3778 END DO
3779 ENDDO
3780 CLOSE(iunit)
3781 ENDIF
3782
3783
3784 CALL calhel2(llow,lupp,depth,ust,vst,heli,cangle)
3785!
3786 IF (iget2 > 0) then
3787!$omp parallel do private(i,j)
3788 DO j=jsta,jend
3789 DO i=ista,iend
3790 grid1(i,j) = heli(i,j,1)
3791 ENDDO
3792 ENDDO
3793 if(grib=='grib2') then
3794 cfld=cfld+1
3795 fld_info(cfld)%ifld=iavblfld(iget1)
3796 fld_info(cfld)%lvl=lvlsxml(1,iget1)
3797!$omp parallel do private(i,j,ii,jj)
3798 do j=1,jend-jsta+1
3799 jj = jsta+j-1
3800 do i=1,iend-ista+1
3801 ii = ista+i-1
3802 datapd(i,j,cfld) = grid1(ii,jj)
3803 enddo
3804 enddo
3805 endif
3806 ENDIF
3807
3808 ENDIF !953
3809
3810
3811 IF (submodelname == 'RTMA') THEN !Start RTMA block
3812
3813!EL field allocation
3814
3815 allocate(eshr(ista_2l:iend_2u,jsta_2l:jend_2u),uvect(ista_2l:iend_2u,jsta_2l:jend_2u),&
3816 vvect(ista_2l:iend_2u,jsta_2l:jend_2u),htsfc(ista_2l:iend_2u,jsta_2l:jend_2u))
3817 allocate(effust(ista_2l:iend_2u,jsta_2l:jend_2u),effvst(ista_2l:iend_2u,jsta_2l:jend_2u),&
3818 esrh(ista_2l:iend_2u,jsta_2l:jend_2u))
3819
3820!
3821 DO j=jsta,jend
3822 DO i=ista,iend
3823 maxthe(i,j)=-h99999
3824 the(i,j)=-h99999
3825 maxthepos(i,j)=0
3826 muq1d(i,j) = 0.
3827 ENDDO
3828 ENDDO
3829 DO l=lm,1,-1
3830
3831 DO j=jsta,jend
3832 DO i=ista,iend
3833 egrid1(i,j) = -h99999
3834 p1d(i,j)=pmid(i,j,l)
3835 t1d(i,j)=t(i,j,l)
3836 q1d(i,j)=q(i,j,l)
3837 ENDDO
3838 ENDDO
3839 CALL calthte(p1d,t1d,q1d,egrid1)
3840 DO j=jsta,jend
3841 DO i=ista,iend
3842 the(i,j)=egrid1(i,j)
3843 IF(the(i,j)>=maxthe(i,j))THEN
3844 maxthe(i,j)=the(i,j)
3845 maxthepos(i,j)=l
3846 muq1d(i,j) = q(i,j,l) ! save the Q of air parcel with max theta-e (MU Parcel)
3847 ENDIF
3848 ENDDO
3849 ENDDO
3850 DO j=jsta,jend
3851 DO i=ista,iend
3852 llow(i,j) = el_base(i,j)
3853 llow_zint(i,j)=zint(i,j,llow(i,j))
3854 ieql_zint(i,j)=zint(i,j,ieql(i,j))
3855 z_temp(i,j)=llow_zint(i,j)+d50*(ieql_zint(i,j)-llow_zint(i,j))
3856 midcal(i,j,l)=abs(zint(i,j,l)-z_temp(i,j))
3857 ENDDO
3858 ENDDO
3859 ENDDO
3860 z_midcal=minloc(midcal,dim=3)
3861
3862!
3863!get surface height
3864 IF(gridtype == 'E')THEN
3865 jvn = 1
3866 jvs = -1
3867 do j=jsta,jend
3868 ive(j) = mod(j,2)
3869 ivw(j) = ive(j)-1
3870 enddo
3871 istart = ista_m
3872 istop = iend_m
3873 jstart = jsta_m
3874 jstop = jend_m
3875 ELSE IF(gridtype == 'B')THEN
3876 jvn = 1
3877 jvs = 0
3878 do j=jsta,jend
3879 ive(j)=1
3880 ivw(j)=0
3881 enddo
3882 istart = ista_m
3883 istop = iend_m
3884 jstart = jsta_m
3885 jstop = jend_m
3886 ELSE
3887 jvn = 0
3888 jvs = 0
3889 do j=jsta,jend
3890 ive(j) = 0
3891 ivw(j) = 0
3892 enddo
3893 istart = ista
3894 istop = iend
3895 jstart = jsta
3896 jstop = jend
3897 END IF
3898
3899 IF(gridtype /= 'A') CALL exch(fis(ista_2l:iend_2u,jsta_2l:jend_2u))
3900 DO j=jstart,jstop
3901 DO i=istart,istop
3902 ie = i+ive(j)
3903 iw = i+ivw(j)
3904 jn = j+jvn
3905 js = j+jvs
3906 IF (gridtype=='B')THEN
3907 htsfc(i,j)=(0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(ie,jn))
3908 ELSE
3909 htsfc(i,j)=(0.25/g)*(fis(iw,j)+fis(ie,j)+fis(i,jn)+fis(i,js))
3910 ENDIF
3911 ENDDO
3912 ENDDO
3913
3914!Height of effbot
3915 IF (iget(979)>0) THEN
3916 grid1=spval
3917 DO j=jsta,jend
3918 DO i=ista,iend
3919 IF(zint(i,j,llow(i,j))<spval.and.htsfc(i,j)<spval)&
3920 grid1(i,j) = zint(i,j,llow(i,j)) - htsfc(i,j)
3921 ENDDO
3922 ENDDO
3923 if(grib=='grib2') then
3924 cfld=cfld+1
3925 fld_info(cfld)%ifld=iavblfld(iget(979))
3926 fld_info(cfld)%lvl=lvlsxml(1,iget(979))
3927!$omp parallel do private(i,j,ii,jj)
3928 do j=1,jend-jsta+1
3929 jj = jsta+j-1
3930 do i=1,iend-ista+1
3931 ii = ista+i-1
3932 datapd(i,j,cfld) = grid1(ii,jj)
3933 enddo
3934 enddo
3935 endif
3936 ENDIF
3937!Height of effbot
3938 IF (iget(980)>0) THEN
3939 grid1=spval
3940 DO j=jsta,jend
3941 DO i=ista,iend
3942 IF(zint(i,j,lupp(i,j))<spval.and.htsfc(i,j)<spval)&
3943 grid1(i,j) = zint(i,j,lupp(i,j)) - htsfc(i,j)
3944 ENDDO
3945 ENDDO
3946 if(grib=='grib2') then
3947 cfld=cfld+1
3948 fld_info(cfld)%ifld=iavblfld(iget(980))
3949 fld_info(cfld)%lvl=lvlsxml(1,iget(980))
3950!$omp parallel do private(i,j,ii,jj)
3951 do j=1,jend-jsta+1
3952 jj = jsta+j-1
3953 do i=1,iend-ista+1
3954 ii = ista+i-1
3955 datapd(i,j,cfld) = grid1(ii,jj)
3956 enddo
3957 enddo
3958 endif
3959 ENDIF
3960
3961!U inflow based to 50% EL shear vector
3962
3963 IF (iget(983)>0) THEN
3964 grid1=spval
3965 DO j=jsta,jend
3966 DO i=ista,iend
3967 IF(llow(i,j)<spval.and.lupp(i,j)<spval) THEN
3968 uvect(i,j)=uh(i,j,z_midcal(i,j))-uh(i,j,llow(i,j))
3969 grid1(i,j)=uvect(i,j)
3970 ENDIF
3971 ENDDO
3972 ENDDO
3973 if(grib=='grib2') then
3974 cfld=cfld+1
3975 fld_info(cfld)%ifld=iavblfld(iget(983))
3976 fld_info(cfld)%lvl=lvlsxml(1,iget(983))
3977!$omp parallel do private(i,j,ii,jj)
3978 do j=1,jend-jsta+1
3979 jj = jsta+j-1
3980 do i=1,iend-ista+1
3981 ii = ista+i-1
3982 datapd(i,j,cfld) = grid1(ii,jj)
3983 enddo
3984 enddo
3985 endif
3986 ENDIF
3987
3988!V inflow based to 50% EL shear vector
3989 IF (iget(984)>0) THEN
3990 grid1=spval
3991 DO j=jsta,jend
3992 DO i=ista,iend
3993 IF(llow(i,j)<spval.and.lupp(i,j)<spval) THEN
3994 vvect(i,j)=vh(i,j,z_midcal(i,j))-vh(i,j,llow(i,j))
3995 grid1(i,j)=vvect(i,j)
3996 ENDIF
3997 ENDDO
3998 ENDDO
3999 if(grib=='grib2') then
4000 cfld=cfld+1
4001 fld_info(cfld)%ifld=iavblfld(iget(984))
4002 fld_info(cfld)%lvl=lvlsxml(1,iget(984))
4003!$omp parallel do private(i,j,ii,jj)
4004 do j=1,jend-jsta+1
4005 jj = jsta+j-1
4006 do i=1,iend-ista+1
4007 ii = ista+i-1
4008 datapd(i,j,cfld) = grid1(ii,jj)
4009 enddo
4010 enddo
4011 endif
4012 ENDIF
4013
4014!Inflow based (ESFC) to (50%) EL shear magnitude
4015 IF (iget(985)>0) THEN
4016 grid1=spval
4017 DO j=jsta,jend
4018 DO i=ista,iend
4019 IF(uvect(i,j)<spval.and.vvect(i,j)<spval) THEN
4020 eshr(i,j)=sqrt((uvect(i,j)**2)+(vvect(i,j))**2)
4021 !effshear
4022 !calc
4023 grid1(i,j)=eshr(i,j) !Effective
4024 ENDIF
4025 ENDDO
4026 ENDDO
4027 if(grib=='grib2') then
4028 cfld=cfld+1
4029 fld_info(cfld)%ifld=iavblfld(iget(985))
4030 fld_info(cfld)%lvl=lvlsxml(1,iget(985))
4031!$omp parallel do private(i,j,ii,jj)
4032 do j=1,jend-jsta+1
4033 jj = jsta+j-1
4034 do i=1,iend-ista+1
4035 ii = ista+i-1
4036 datapd(i,j,cfld) = grid1(ii,jj)
4037 enddo
4038 enddo
4039 endif
4040 ENDIF
4041
4042! Effective Helicity
4043!$omp parallel do private(i,j)
4044 DO j=jsta,jend
4045 DO i=ista,iend
4046 llow(i,j) = el_base(i,j)
4047 lupp(i,j) = el_tops(i,j)
4048 ENDDO
4049 ENDDO
4050
4051 CALL calhel3(llow,lupp,effust,effvst,esrh)
4052
4053!U Bunkers Effective right motion
4054!
4055
4056 IF (iget(986)>0) THEN
4057 grid1=spval
4058 DO j=jsta,jend
4059 DO i=ista,iend
4060 IF(llow(i,j)<spval.and.lupp(i,j)<spval)&
4061 grid1(i,j)=effust(i,j)
4062 ENDDO
4063 ENDDO
4064 if(grib=='grib2') then
4065 cfld=cfld+1
4066 fld_info(cfld)%ifld=iavblfld(iget(986))
4067 fld_info(cfld)%lvl=lvlsxml(1,iget(986))
4068! $omp parallel do private(i,j,ii,jj)
4069 do j=1,jend-jsta+1
4070 jj = jsta+j-1
4071 do i=1,iend-ista+1
4072 ii = ista+i-1
4073 datapd(i,j,cfld) = grid1(ii,jj)
4074 enddo
4075 enddo
4076 endif
4077 ENDIF
4078
4079!V Bunkers Effective right motion
4080 IF (iget(987)>0) THEN
4081 grid1=spval
4082 DO j=jsta,jend
4083 DO i=ista,iend
4084 IF(llow(i,j)<spval.and.lupp(i,j)<spval)&
4085 grid1(i,j)=effvst(i,j)
4086 ENDDO
4087 ENDDO
4088 if(grib=='grib2') then
4089 cfld=cfld+1
4090 fld_info(cfld)%ifld=iavblfld(iget(987))
4091 fld_info(cfld)%lvl=lvlsxml(1,iget(987))
4092! $omp parallel do private(i,j,ii,jj)
4093 do j=1,jend-jsta+1
4094 jj = jsta+j-1
4095 do i=1,iend-ista+1
4096 ii = ista+i-1
4097 datapd(i,j,cfld) = grid1(ii,jj)
4098 enddo
4099 enddo
4100 endif
4101 ENDIF
4102
4103!Effective layer helicity
4104 IF (iget(988)>0) THEN
4105 grid1=spval
4106 DO j=jsta,jend
4107 DO i=ista,iend
4108 IF(llow(i,j)<spval.and.lupp(i,j)<spval)&
4109 grid1(i,j)=esrh(i,j)
4110 ENDDO
4111 ENDDO
4112 if(grib=='grib2') then
4113 cfld=cfld+1
4114 fld_info(cfld)%ifld=iavblfld(iget(988))
4115 fld_info(cfld)%lvl=lvlsxml(1,iget(988))
4116! $omp parallel do private(i,j,ii,jj)
4117 do j=1,jend-jsta+1
4118 jj = jsta+j-1
4119 do i=1,iend-ista+1
4120 ii = ista+i-1
4121 datapd(i,j,cfld) = grid1(ii,jj)
4122 enddo
4123 enddo
4124 endif
4125 ENDIF
4126
4127!Effective Layer Tornado Parameter
4128 IF (iget(989)>0) THEN
4129 DO j=jsta,jend
4130 DO i=ista,iend
4131 IF (mllcl(i,j)>d2000) THEN
4132 mllcltmp=d00
4133 ELSEIF (mllcl(i,j)<d1000) THEN
4134 mllcltmp=1.0
4135 ELSE
4136 mllcltmp=((d2000-mllcl(i,j))/d1000)
4137 ENDIF
4138 IF (eshr(i,j)<12.5) THEN
4139 eshrtmp=d00
4140 ELSEIF (eshr(i,j)>30.0) THEN
4141 eshrtmp=1.5
4142 ELSE
4143 eshrtmp=(eshr(i,j)/20.)
4144 ENDIF
4145 IF (mlcin(i,j)>-50.) THEN
4146 mlcintmp=1.0
4147 ELSEIF (mlcin(i,j)<-200.) THEN
4148 mlcintmp=d00
4149 ELSE
4150 mlcintmp=(200.+mlcin(i,j))/150.
4151 ENDIF
4152 stp=(mlcape(i,j)/d1500)*mllcltmp*(esrh(i,j)/150.)*&
4153 eshrtmp*mlcintmp
4154 grid1(i,j) = spval
4155 IF(llow(i,j)<spval.and.lupp(i,j)<spval) THEN
4156 IF (stp>0) THEN
4157 grid1(i,j)=stp
4158 ELSE
4159 grid1(i,j)=d00
4160 ENDIF
4161 ENDIF
4162 ENDDO
4163 ENDDO
4164 if(grib=='grib2') then
4165 cfld=cfld+1
4166 fld_info(cfld)%ifld=iavblfld(iget(989))
4167 fld_info(cfld)%lvl=lvlsxml(1,iget(989))
4168! $omp parallel do private(i,j,ii,jj)
4169 do j=1,jend-jsta+1
4170 jj = jsta+j-1
4171 do i=1,iend-ista+1
4172 ii = ista+i-1
4173 datapd(i,j,cfld) = grid1(ii,jj)
4174 enddo
4175 enddo
4176 endif
4177 ENDIF
4178
4179!Fixed Layer Tornado Parameter
4180 IF (iget(990)>0) THEN
4181 DO j=jsta,jend
4182 DO i=ista,iend
4183 llmh = nint(lmh(i,j))
4184 p1d(i,j) = pmid(i,j,llmh)
4185 t1d(i,j) = t(i,j,llmh)
4186 q1d(i,j) = q(i,j,llmh)
4187 ENDDO
4188 ENDDO
4189 CALL callcl(p1d,t1d,q1d,egrid1,egrid2)
4190 DO j=jsta,jend
4191 DO i=ista,iend
4192 slcl(i,j)=egrid2(i,j)
4193 ENDDO
4194 ENDDO
4195 itype = 1
4196 dpbnd = 10.e2
4197 dummy = 0.
4198 idummy = 0
4199 CALL calcape(itype,dpbnd,dummy,dummy,dummy,&
4200 idummy,egrid1,egrid2,&
4201 egrid3,dummy,dummy)
4202
4203 DO j=jsta,jend
4204 DO i=ista,iend
4205 IF (slcl(i,j)>d2000) THEN
4206 slcltmp=d00
4207 ELSEIF (slcl(i,j)<=d1000) THEN
4208 slcltmp=1.0
4209 ELSE
4210 slcltmp=((d2000-slcl(i,j))/d1000)
4211 ENDIF
4212 IF (fshr(i,j)<12.5) THEN
4213 fshrtmp=d00
4214 ELSEIF (fshr(i,j)>30.0) THEN
4215 fshrtmp=1.5
4216 ELSE
4217 fshrtmp=(fshr(i,j)/20.)
4218 ENDIF
4219 IF (egrid2(i,j)>-50.) THEN
4220 scintmp=1.0
4221 ELSEIF (egrid2(i,j)<-200.) THEN
4222 scintmp=d00
4223 ELSE
4224 scintmp=((200.+egrid2(i,j)/150.))
4225 ENDIF
4226 stp=(egrid1(i,j)/d1500)*slcltmp*(heli(i,j,2)/150.)*&
4227 fshrtmp*scintmp
4228 grid1(i,j) = spval
4229 IF(t1d(i,j) < spval) THEN
4230 IF (stp>0) THEN
4231 grid1(i,j)=stp
4232 ELSE
4233 grid1(i,j)=d00
4234 ENDIF
4235 ENDIF
4236 ENDDO
4237 ENDDO
4238 if(grib=='grib2') then
4239 cfld=cfld+1
4240 fld_info(cfld)%ifld=iavblfld(iget(990))
4241 fld_info(cfld)%lvl=lvlsxml(1,iget(990))
4242! $omp parallel do private(i,j,ii,jj)
4243 do j=1,jend-jsta+1
4244 jj = jsta+j-1
4245 do i=1,iend-ista+1
4246 ii = ista+i-1
4247 datapd(i,j,cfld) = grid1(ii,jj)
4248 enddo
4249 enddo
4250 endif
4251 ENDIF
4252
4253!Effective Layer Supercell Parameter
4254 IF (iget(991)>0) THEN
4255 DO j=jsta,jend
4256 DO i=ista,iend
4257 IF (eshr(i,j)<10.) THEN
4258 eshrtmp=d00
4259 ELSEIF (eshr(i,j)>20.0) THEN
4260 eshrtmp=1
4261 ELSE
4262 eshrtmp=(eshr(i,j)/20.)
4263 ENDIF
4264 IF (mucin(i,j)>-40.) THEN
4265 mucintmp=1.0
4266 ELSE
4267 mucintmp=(-40./mucin(i,j))
4268 ENDIF
4269 stp=(mucape(i,j)/d1000)*(esrh(i,j)/50.)*&
4270 eshrtmp*mucintmp
4271 grid1(i,j) = spval
4272 IF(t1d(i,j) < spval) THEN
4273 IF (stp>0) THEN
4274 grid1(i,j)=stp
4275 ELSE
4276 grid1(i,j)=d00
4277 ENDIF
4278 ENDIF
4279 ENDDO
4280 ENDDO
4281 if(grib=='grib2') then
4282 cfld=cfld+1
4283 fld_info(cfld)%ifld=iavblfld(iget(991))
4284 fld_info(cfld)%lvl=lvlsxml(1,iget(991))
4285! $omp parallel do private(i,j,ii,jj)
4286 do j=1,jend-jsta+1
4287 jj = jsta+j-1
4288 do i=1,iend-ista+1
4289 ii = ista+i-1
4290 datapd(i,j,cfld) = grid1(ii,jj)
4291 enddo
4292 enddo
4293 endif
4294 ENDIF
4295
4296!Mixed Layer (100 mb) Virtual LFC
4297
4298 IF (iget(992)>0) THEN
4299!$omp parallel do private(i,j)
4300 DO j=jsta,jend
4301 DO i=ista,iend
4302 egrid1(i,j) = -h99999
4303 egrid2(i,j) = -h99999
4304 egrid3(i,j) = -h99999
4305 egrid4(i,j) = -h99999
4306 egrid5(i,j) = -h99999
4307 egrid6(i,j) = -h99999
4308 egrid7(i,j) = -h99999
4309 egrid8(i,j) = -h99999
4310 lb2(i,j) = (lvlbnd(i,j,1) + lvlbnd(i,j,2) + &
4311 lvlbnd(i,j,3))/3
4312 p1d(i,j) = (pbnd(i,j,1) + pbnd(i,j,2) + pbnd(i,j,3))/3
4313 t1d(i,j) = (tvirtual(tbnd(i,j,1),qbnd(i,j,1)) + &
4314 tvirtual(tbnd(i,j,2),qbnd(i,j,2)) + &
4315 tvirtual(tbnd(i,j,3),qbnd(i,j,3)))/3
4316 q1d(i,j) = (qbnd(i,j,1) + qbnd(i,j,2) + qbnd(i,j,3))/3
4317 ENDDO
4318 ENDDO
4319
4320 dpbnd = 0.
4321 itype = 2
4322! EGRID3 is Virtual LFC
4323 CALL calcape2(itype,dpbnd,p1d,t1d,q1d,lb2, &
4324 egrid1,egrid2,egrid3,egrid4,egrid5, &
4325 egrid6,egrid7,egrid8)
4326
4327 grid1=spval
4328 DO j=jsta,jend
4329 DO i=ista,iend
4330 IF(t1d(i,j) < spval) grid1(i,j) = egrid3(i,j)
4331 ENDDO
4332 ENDDO
4333 CALL bound(grid1,d00,h99999)
4334 if(grib=='grib2') then
4335 cfld=cfld+1
4336 fld_info(cfld)%ifld=iavblfld(iget(992))
4337 fld_info(cfld)%lvl=lvlsxml(1,iget(992))
4338!$omp parallel do private(i,j,ii,jj)
4339 do j=1,jend-jsta+1
4340 jj = jsta+j-1
4341 do i=1,iend-ista+1
4342 ii = ista+i-1
4343 datapd(i,j,cfld) = grid1(ii,jj)
4344 enddo
4345 enddo
4346 endif
4347 ENDIF !992
4348
4349
4350 IF (iget(763)>0) THEN
4351!$omp parallel do private(i,j)
4352! EGRID3 is Virtual LFC
4353 DO j=jsta,jend
4354 DO i=ista,iend
4355 grid1(i,j) = q1d(i,j)
4356 ENDDO
4357 ENDDO
4358 if(grib=='grib2') then
4359 cfld=cfld+1
4360 fld_info(cfld)%ifld=iavblfld(iget(763))
4361 fld_info(cfld)%lvl=lvlsxml(1,iget(763))
4362!$omp parallel do private(i,j,ii,jj)
4363 do j=1,jend-jsta+1
4364 jj = jsta+j-1
4365 do i=1,iend-ista+1
4366 ii = ista+i-1
4367 datapd(i,j,cfld) = grid1(ii,jj)
4368 enddo
4369 enddo
4370 endif
4371 ENDIF
4372
4373!Hail parameter
4374 IF (iget(993)>0) THEN
4375 grid1=spval
4376 DO j=jsta,jend
4377 DO i=ista,iend
4378 lapse=-((t700(i,j)-t500(i,j))/((z700(i,j)-z500(i,j))))
4379 ship=(mucape(i,j)*d1000*muq1d(i,j)*lapse*(t500(i,j)-k2c)*fshr(i,j))/hconst
4380 IF (mucape(i,j)<1300.)THEN
4381 ship=ship*(mucape(i,j)/1300.)
4382 ENDIF
4383 IF (lapse < 5.8)THEN
4384 ship=ship*(lapse/5.8)
4385 ENDIF
4386 IF (freezelvl(i,j) < 2400.)THEN
4387 ship=ship*(freezelvl(i,j)/2400.)
4388 ENDIF
4389 grid1(i,j)=ship
4390 ENDDO
4391 ENDDO
4392 if(grib=='grib2') then
4393 cfld=cfld+1
4394 fld_info(cfld)%ifld=iavblfld(iget(993))
4395 fld_info(cfld)%lvl=lvlsxml(1,iget(993))
4396! $omp parallel do private(i,j,ii,jj)
4397 do j=1,jend-jsta+1
4398 jj = jsta+j-1
4399 do i=1,iend-ista+1
4400 ii = ista+i-1
4401 datapd(i,j,cfld) = grid1(ii,jj)
4402 enddo
4403 enddo
4404 endif
4405 ENDIF
4406
4407 ENDIF !END RTMA BLOCK
4408
4409
4410! Critical Angle
4411
4412 IF (iget(957)>0) THEN
4413 grid1=spval
4414!$omp parallel do private(i,j)
4415 DO j=jsta,jend
4416 DO i=ista,iend
4417 IF(t1d(i,j) < spval ) grid1(i,j) = cangle(i,j)
4418 ! IF(EGRID1(I,J)<100. .OR. EGRID2(I,J)>-250.) THEN
4419 ! GRID1(I,J) = 0.
4420 ! ENDIF
4421 ENDDO
4422 ENDDO
4423 if(grib=='grib2') then
4424 cfld=cfld+1
4425 fld_info(cfld)%ifld=iavblfld(iget(957))
4426 fld_info(cfld)%lvl=lvlsxml(1,iget(957))
4427!$omp parallel do private(i,j,ii,jj)
4428 do j=1,jend-jsta+1
4429 jj = jsta+j-1
4430 do i=1,iend-ista+1
4431 ii = ista+i-1
4432 datapd(i,j,cfld) = grid1(ii,jj)
4433 enddo
4434 enddo
4435 endif
4436 ENDIF !957
4437
4438! Dendritic Layer Depth, -17C < T < -12C
4439
4440 IF (iget(955)>0) THEN
4441 grid1=spval
4442!$omp parallel do private(i,j)
4443 DO j=jsta,jend
4444 DO i=ista,iend
4445 IF(t1d(i,j) < spval ) grid1(i,j) = egrid7(i,j)
4446 ENDDO
4447 ENDDO
4448 CALL bound(grid1,d00,h99999)
4449 if(grib=='grib2') then
4450 cfld=cfld+1
4451 fld_info(cfld)%ifld=iavblfld(iget(955))
4452 fld_info(cfld)%lvl=lvlsxml(1,iget(955))
4453!$omp parallel do private(i,j,ii,jj)
4454 do j=1,jend-jsta+1
4455 jj = jsta+j-1
4456 do i=1,iend-ista+1
4457 ii = ista+i-1
4458 datapd(i,j,cfld) = grid1(ii,jj)
4459 enddo
4460 enddo
4461 endif
4462 ENDIF !955
4463
4464! Enhanced Stretching Potential
4465
4466 IF (iget(956)>0) THEN
4467 grid1=spval
4468!$omp parallel do private(i,j)
4469 DO j=jsta,jend
4470 DO i=ista,iend
4471 IF(t1d(i,j) < spval ) grid1(i,j) = egrid8(i,j)
4472 ENDDO
4473 ENDDO
4474 CALL bound(grid1,d00,h99999)
4475 if(grib=='grib2') then
4476 cfld=cfld+1
4477 fld_info(cfld)%ifld=iavblfld(iget(956))
4478 fld_info(cfld)%lvl=lvlsxml(1,iget(956))
4479!$omp parallel do private(i,j,ii,jj)
4480 do j=1,jend-jsta+1
4481 jj = jsta+j-1
4482 do i=1,iend-ista+1
4483 ii = ista+i-1
4484 datapd(i,j,cfld) = grid1(ii,jj)
4485 enddo
4486 enddo
4487 endif
4488 ENDIF !956
4489
4490! Downdraft CAPE
4491
4492! ITYPE = 1
4493! DO J=JSTA,JEND
4494! DO I=ISTA,IEND
4495! LB2(I,J) = (LVLBND(I,J,1) + LVLBND(I,J,2) + &
4496! LVLBND(I,J,3))/3
4497! P1D(I,J) = (PBND(I,J,1) + PBND(I,J,2) + PBND(I,J,3))/3
4498! T1D(I,J) = (TBND(I,J,1) + TBND(I,J,2) + TBND(I,J,3))/3
4499! Q1D(I,J) = (QBND(I,J,1) + QBND(I,J,2) + QBND(I,J,3))/3
4500! ENDDO
4501! ENDDO
4502
4503! DPBND = 400.E2
4504! CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2, &
4505! EGRID1,EGRID2,EGRID3,EGRID4,EGRID5, &
4506! EGRID6,EGRID7,EGRID8)
4507
4508 IF (iget(954)>0) THEN
4509 grid1 = spval
4510!$omp parallel do private(i,j)
4511 DO j=jsta,jend
4512 DO i=ista,iend
4513 IF(t1d(i,j) < spval) grid1(i,j) = -egrid6(i,j)
4514 ENDDO
4515 ENDDO
4516 CALL bound(grid1,d00,h99999)
4517 if(grib=='grib2') then
4518 cfld=cfld+1
4519 fld_info(cfld)%ifld=iavblfld(iget(954))
4520 fld_info(cfld)%lvl=lvlsxml(1,iget(954))
4521!$omp parallel do private(i,j,ii,jj)
4522 do j=1,jend-jsta+1
4523 jj = jsta+j-1
4524 do i=1,iend-ista+1
4525 ii = ista+i-1
4526 datapd(i,j,cfld) = grid1(ii,jj)
4527 enddo
4528 enddo
4529 endif
4530
4531 ENDIF !954
4532
4533 if (allocated(ushr1)) deallocate(ushr1)
4534 if (allocated(vshr1)) deallocate(vshr1)
4535 if (allocated(ushr6)) deallocate(ushr6)
4536 if (allocated(vshr6)) deallocate(vshr6)
4537 if (allocated(ust)) deallocate(ust)
4538 if (allocated(vst)) deallocate(vst)
4539 if (allocated(heli)) deallocate(heli)
4540 if (allocated(llow)) deallocate(llow)
4541 if (allocated(lupp)) deallocate(lupp)
4542 if (allocated(cangle))deallocate(cangle)
4543 if (allocated(effust))deallocate(effust)
4544 if (allocated(effvst))deallocate(effvst)
4545 if (allocated(eshr)) deallocate(eshr)
4546 if (allocated(uvect)) deallocate(uvect)
4547 if (allocated(vvect)) deallocate(vvect)
4548 if (allocated(esrh)) deallocate(esrh)
4549 if (allocated(htsfc)) deallocate(htsfc)
4550 if (allocated(fshr)) deallocate(fshr)
4551 if (allocated(llow_zint)) deallocate(llow_zint)
4552 if (allocated(ieql_zint)) deallocate(ieql_zint)
4553 if (allocated(z_temp)) deallocate(z_temp)
4554 if (allocated(midcal)) deallocate(midcal)
4555 if (allocated(z_midcal)) deallocate(z_midcal)
4556 if (allocated(el_base)) deallocate(el_base)
4557 if (allocated(el_tops)) deallocate(el_tops)
4558
4559 ENDIF
4560
4561 if (allocated(pbnd)) deallocate(pbnd)
4562 if (allocated(tbnd)) deallocate(tbnd)
4563 if (allocated(qbnd)) deallocate(qbnd)
4564 if (allocated(ubnd)) deallocate(ubnd)
4565 if (allocated(vbnd)) deallocate(vbnd)
4566 if (allocated(rhbnd)) deallocate(rhbnd)
4567 if (allocated(wbnd)) deallocate(wbnd)
4568 if (allocated(lvlbnd)) deallocate(lvlbnd)
4569 if (allocated(lb2)) deallocate(lb2)
4570
4571!
4572!
4573! RELATIVE HUMIDITY WITH RESPECT TO PRECIPITABLE WATER
4574 IF (iget(749)>0) THEN
4575 CALL calrh_pw(grid1(ista:iend,jsta:jend))
4576 if(grib=='grib2') then
4577 cfld=cfld+1
4578 fld_info(cfld)%ifld=iavblfld(iget(749))
4579!$omp parallel do private(i,j,ii,jj)
4580 do j=1,jend-jsta+1
4581 jj = jsta+j-1
4582 do i=1,iend-ista+1
4583 ii = ista+i-1
4584 datapd(i,j,cfld) = grid1(ii,jj)
4585 enddo
4586 enddo
4587 endif
4588 ENDIF
4589
4590
4591! END OF ROUTINE.
4592!
4593 RETURN
4594 CONTAINS
4595
4596 subroutine allocate_cape_arrays
4597 if(.not.allocated(omgbnd)) allocate(omgbnd(ista:iend,jsta:jend,nbnd))
4598 if(.not.allocated(pwtbnd)) allocate(pwtbnd(ista:iend,jsta:jend,nbnd))
4599 if(.not.allocated(qcnvbnd)) allocate(qcnvbnd(ista:iend,jsta:jend,nbnd))
4600 if(.not.allocated(lvlbnd)) allocate(lvlbnd(ista:iend,jsta:jend,nbnd))
4601 if(.not.allocated(lb2)) allocate(lb2(ista:iend,jsta:jend))
4602 end subroutine allocate_cape_arrays
4603
4604 END SUBROUTINE miscln
bndlyr
subroutine bndlyr(pbnd, tbnd, qbnd, rhbnd, ubnd, vbnd, wbnd, omgbnd, pwtbnd, qcnvbnd, lvlbnd)
Computes constant mass mean fields.
Definition BNDLYR.f:50
bound
subroutine bound(fld, fmin, fmax)
This routine bounds data in the passed array FLD (im x jm elements long) and clips data values such t...
Definition BOUND.f:35
caldwp
subroutine caldwp(p1d, q1d, tdwp, t1d)
Computes dewpoint from P, T, and Q.
Definition CALDWP.f:28
calhel2
subroutine calhel2(llow, lupp, depth, ust, vst, heli, cangle)
Subroutine that computes storm relative helicity.
Definition CALHEL2.f:58
calhel3
subroutine calhel3(llow, lupp, ust, vst, heli)
Subroutine that computes storm relative helicity.
Definition CALHEL3.f:57
calhel
subroutine calhel(depth, ust, vst, heli, ushr1, vshr1, ushr6, vshr6)
This routine computes estimated storm motion and storm-relative environmental helicity.
Definition CALHEL.f:44
callcl
subroutine callcl(p1d, t1d, q1d, plcl, zlcl)
Subroutine that computes the lifting condensation level (LCL) height (above ground level) and pressur...
Definition CALLCL.f:42
calpot
subroutine calpot(p1d, t1d, theta)
Subroutine that computes potential temperature.
Definition CALPOT.f:29
calthte
subroutine calthte(p1d, t1d, q1d, thte)
Subroutine that computes Theta-E.
Definition CALTHTE.f:32
calupdhel
subroutine calupdhel(updhel)
Subroutine that computes the updraft helicity.
Definition CALUPDHEL.f:24
fdlvl_mass
subroutine fdlvl_mass(itype, nfd, ptfd, htfd, nin, qin, qtype, qfd)
Computes FD level for mass variables.
Definition FDLVL.f:785
frzlvl2
subroutine frzlvl2(isotherm, zfrz, rhfrz, pfrzl)
Subroutine that computes FRZING LVL, Z and RH.
Definition FRZLVL2.f:56
frzlvl
subroutine frzlvl(zfrz, rhfrz, pfrzl)
FRZLVL() Subroutine that computes FRZING LVL, Z and RH.
Definition FRZLVL.f:49
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition GFSPOST.F:460
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition GFSPOST.F:373
lfmfld
subroutine lfmfld(rh3310, rh6610, rh3366, pw3310)
LFMFLD() computes layer mean LFM fields.
Definition LFMFLD.f:53
lfmfld_gfs
subroutine lfmfld_gfs(rh4410, rh7294, rh4472, rh3310)
LFMFLD_GFS() computes layer mean LFM fields.
Definition LFMFLD_GFS.f:55
miscln
subroutine miscln
SUBPROGRAM: MISCLN POSTS MISCELLANEOUS FIELDS PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-20.
Definition MISCLN.f:90
ngmfld
subroutine ngmfld(rh4710, rh4796, rh1847, rh8498, qm8510)
ngmfld() computes layer mean NGM fields
Definition NGMFLD.f:60
otlft
subroutine otlft(pbnd, tbnd, qbnd, slindx)
otlft() computes lifted index.
Definition OTLFT.f:45
sclfld
subroutine sclfld(fld, scale, imo, jmo)
sclfld() scale array element by constant.
Definition SCLFLD.f:34