UPP  V11.0.0
 All Data Structures Files Functions Pages
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
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition: GFSPOST.F:431
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
grib2_module
Definition: grib2_module.f:1
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:345
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition: GFSPOST.F:344
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27