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