UPP  V11.0.0
 All Data Structures Files Functions Pages
SURFCE.f
Go to the documentation of this file.
1 
2 !
63  SUBROUTINE surfce
64 
65 !
66 !
67 ! INCLUDE GRID DIMENSIONS. SET/DERIVE OTHER PARAMETERS.
68 !
69  use vrbls4d, only: smoke
70  use vrbls3d, only: zint, pint, t, pmid, q, f_rimef
71  use vrbls2d, only: ths, qs, qvg, qv2m, tsnow, tg, smstav, smstot, &
72  cmc, sno, snoavg, psfcavg, t10avg, snonc, ivgtyp, &
73  si, potevp, dzice, qwbs, vegfrc, isltyp, pshltr, &
74  tshltr, qshltr, mrshltr, maxtshltr, mintshltr, &
75  maxrhshltr, minrhshltr, u10, psfcavg, v10, u10max, &
76  v10max, th10, t10m, q10, wspd10max, &
77  wspd10umax, wspd10vmax, prec, sr, &
78  cprate, avgcprate, avgprec, acprec, cuprec, ancprc, &
79  lspa, acsnow, acsnom, snowfall,ssroff, bgroff, &
80  runoff, pcp_bucket, rainnc_bucket, snow_bucket, &
81  snownc, tmax, graup_bucket, graupelnc, qrmax, sfclhx,&
82  rainc_bucket, sfcshx, subshx, snopcx, sfcuvx, &
83  sfcvx, smcwlt, suntime, pd, sfcux, sfcuxi, sfcvxi, sfcevp, z0, &
84  ustar, mdltaux, mdltauy, gtaux, gtauy, twbs, &
85  sfcexc, grnflx, islope, czmean, czen, rswin,akhsavg ,&
86  akmsavg, u10h, v10h,snfden,sndepac,qvl1, &
87  spduv10mean,swradmean,swnormmean,prate_max,fprate_max &
88  ,fieldcapa,edir,ecan,etrans,esnow,u10mean,v10mean, &
89  avgedir,avgecan,avgetrans,avgesnow,acgraup,acfrain, &
90  acond,maxqshltr,minqshltr,avgpotevp,avgprec_cont, &
91  avgcprate_cont,sst,pcp_bucket1,rainnc_bucket1, &
92  snow_bucket1, rainc_bucket1, graup_bucket1, &
93  shdmin, shdmax, lai, ch10,cd10,landfrac,paha,pahi, &
94  tecan,tetran,tedir,twa
95  use soil, only: stc, sllevel, sldpth, smc, sh2o
96  use masks, only: lmh, sm, sice, htm, gdlat, gdlon
97  use physcons_post,only: con_eps, con_epsm1
98  use params_mod, only: p1000, capa, h1m12, pq0, a2,a3, a4, h1, d00, d01,&
99  eps, oneps, d001, h99999, h100, small, h10e5, &
100  elocp, g, xlai, tfrz, rd
101  use ctlblk_mod, only: jsta, jend, lm, spval, grib, cfld, fld_info, &
102  datapd, nsoil, isf_surface_physics, tprec, ifmin,&
103  modelname, tmaxmin, pthresh, dtq2, dt, nphs, &
104  ifhr, prec_acc_dt, sdat, ihrst, jsta_2l, jend_2u,&
105  lp1, imp_physics, me, asrfc, tsrfc, pt, pdtop, &
106  mpi_comm_comp, im, jm, prec_acc_dt1, &
107  ista, iend, ista_2l, iend_2u
108  use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
109  use grib2_module, only: read_grib2_head, read_grib2_sngle
110  use upp_physics, only: fpvsnew, calrh
111 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112  implicit none
113 !
114  include "mpif.h"
115 !
116 ! IN NGM SUBROUTINE OUTPUT WE FIND THE FOLLOWING COMMENT.
117 ! "IF THE FOLLOWING THRESHOLD VALUES ARE CHANGED, CONTACT
118 ! TDL/SYNOPTIC-SCALE TECHNIQUES BRANCH (PAUL DALLAVALLE
119 ! AND JOHN JENSENIUS). THEY MAY BE USING IT IN ONE OF
120 ! THEIR PACKING CODES." THE THRESHOLD VALUE IS 0.01 INCH
121 ! OR 2.54E-4 METER. PRECIPITATION VALUES LESS THAN THIS
122 ! THRESHOLD ARE SET TO MINUS ONE TIMES THIS THRESHOLD.
123  real,PARAMETER :: ptrace = 0.000254e0
124 !
125 ! SET CELCIUS TO KELVIN AND SECOND TO HOUR CONVERSION.
126  integer,parameter :: nalg=5, nosoiltype=9
127  real, PARAMETER :: c2k = 273.15, sec2hr = 1./3600.
128 !
129 ! DECLARE VARIABLES.
130 !
131  integer, dimension(ista:iend,jsta:jend) :: nroots, iwx1
132  real, allocatable, dimension(:,:) :: zsfc, psfc, tsfc, qsfc, &
133  rhsfc, thsfc, dwpsfc, p1d, &
134  t1d, q1d, zwet, &
135  smcdry, smcmax,doms, domr, &
136  domip, domzr, rsmin, smcref,&
137  rcq, rct, rcsoil, gc, rcs
138 
139  real, dimension(ista:iend,jsta:jend) :: evp
140  real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: egrid1, egrid2
141  real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: grid2
142  real, dimension(im,jm) :: grid1
143  real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: iceg
144 ! , ua, va
145  real, allocatable, dimension(:,:,:) :: sleet, rain, freezr, snow
146 ! real, dimension(im,jm,nalg) :: sleet, rain, freezr, snow
147 
148 !GSD
149  REAL totprcp, snowratio,t2,rainl
150 
151 !
152  integer i,j,iwx,itmaxmin,ifincr,isvalue,ii,jj, &
153  itprec,itsrfc,l,ls,iveg,llmh, &
154  ivg,irtn,iseed, icat, cnt_snowratio(10),icnt_snow_rain_mixed
155 
156  real rdtphs,tlow,tsfck,qsat,dtop,dbot,sneqv,rrnum,sfcprs,sfcq, &
157  rc,sfctmp,sncovr,factrs,solar, s,tk,tl,w,t2c,dlt,ape, &
158  qv,e,dwpt,dum1,dum2,dum3,dum1s,dum3s,dum21,dum216,es
159 
160  character(len=256) :: ffgfile
161  character(len=256) :: arifile
162 
163  logical file_exists
164 
165  logical, parameter :: debugprint = .false.
166 
167 
168 !****************************************************************************
169 !
170 ! START SURFCE.
171 !
172 !
173 !*** BLOCK 1. SURFACE BASED FIELDS.
174 !
175 ! IF ANY OF THE FOLLOWING "SURFACE" FIELDS ARE REQUESTED,
176 ! WE NEED TO COMPUTE THE FIELDS FIRST.
177 !
178  IF ( (iget(024)>0).OR.(iget(025)>0).OR. &
179  (iget(026)>0).OR.(iget(027)>0).OR. &
180  (iget(028)>0).OR.(iget(029)>0).OR. &
181  (iget(154)>0).OR. &
182  (iget(034)>0).OR.(iget(076)>0) ) THEN
183 !
184  allocate(zsfc(ista:iend,jsta:jend), psfc(ista:iend,jsta:jend), tsfc(ista:iend,jsta:jend)&
185  ,rhsfc(ista:iend,jsta:jend), thsfc(ista:iend,jsta:jend), qsfc(ista:iend,jsta:jend))
186 !$omp parallel do private(i,j,tsfck,qsat,es)
187  DO j=jsta,jend
188  DO i=ista,iend
189 !
190 ! SCALE ARRAY FIS BY GI TO GET SURFACE HEIGHT.
191 ! ZSFC(I,J)=FIS(I,J)*GI
192 
193 ! dong add missing value for zsfc
194  zsfc(i,j) = spval
195  IF(zint(i,j,lm+1) < spval) &
196  zsfc(i,j) = zint(i,j,lm+1)
197  psfc(i,j) = pint(i,j,nint(lmh(i,j))+1) ! SURFACE PRESSURE.
198 !
199 ! SURFACE (SKIN) POTENTIAL TEMPERATURE AND TEMPERATURE.
200  thsfc(i,j) = ths(i,j)
201  tsfc(i,j) = spval
202  IF(thsfc(i,j) /= spval .and. psfc(i,j) /= spval) &
203  tsfc(i,j) = thsfc(i,j)*(psfc(i,j)/p1000)**capa
204 !
205 ! SURFACE SPECIFIC HUMIDITY, RELATIVE HUMIDITY, AND DEWPOINT.
206 ! ADJUST SPECIFIC HUMIDITY IF RELATIVE HUMIDITY EXCEEDS 0.1 OR 1.0.
207 
208 ! dong spfh sfc set missing value
209  qsfc(i,j) = spval
210  rhsfc(i,j) = spval
211  evp(i,j) = spval
212  IF(tsfc(i,j) < spval) then
213  IF(qs(i,j)<spval) qsfc(i,j) = max(h1m12,qs(i,j))
214  tsfck = tsfc(i,j)
215 
216  IF(modelname == 'RAPR') THEN
217  qsat = max(0.0001,pq0/psfc(i,j)*exp(a2*(tsfck-a3)/(tsfck-a4)))
218  elseif (modelname == 'GFS') then
219  es = fpvsnew(tsfck)
220  qsat = con_eps*es/(psfc(i,j)+con_epsm1*es)
221  ELSE
222  qsat = pq0/psfc(i,j)*exp(a2*(tsfck-a3)/(tsfck-a4))
223  ENDIF
224  rhsfc(i,j) = max(d01, min(h1,qsfc(i,j)/qsat))
225 
226  qsfc(i,j) = rhsfc(i,j)*qsat
227  rhsfc(i,j) = rhsfc(i,j) * 100.0
228  evp(i,j) = d001*psfc(i,j)*qsfc(i,j)/(eps+oneps*qsfc(i,j))
229  END IF !end TSFC
230 !
231 !mp ACCUMULATED NON-CONVECTIVE PRECIP.
232 !mp IF(IGET(034)>0)THEN
233 !mp IF(LVLS(1,IGET(034))>0)THEN
234 
235 ! ACCUMULATED PRECIP (convective + non-convective)
236 ! IF(IGET(087) > 0)THEN
237 ! IF(LVLS(1,IGET(087)) > 0)THEN
238 ! write(6,*) 'acprec, ancprc, cuprec: ', ANCPRC(I,J)+CUPREC(I,J),
239 ! + ANCPRC(I,J),CUPREC(I,J)
240 ! ACPREC(I,J) = ANCPRC(I,J) + CUPREC(I,J) ???????
241 ! ENDIF
242 ! ENDIF
243 
244  ENDDO
245  ENDDO
246 !
247 ! INTERPOLATE/OUTPUT REQUESTED SURFACE FIELDS.
248 !
249 ! SURFACE PRESSURE.
250  IF (iget(024)>0) THEN
251  if(grib == 'grib2') then
252  cfld = cfld+1
253  fld_info(cfld)%ifld = iavblfld(iget(024))
254 !$omp parallel do private(i,j,ii,jj)
255  do j=1,jend-jsta+1
256  jj = jsta+j-1
257  do i=1,iend-ista+1
258  ii = ista+i-1
259  datapd(i,j,cfld) = psfc(ii,jj)
260  enddo
261  enddo
262  endif
263  ENDIF
264 !
265 ! SURFACE HEIGHT.
266  IF (iget(025)>0) THEN
267 !! CALL BOUND(GRID1,D00,H99999)
268  if(grib == 'grib2') then
269  cfld=cfld+1
270  fld_info(cfld)%ifld = iavblfld(iget(025))
271 !$omp parallel do private(i,j,ii,jj)
272  do j=1,jend-jsta+1
273  jj = jsta+j-1
274  do i=1,iend-ista+1
275  ii = ista+i-1
276  datapd(i,j,cfld) = zsfc(ii,jj)
277  enddo
278  enddo
279  endif
280  ENDIF
281  if (allocated(zsfc)) deallocate(zsfc)
282  if (allocated(psfc)) deallocate(psfc)
283 !
284 ! SURFACE (SKIN) TEMPERATURE.
285  IF (iget(026)>0) THEN
286  if(grib == 'grib2') then
287  cfld = cfld+1
288  fld_info(cfld)%ifld = iavblfld(iget(026))
289 !$omp parallel do private(i,j,ii,jj)
290  do j=1,jend-jsta+1
291  jj = jsta+j-1
292  do i=1,iend-ista+1
293  ii = ista+i-1
294  datapd(i,j,cfld) = tsfc(ii,jj)
295  enddo
296  enddo
297  endif
298  ENDIF
299  if (allocated(tsfc)) deallocate(tsfc)
300 !
301 ! SURFACE (SKIN) POTENTIAL TEMPERATURE.
302  IF (iget(027)>0) THEN
303  if(grib=='grib2') then
304  cfld=cfld+1
305  fld_info(cfld)%ifld=iavblfld(iget(027))
306 !$omp parallel do private(i,j,ii,jj)
307  do j=1,jend-jsta+1
308  jj = jsta+j-1
309  do i=1,iend-ista+1
310  ii = ista+i-1
311  datapd(i,j,cfld) = thsfc(ii,jj)
312  enddo
313  enddo
314  endif
315  ENDIF
316  if (allocated(thsfc)) deallocate(thsfc)
317 !
318 ! SURFACE SPECIFIC HUMIDITY.
319  IF (iget(028)>0) THEN
320  !CALL BOUND(GRID1,H1M12,H99999)
321  if(grib=='grib2') then
322  cfld=cfld+1
323  fld_info(cfld)%ifld=iavblfld(iget(028))
324 !$omp parallel do private(i,j,ii,jj)
325  do j=1,jend-jsta+1
326  jj = jsta+j-1
327  do i=1,iend-ista+1
328  ii = ista+i-1
329  datapd(i,j,cfld) = qsfc(ii,jj)
330  enddo
331  enddo
332  endif
333  ENDIF
334  if (allocated(qsfc)) deallocate(qsfc)
335 !
336 ! SURFACE DEWPOINT TEMPERATURE.
337  IF (iget(029)>0) THEN
338  allocate(dwpsfc(ista:iend,jsta:jend))
339  CALL dewpoint(evp,dwpsfc)
340  if(grib=='grib2') then
341  cfld=cfld+1
342  fld_info(cfld)%ifld=iavblfld(iget(029))
343 !$omp parallel do private(i,j,ii,jj)
344  do j=1,jend-jsta+1
345  jj = jsta+j-1
346  do i=1,iend-ista+1
347  ii = ista+i-1
348  datapd(i,j,cfld) = dwpsfc(ii,jj)
349  enddo
350  enddo
351  endif
352  if (allocated(dwpsfc)) deallocate(dwpsfc)
353  ENDIF
354 !
355 ! SURFACE RELATIVE HUMIDITY.
356  IF (iget(076)>0) THEN
357  CALL bound(rhsfc,h1,h100)
358  if(grib=='grib2') then
359  cfld=cfld+1
360  fld_info(cfld)%ifld=iavblfld(iget(076))
361 !$omp parallel do private(i,j,ii,jj)
362  do j=1,jend-jsta+1
363  jj = jsta+j-1
364  do i=1,iend-ista+1
365  ii = ista+i-1
366  datapd(i,j,cfld) = rhsfc(ii,jj)
367  enddo
368  enddo
369  endif
370  ENDIF
371  if (allocated(rhsfc)) deallocate(rhsfc)
372 !
373  ENDIF
374 
375 ! ADDITIONAL SURFACE-SOIL LEVEL FIELDS.
376 !
377 ! SURFACE MIXING RATIO
378  IF (iget(762)>0) THEN
379  if(grib=='grib2') then
380  cfld=cfld+1
381  fld_info(cfld)%ifld=iavblfld(iget(762))
382 !$omp parallel do private(i,j,ii,jj)
383  do j=1,jend-jsta+1
384  jj = jsta+j-1
385  do i=1,iend-ista+1
386  ii = ista+i-1
387  datapd(i,j,cfld) = qvg(ii,jj)
388  enddo
389  enddo
390  endif
391  ENDIF
392 !
393 
394 ! SHELTER MIXING RATIO
395  IF (iget(760)>0) THEN
396  if(grib=='grib2') then
397  cfld=cfld+1
398  fld_info(cfld)%ifld=iavblfld(iget(760))
399 !$omp parallel do private(i,j,ii,jj)
400  do j=1,jend-jsta+1
401  jj = jsta+j-1
402  do i=1,iend-ista+1
403  ii = ista+i-1
404  datapd(i,j,cfld) = qv2m(ii,jj)
405  enddo
406  enddo
407  endif
408  ENDIF
409 
410 ! SNOW TEMERATURE
411  IF (iget(761)>0) THEN
412  if(grib=='grib2') then
413  cfld=cfld+1
414  fld_info(cfld)%ifld=iavblfld(iget(761))
415 !$omp parallel do private(i,j,ii,jj)
416  do j=1,jend-jsta+1
417  jj = jsta+j-1
418  do i=1,iend-ista+1
419  ii = ista+i-1
420  datapd(i,j,cfld) = tsnow(ii,jj)
421  enddo
422  enddo
423  endif
424  ENDIF
425 
426 ! DENSITY OF SNOWFALL
427  IF (iget(724)>0) THEN
428  if(grib=='grib2') then
429  cfld=cfld+1
430  fld_info(cfld)%ifld=iavblfld(iget(724))
431 !$omp parallel do private(i,j,ii,jj)
432  do j=1,jend-jsta+1
433  jj = jsta+j-1
434  do i=1,iend-ista+1
435  ii = ista+i-1
436  datapd(i,j,cfld) = snfden(ii,jj)
437  enddo
438  enddo
439  endif
440  ENDIF
441 
442 ! ACCUMULATED DEPTH OF SNOWFALL
443  IF (iget(725)>0) THEN
444  id(1:25) = 0
445  itprec = nint(tprec)
446 !mp
447  IF(itprec /= 0) THEN
448  ifincr = mod(ifhr,itprec)
449  IF(ifmin >= 1)ifincr = mod(ifhr*60+ifmin,itprec*60)
450  ELSE
451  ifincr = 0
452  ENDIF
453 !mp
454  id(18) = 0
455  id(19) = ifhr
456  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
457  id(20) = 4
458  IF (ifincr==0) THEN
459  id(18) = ifhr-itprec
460  ELSE
461  id(18) = ifhr-ifincr
462  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
463  ENDIF
464  IF (id(18)<0) id(18) = 0
465  if(grib=='grib2') then
466  cfld=cfld+1
467  fld_info(cfld)%ifld=iavblfld(iget(725))
468  fld_info(cfld)%ntrange=1
469  fld_info(cfld)%tinvstat=ifhr-id(18)
470 !$omp parallel do private(i,j,ii,jj)
471  do j=1,jend-jsta+1
472  jj = jsta+j-1
473  do i=1,iend-ista+1
474  ii = ista+i-1
475  datapd(i,j,cfld) = sndepac(ii,jj)
476  enddo
477  enddo
478  endif
479  ENDIF
480 
481 !
482 ! ADDITIONAL SURFACE-SOIL LEVEL FIELDS.
483 !
484 ! print *,'in surf,nsoil=',nsoil,'iget(116)=',iget(116), &
485 ! 'lvls(116)=',LVLS(1:4,IGET(116)), &
486 ! 'sf_sfc_phys=',iSF_SURFACE_PHYSICS
487 
488  DO l=1,nsoil
489 ! SOIL TEMPERATURE.
490  IF (iget(116)>0) THEN
491  IF (lvls(l,iget(116))>0) THEN
492  IF(isf_surface_physics==3)THEN
493  if(grib=='grib2') then
494  cfld=cfld+1
495  fld_info(cfld)%ifld=iavblfld(iget(116))
496  fld_info(cfld)%lvl=lvlsxml(l,iget(116))
497 !$omp parallel do private(i,j,ii,jj)
498  do j=1,jend-jsta+1
499  jj = jsta+j-1
500  do i=1,iend-ista+1
501  ii = ista+i-1
502  datapd(i,j,cfld) = stc(ii,jj,l)
503  enddo
504  enddo
505  endif
506 
507  ELSE
508 
509  dtop = 0.
510  DO ls=1,l-1
511  dtop = dtop + sldpth(ls)
512  ENDDO
513  dbot = dtop + sldpth(l)
514  if(grib=='grib2') then
515  cfld=cfld+1
516  fld_info(cfld)%ifld=iavblfld(iget(116))
517  fld_info(cfld)%lvl=lvlsxml(l,iget(116))
518 !$omp parallel do private(i,j,ii,jj)
519  do j=1,jend-jsta+1
520  jj = jsta+j-1
521  do i=1,iend-ista+1
522  ii = ista+i-1
523  datapd(i,j,cfld) = stc(ii,jj,l)
524  enddo
525  enddo
526  endif
527 
528  ENDIF
529  ENDIF
530  ENDIF
531 !
532 ! SOIL MOISTURE.
533  IF (iget(117)>0) THEN
534  IF (lvls(l,iget(117))>0) THEN
535  IF(isf_surface_physics==3)THEN
536  if(grib=='grib2') then
537  cfld=cfld+1
538  fld_info(cfld)%ifld=iavblfld(iget(117))
539  fld_info(cfld)%lvl=lvlsxml(l,iget(117))
540 !$omp parallel do private(i,j,ii,jj)
541  do j=1,jend-jsta+1
542  jj = jsta+j-1
543  do i=1,iend-ista+1
544  ii = ista+i-1
545  datapd(i,j,cfld) = smc(ii,jj,l)
546  enddo
547  enddo
548  endif
549  ELSE
550  dtop = 0.
551  DO ls=1,l-1
552  dtop = dtop + sldpth(ls)
553  ENDDO
554  dbot = dtop + sldpth(l)
555  if(grib=='grib2') then
556  cfld=cfld+1
557  fld_info(cfld)%ifld=iavblfld(iget(117))
558  fld_info(cfld)%lvl=lvlsxml(l,iget(117))
559 !$omp parallel do private(i,j,ii,jj)
560  do j=1,jend-jsta+1
561  jj = jsta+j-1
562  do i=1,iend-ista+1
563  ii = ista+i-1
564  datapd(i,j,cfld) = smc(ii,jj,l)
565  enddo
566  enddo
567  endif
568  ENDIF
569  ENDIF
570  ENDIF
571 ! ADD LIQUID SOIL MOISTURE
572  IF (iget(225)>0) THEN
573  IF (lvls(l,iget(225))>0) THEN
574  IF(isf_surface_physics==3)THEN
575  if(grib=='grib2') then
576  cfld=cfld+1
577  fld_info(cfld)%ifld=iavblfld(iget(225))
578  fld_info(cfld)%lvl=lvlsxml(l,iget(225))
579 !$omp parallel do private(i,j,ii,jj)
580  do j=1,jend-jsta+1
581  jj = jsta+j-1
582  do i=1,iend-ista+1
583  ii = ista+i-1
584  datapd(i,j,cfld) = sh2o(ii,jj,l)
585  enddo
586  enddo
587  endif
588  ELSE
589  dtop = 0.
590  DO ls=1,l-1
591  dtop = dtop + sldpth(ls)
592  ENDDO
593  dbot = dtop + sldpth(l)
594  if(grib=='grib2') then
595  cfld=cfld+1
596  fld_info(cfld)%ifld=iavblfld(iget(225))
597  fld_info(cfld)%lvl=lvlsxml(l,iget(225))
598 !$omp parallel do private(i,j,ii,jj)
599  do j=1,jend-jsta+1
600  jj = jsta+j-1
601  do i=1,iend-ista+1
602  ii = ista+i-1
603  datapd(i,j,cfld) = sh2o(ii,jj,l)
604  enddo
605  enddo
606  endif
607  ENDIF
608  ENDIF
609  ENDIF
610  ENDDO ! END OF NSOIL LOOP
611 ! -----------------
612 !
613 ! BOTTOM SOIL TEMPERATURE.
614  IF (iget(115)>0.or.iget(571)>0) THEN
615  if(iget(115)>0) then
616  if(grib=='grib2') then
617  cfld=cfld+1
618  fld_info(cfld)%ifld=iavblfld(iget(115))
619 !$omp parallel do private(i,j,ii,jj)
620  do j=1,jend-jsta+1
621  jj = jsta+j-1
622  do i=1,iend-ista+1
623  ii = ista+i-1
624  datapd(i,j,cfld) = tg(ii,jj)
625  enddo
626  enddo
627  endif
628  endif
629  if(iget(571)>0.and.grib=='grib2') then
630  cfld=cfld+1
631  fld_info(cfld)%ifld=iavblfld(iget(571))
632 !$omp parallel do private(i,j,ii,jj)
633  do j=1,jend-jsta+1
634  jj = jsta+j-1
635  do i=1,iend-ista+1
636  ii = ista+i-1
637  datapd(i,j,cfld) = tg(ii,jj)
638  enddo
639  enddo
640  endif
641  ENDIF
642 !
643 ! SOIL MOISTURE AVAILABILITY
644  IF (iget(171)>0) THEN
645 !!$omp parallel do private(i,j)
646  DO j=jsta,jend
647  DO i=ista,iend
648  IF(smstav(i,j) /= spval)THEN
649  grid1(i,j) = smstav(i,j)*100.
650  ELSE
651  grid1(i,j) = 0.
652  ENDIF
653  ENDDO
654  ENDDO
655  if(grib=='grib2') then
656  cfld=cfld+1
657  fld_info(cfld)%ifld=iavblfld(iget(171))
658 !$omp parallel do private(i,j,ii,jj)
659  do j=1,jend-jsta+1
660  jj = jsta+j-1
661  do i=1,iend-ista+1
662  ii = ista+i-1
663  datapd(i,j,cfld) = grid1(ii,jj)
664  enddo
665  enddo
666  endif
667  ENDIF
668 !
669 ! TOTAL SOIL MOISTURE
670  IF (iget(036)>0) THEN
671 !$omp parallel do private(i,j)
672  DO j=jsta,jend
673  DO i=ista,iend
674  IF(smstot(i,j)/=spval) THEN
675  IF(sm(i,j) > small .AND. sice(i,j) < small) THEN
676  grid1(i,j) = 1000.0 ! TEMPORY FIX TO MAKE SURE SMSTOT=1 FOR WATER
677  ELSE
678  grid1(i,j) = smstot(i,j)
679  END IF
680  ELSE
681  grid1(i,j) = 1000.0
682  ENDIF
683  ENDDO
684  ENDDO
685  if(grib=='grib2') then
686  cfld=cfld+1
687  fld_info(cfld)%ifld=iavblfld(iget(036))
688 !$omp parallel do private(i,j,ii,jj)
689  do j=1,jend-jsta+1
690  jj = jsta+j-1
691  do i=1,iend-ista+1
692  ii = ista+i-1
693  datapd(i,j,cfld) = grid1(ii,jj)
694  enddo
695  enddo
696  endif
697  ENDIF
698 !
699 ! PLANT CANOPY SURFACE WATER.
700  IF ( iget(118)>0 ) THEN
701  IF(modelname == 'RAPR') THEN
702 !$omp parallel do private(i,j)
703  DO j=jsta,jend
704  DO i=ista,iend
705  IF(cmc(i,j) /= spval) then
706  grid1(i,j) = cmc(i,j)
707  else
708  grid1(i,j) = spval
709  endif
710  ENDDO
711  ENDDO
712  else
713 !$omp parallel do private(i,j)
714  DO j=jsta,jend
715  DO i=ista,iend
716  IF(cmc(i,j) /= spval) then
717  grid1(i,j) = cmc(i,j)*1000.
718  else
719  grid1(i,j) = spval
720  endif
721  ENDDO
722  ENDDO
723  endif
724  if(grib=='grib2') then
725  cfld=cfld+1
726  fld_info(cfld)%ifld=iavblfld(iget(118))
727 !$omp parallel do private(i,j,ii,jj)
728  do j=1,jend-jsta+1
729  jj = jsta+j-1
730  do i=1,iend-ista+1
731  ii = ista+i-1
732  datapd(i,j,cfld) = grid1(ii,jj)
733  enddo
734  enddo
735  endif
736  ENDIF
737 !
738 ! SNOW WATER EQUIVALENT.
739  IF ( iget(119)>0 ) THEN
740 ! GRID1 = SPVAL
741  if(grib=='grib2') then
742  cfld=cfld+1
743  fld_info(cfld)%ifld=iavblfld(iget(119))
744 !$omp parallel do private(i,j,ii,jj)
745  do j=1,jend-jsta+1
746  jj = jsta+j-1
747  do i=1,iend-ista+1
748  ii = ista+i-1
749  datapd(i,j,cfld) = sno(ii,jj)
750  enddo
751  enddo
752  endiF
753  ENDIF
754 !
755 ! Time averaged percent SNOW COVER (for AQ)
756  IF ( iget(500)>0 ) THEN
757 ! GRID1=SPVAL
758 !$omp parallel do private(i,j)
759  DO j=jsta,jend
760  DO i=ista,iend
761 ! GRID1(I,J) = 100.*SNOAVG(I,J)
762  grid1(i,j) = snoavg(i,j)
763  if (snoavg(i,j) /= spval) grid1(i,j) = 100.*snoavg(i,j)
764  ENDDO
765  ENDDO
766  CALL bound(grid1,d00,h100)
767  id(1:25) = 0
768  itsrfc = nint(tsrfc)
769  IF(itsrfc /= 0) then
770  ifincr = mod(ifhr,itsrfc)
771  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
772  ELSE
773  ifincr = 0
774  endif
775  id(19) = ifhr
776  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
777  id(20) = 3
778  IF (ifincr==0) THEN
779  id(18) = ifhr-itsrfc
780  ELSE
781  id(18) = ifhr-ifincr
782  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
783  ENDIF
784  IF (id(18)<0) id(18) = 0
785  if(grib=='grib2') then
786  cfld=cfld+1
787  fld_info(cfld)%ifld=iavblfld(iget(500))
788  if(itsrfc>0) then
789  fld_info(cfld)%ntrange=1
790  else
791  fld_info(cfld)%ntrange=0
792  endif
793  fld_info(cfld)%tinvstat=ifhr-id(18)
794  ! fld_info(cfld)%ntrange=IFHR-ID(18)
795  ! fld_info(cfld)%tinvstat=1
796 !$omp parallel do private(i,j,ii,jj)
797  do j=1,jend-jsta+1
798  jj = jsta+j-1
799  do i=1,iend-ista+1
800  ii = ista+i-1
801  datapd(i,j,cfld) = grid1(ii,jj)
802  enddo
803  enddo
804  endif
805  ENDIF
806 
807 ! Time averaged surface pressure (for AQ)
808  IF ( iget(501)>0 ) THEN
809 ! GRID1 = SPVAL
810  id(1:25) = 0
811  id(19) = ifhr
812  IF (ifhr==0) THEN
813  id(18) = 0
814  ELSE
815  id(18) = ifhr - 1
816  ENDIF
817  id(20) = 3
818  if(grib=='grib2') then
819  cfld=cfld+1
820  fld_info(cfld)%ifld=iavblfld(iget(501))
821  fld_info(cfld)%ntrange=ifhr-id(18)
822  fld_info(cfld)%tinvstat=1
823 !$omp parallel do private(i,j,ii,jj)
824  do j=1,jend-jsta+1
825  jj = jsta+j-1
826  do i=1,iend-ista+1
827  ii = ista+i-1
828  datapd(i,j,cfld) = psfcavg(ii,jj)
829  enddo
830  enddo
831  endif
832  ENDIF
833 
834 ! Time averaged 10 m temperature (for AQ)
835  IF ( iget(502)>0 ) THEN
836 ! GRID1 = SPVAL
837  id(1:25) = 0
838  id(19) = ifhr
839  IF (ifhr==0) THEN
840  id(18) = 0
841  ELSE
842  id(18) = ifhr - 1
843  ENDIF
844  id(20) = 3
845  isvalue = 10
846  id(10) = mod(isvalue/256,256)
847  id(11) = mod(isvalue,256)
848  if(grib=='grib2') then
849  cfld=cfld+1
850  fld_info(cfld)%ifld=iavblfld(iget(502))
851  fld_info(cfld)%ntrange=ifhr-id(18)
852  fld_info(cfld)%tinvstat=1
853 !$omp parallel do private(i,j,ii,jj)
854  do j=1,jend-jsta+1
855  jj = jsta+j-1
856  do i=1,iend-ista+1
857  ii = ista+i-1
858  datapd(i,j,cfld) = t10avg(ii,jj)
859  enddo
860  enddo
861  endif
862  ENDIF
863 !
864 ! ACM GRID SCALE SNOW AND ICE
865  IF ( iget(244)>0 ) THEN
866 !$omp parallel do private(i,j)
867  DO j=jsta,jend
868  DO i=ista,iend
869  grid1(i,j) = snonc(i,j)
870  ENDDO
871  ENDDO
872  id(1:25) = 0
873  itprec = nint(tprec)
874 !mp
875  if (itprec /= 0) then
876  ifincr = mod(ifhr,itprec)
877  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
878  else
879  ifincr = 0
880  endif
881 !mp
882  id(18) = 0
883  id(19) = ifhr
884  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
885  id(20) = 4
886  IF (ifincr==0) THEN
887  id(18) = ifhr-itprec
888  ELSE
889  id(18) = ifhr-ifincr
890  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
891  ENDIF
892  IF (id(18)<0) id(18) = 0
893 
894  if(grib=='grib2') then
895  cfld=cfld+1
896  fld_info(cfld)%ifld=iavblfld(iget(244))
897  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
898  endif
899  ENDIF
900 !
901 ! PERCENT SNOW COVER.
902  IF ( iget(120)>0 ) THEN
903  grid1=spval
904  DO j=jsta,jend
905  DO i=ista,iend
906 ! GRID1(I,J)=PCTSNO(I,J)
907  IF ( sno(i,j) /= spval ) THEN
908  sneqv = sno(i,j)
909  iveg = ivgtyp(i,j)
910  IF(iveg==0)iveg=7
911  CALL snfrac(sneqv,iveg,sncovr)
912  grid1(i,j) = sncovr*100.
913  ENDIF
914  ENDDO
915  ENDDO
916  CALL bound(grid1,d00,h100)
917  if(grib=='grib2') then
918  cfld=cfld+1
919  fld_info(cfld)%ifld=iavblfld(iget(120))
920 !$omp parallel do private(i,j,ii,jj)
921  do j=1,jend-jsta+1
922  jj = jsta+j-1
923  do i=1,iend-ista+1
924  ii = ista+i-1
925  datapd(i,j,cfld) = grid1(ii,jj)
926  enddo
927  enddo
928  endif
929  ENDIF
930 ! ADD SNOW DEPTH
931  IF ( iget(224)>0 ) THEN
932  ii = (ista+iend)/2
933  jj = (jsta+jend)/2
934 ! GRID1=SPVAL
935 !$omp parallel do private(i,j)
936  DO j=jsta,jend
937  DO i=ista,iend
938  grid1(i,j) = spval
939  IF(si(i,j) /= spval) grid1(i,j) = si(i,j)*0.001 ! SI comes out of WRF in mm
940  ENDDO
941  ENDDO
942 ! print*,'sample snow depth in GRIBIT= ',si(ii,jj)
943  if(grib=='grib2') then
944  cfld=cfld+1
945  fld_info(cfld)%ifld=iavblfld(iget(224))
946 !$omp parallel do private(i,j,ii,jj)
947  do j=1,jend-jsta+1
948  jj = jsta+j-1
949  do i=1,iend-ista+1
950  ii = ista+i-1
951  datapd(i,j,cfld) = grid1(ii,jj)
952  enddo
953  enddo
954  endif
955  ENDIF
956 ! ADD POTENTIAL EVAPORATION
957  IF ( iget(242)>0 ) THEN
958  if(grib=='grib2') then
959  cfld=cfld+1
960  fld_info(cfld)%ifld=iavblfld(iget(242))
961 !$omp parallel do private(i,j,ii,jj)
962  do j=1,jend-jsta+1
963  jj = jsta+j-1
964  do i=1,iend-ista+1
965  ii = ista+i-1
966  datapd(i,j,cfld) = potevp(ii,jj)
967  enddo
968  enddo
969  endif
970  ENDIF
971 ! ADD ICE THICKNESS
972  IF ( iget(349)>0 ) THEN
973  if(grib=='grib2') then
974  cfld=cfld+1
975  fld_info(cfld)%ifld=iavblfld(iget(349))
976 !$omp parallel do private(i,j,ii,jj)
977  do j=1,jend-jsta+1
978  jj = jsta+j-1
979  do i=1,iend-ista+1
980  ii = ista+i-1
981  datapd(i,j,cfld) = dzice(ii,jj)
982  enddo
983  enddo
984  endif
985  ENDIF
986 
987 ! ADD EC,EDIR,ETRANS,ESNOW,SMCDRY,SMCMAX
988 ! ONLY OUTPUT NEW LSM FIELDS FOR NMM AND ARW BECAUSE RSM USES OLD SOIL TYPES
989  IF (modelname == 'NCAR'.OR. modelname == 'NMM' &
990  .OR. modelname == 'FV3R' .OR. modelname == 'RAPR') THEN
991 ! write(0,*)'in surf,isltyp=',maxval(isltyp(1:im,jsta:jend)), &
992 ! minval(isltyp(1:im,jsta:jend)),'qwbs=',maxval(qwbs(1:im,jsta:jend)), &
993 ! minval(qwbs(1:im,jsta:jend)),'potsvp=',maxval(potevp(1:im,jsta:jend)), &
994 ! minval(potevp(1:im,jsta:jend)),'sno=',maxval(sno(1:im,jsta:jend)), &
995 ! minval(sno(1:im,jsta:jend)),'vegfrc=',maxval(vegfrc(1:im,jsta:jend)), &
996 ! minval(vegfrc(1:im,jsta:jend)), 'sh2o=',maxval(sh2o(1:im,jsta:jend,1)), &
997 ! minval(sh2o(1:im,jsta:jend,1)),'cmc=',maxval(cmc(1:im,jsta:jend)), &
998 ! minval(cmc(1:im,jsta:jend))
999  IF ( iget(228)>0 .OR. iget(229)>0 &
1000  .OR.iget(230)>0 .OR. iget(231)>0 &
1001  .OR.iget(232)>0 .OR. iget(233)>0) THEN
1002 
1003  allocate(smcdry(ista:iend,jsta:jend), &
1004  smcmax(ista:iend,jsta:jend))
1005  DO j=jsta,jend
1006  DO i=ista,iend
1007 ! ----------------------------------------------------------------------
1008 ! IF(QWBS(I,J)>0.001)print*,'NONZERO QWBS',i,j,QWBS(I,J)
1009 ! IF(abs(SM(I,J)-0.)<1.0E-5)THEN
1010 ! WRF ARW has no POTEVP field. So has to block out RAPR
1011  IF( (modelname/='RAPR') .AND. (abs(sm(i,j)-0.) < 1.0e-5) .AND. &
1012  & (abs(sice(i,j)-0.) < 1.0e-5) ) THEN
1013  CALL etcalc(qwbs(i,j),potevp(i,j),sno(i,j),vegfrc(i,j) &
1014  & , isltyp(i,j),sh2o(i,j,1:1),cmc(i,j) &
1015  & , ecan(i,j),edir(i,j),etrans(i,j),esnow(i,j) &
1016  & , smcdry(i,j),smcmax(i,j) )
1017  ELSE
1018  ecan(i,j) = 0.
1019  edir(i,j) = 0.
1020  etrans(i,j) = 0.
1021  esnow(i,j) = 0.
1022  smcdry(i,j) = 0.
1023  smcmax(i,j) = 0.
1024  ENDIF
1025  ENDDO
1026  ENDDO
1027 
1028  IF ( iget(228)>0 )THEN
1029  if(grib=='grib2') then
1030  cfld=cfld+1
1031  fld_info(cfld)%ifld=iavblfld(iget(228))
1032 !$omp parallel do private(i,j,ii,jj)
1033  do j=1,jend-jsta+1
1034  jj = jsta+j-1
1035  do i=1,iend-ista+1
1036  ii = ista+i-1
1037  datapd(i,j,cfld) = ecan(ii,jj)
1038  enddo
1039  enddo
1040  endiF
1041  ENDIF
1042 
1043  IF ( iget(229)>0 )THEN
1044  if(grib=='grib2') then
1045  cfld=cfld+1
1046  fld_info(cfld)%ifld=iavblfld(iget(229))
1047 !$omp parallel do private(i,j,ii,jj)
1048  do j=1,jend-jsta+1
1049  jj = jsta+j-1
1050  do i=1,iend-ista+1
1051  ii = ista+i-1
1052  datapd(i,j,cfld) = edir(ii,jj)
1053  enddo
1054  enddo
1055  endif
1056  ENDIF
1057 
1058  IF ( iget(230)>0 )THEN
1059  if(grib=='grib2') then
1060  cfld=cfld+1
1061  fld_info(cfld)%ifld=iavblfld(iget(230))
1062  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = etrans(ista:iend,jsta:jend)
1063  endif
1064  ENDIF
1065 
1066  IF ( iget(231)>0 )THEN
1067  if(grib=='grib2') then
1068  cfld=cfld+1
1069  fld_info(cfld)%ifld=iavblfld(iget(231))
1070  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = esnow(ista:iend,jsta:jend)
1071  endif
1072  ENDIF
1073 
1074  IF ( iget(232)>0 )THEN
1075  if(grib=='grib2') then
1076  cfld=cfld+1
1077  fld_info(cfld)%ifld=iavblfld(iget(232))
1078 !$omp parallel do private(i,j,ii,jj)
1079  do j=1,jend-jsta+1
1080  jj = jsta+j-1
1081  do i=1,iend-ista+1
1082  ii = ista+i-1
1083  datapd(i,j,cfld) = smcdry(ii,jj)
1084  enddo
1085  enddo
1086  endif
1087  ENDIF
1088 
1089  IF ( iget(233)>0 )THEN
1090  if(grib=='grib2') then
1091  cfld=cfld+1
1092  fld_info(cfld)%ifld=iavblfld(iget(233))
1093 !$omp parallel do private(i,j,ii,jj)
1094  do j=1,jend-jsta+1
1095  jj = jsta+j-1
1096  do i=1,iend-ista+1
1097  ii = ista+i-1
1098  datapd(i,j,cfld) = smcmax(ii,jj)
1099  enddo
1100  enddo
1101  endif
1102  ENDIF
1103 
1104  ENDIF
1105 ! if (allocated(ecan)) deallocate(ecan)
1106 ! if (allocated(edir)) deallocate(edir)
1107 ! if (allocated(etrans)) deallocate(etrans)
1108 ! if (allocated(esnow)) deallocate(esnow)
1109  if (allocated(smcdry)) deallocate(smcdry)
1110  if (allocated(smcmax)) deallocate(smcmax)
1111 
1112  END IF ! endif for ncar and nmm options
1113 
1114  IF ( iget(512)>0 )THEN
1115  if(grib=='grib2') then
1116  cfld=cfld+1
1117  fld_info(cfld)%ifld=iavblfld(iget(512))
1118 !$omp parallel do private(i,j,ii,jj)
1119  do j=1,jend-jsta+1
1120  jj = jsta+j-1
1121  do i=1,iend-ista+1
1122  ii = ista+i-1
1123  datapd(i,j,cfld) = acond(ii,jj)
1124  enddo
1125  enddo
1126  endiF
1127  ENDIF
1128 
1129  IF ( iget(513)>0 )THEN
1130  id(1:25) = 0
1131  itsrfc = nint(tsrfc)
1132  IF(itsrfc /= 0) then
1133  ifincr = mod(ifhr,itsrfc)
1134  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
1135  ELSE
1136  ifincr = 0
1137  endif
1138  id(19) = ifhr
1139  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1140  id(20) = 3
1141  IF (ifincr==0) THEN
1142  id(18) = ifhr-itsrfc
1143  ELSE
1144  id(18) = ifhr-ifincr
1145  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1146  ENDIF
1147  IF (id(18)<0) id(18) = 0
1148  if(grib=='grib2') then
1149  cfld=cfld+1
1150  fld_info(cfld)%ifld=iavblfld(iget(513))
1151  if(itsrfc>0) then
1152  fld_info(cfld)%ntrange=1
1153  else
1154  fld_info(cfld)%ntrange=0
1155  endif
1156  fld_info(cfld)%tinvstat=ifhr-id(18)
1157 !$omp parallel do private(i,j,ii,jj)
1158  do j=1,jend-jsta+1
1159  jj = jsta+j-1
1160  do i=1,iend-ista+1
1161  ii = ista+i-1
1162  datapd(i,j,cfld) = avgecan(ii,jj)
1163  enddo
1164  enddo
1165  endiF
1166  ENDIF
1167 
1168  IF ( iget(514)>0 )THEN
1169  id(1:25) = 0
1170  itsrfc = nint(tsrfc)
1171  IF(itsrfc /= 0) then
1172  ifincr = mod(ifhr,itsrfc)
1173  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
1174  ELSE
1175  ifincr = 0
1176  endif
1177  id(19) = ifhr
1178  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1179  id(20) = 3
1180  IF (ifincr==0) THEN
1181  id(18) = ifhr-itsrfc
1182  ELSE
1183  id(18) = ifhr-ifincr
1184  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1185  ENDIF
1186  IF (id(18)<0) id(18) = 0
1187  if(grib=='grib2') then
1188  cfld=cfld+1
1189  fld_info(cfld)%ifld=iavblfld(iget(514))
1190  if(itsrfc>0) then
1191  fld_info(cfld)%ntrange=1
1192  else
1193  fld_info(cfld)%ntrange=0
1194  endif
1195  fld_info(cfld)%tinvstat=ifhr-id(18)
1196 !$omp parallel do private(i,j,ii,jj)
1197  do j=1,jend-jsta+1
1198  jj = jsta+j-1
1199  do i=1,iend-ista+1
1200  ii = ista+i-1
1201  datapd(i,j,cfld) = avgedir(ii,jj)
1202  enddo
1203  enddo
1204  endif
1205  ENDIF
1206 
1207  IF ( iget(515)>0 )THEN
1208  id(1:25) = 0
1209  itsrfc = nint(tsrfc)
1210  IF(itsrfc /= 0) then
1211  ifincr = mod(ifhr,itsrfc)
1212  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
1213  ELSE
1214  ifincr = 0
1215  endif
1216  id(19) = ifhr
1217  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1218  id(20) = 3
1219  IF (ifincr==0) THEN
1220  id(18) = ifhr-itsrfc
1221  ELSE
1222  id(18) = ifhr-ifincr
1223  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1224  ENDIF
1225  IF (id(18)<0) id(18) = 0
1226  if(grib=='grib2') then
1227  cfld=cfld+1
1228  fld_info(cfld)%ifld=iavblfld(iget(515))
1229  if(itsrfc>0) then
1230  fld_info(cfld)%ntrange=1
1231  else
1232  fld_info(cfld)%ntrange=0
1233  endif
1234  fld_info(cfld)%tinvstat=ifhr-id(18)
1235  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = avgetrans(ista:iend,jsta:jend)
1236  endif
1237  ENDIF
1238 
1239  IF ( iget(516)>0 )THEN
1240  id(1:25) = 0
1241  itsrfc = nint(tsrfc)
1242  IF(itsrfc /= 0) then
1243  ifincr = mod(ifhr,itsrfc)
1244  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
1245  ELSE
1246  ifincr = 0
1247  endif
1248  id(19) = ifhr
1249  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1250  id(20) = 3
1251  IF (ifincr==0) THEN
1252  id(18) = ifhr-itsrfc
1253  ELSE
1254  id(18) = ifhr-ifincr
1255  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1256  ENDIF
1257  IF (id(18)<0) id(18) = 0
1258  if(grib=='grib2') then
1259  cfld=cfld+1
1260  fld_info(cfld)%ifld=iavblfld(iget(516))
1261  if(itsrfc>0) then
1262  fld_info(cfld)%ntrange=1
1263  else
1264  fld_info(cfld)%ntrange=0
1265  endif
1266  fld_info(cfld)%tinvstat=ifhr-id(18)
1267  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = avgesnow(ista:iend,jsta:jend)
1268  endif
1269  ENDIF
1270 
1271  IF ( iget(996)>0 )THEN
1272  if(grib=='grib2') then
1273  cfld=cfld+1
1274  fld_info(cfld)%ifld=iavblfld(iget(996))
1275 !$omp parallel do private(i,j,ii,jj)
1276  do j=1,jend-jsta+1
1277  jj = jsta+j-1
1278  do i=1,iend-ista+1
1279  ii = ista+i-1
1280  datapd(i,j,cfld) = landfrac(ii,jj)
1281  enddo
1282  enddo
1283  endif
1284  ENDIF
1285 
1286  IF ( iget(997)>0 )THEN
1287  if(grib=='grib2') then
1288  cfld=cfld+1
1289  fld_info(cfld)%ifld=iavblfld(iget(997))
1290 !$omp parallel do private(i,j,ii,jj)
1291  do j=1,jend-jsta+1
1292  jj = jsta+j-1
1293  do i=1,iend-ista+1
1294  ii = ista+i-1
1295  datapd(i,j,cfld) = pahi(ii,jj)
1296  enddo
1297  enddo
1298  endif
1299  ENDIF
1300 
1301  IF ( iget(998)>0 )THEN
1302  if(grib=='grib2') then
1303  cfld=cfld+1
1304  fld_info(cfld)%ifld=iavblfld(iget(998))
1305 !$omp parallel do private(i,j,ii,jj)
1306  do j=1,jend-jsta+1
1307  jj = jsta+j-1
1308  do i=1,iend-ista+1
1309  ii = ista+i-1
1310  datapd(i,j,cfld) = twa(ii,jj)
1311  enddo
1312  enddo
1313  endif
1314  ENDIF
1315 
1316  IF ( iget(999)>0 )THEN
1317 !$omp parallel do private(i,j)
1318  DO j=jsta,jend
1319  DO i=ista,iend
1320  grid1(i,j) = tecan(i,j)
1321  ENDDO
1322  ENDDO
1323  id(1:25) = 0
1324  itprec = nint(tprec)
1325  if (itprec /= 0) then
1326  ifincr = mod(ifhr,itprec)
1327  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
1328  else
1329  ifincr = 0
1330  endif
1331  id(18) = 0
1332  id(19) = ifhr
1333  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1334  id(20) = 4
1335  IF (ifincr==0) THEN
1336  id(18) = ifhr-itprec
1337  ELSE
1338  id(18) = ifhr-ifincr
1339  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1340  ENDIF
1341  IF (id(18)<0) id(18) = 0
1342  if(grib=='grib2') then
1343  cfld=cfld+1
1344  fld_info(cfld)%ifld=iavblfld(iget(999))
1345  fld_info(cfld)%ntrange=1
1346  fld_info(cfld)%tinvstat=ifhr-id(18)
1347 !$omp parallel do private(i,j,ii,jj)
1348  do j=1,jend-jsta+1
1349  jj = jsta+j-1
1350  do i=1,iend-ista+1
1351  ii = ista+i-1
1352  datapd(i,j,cfld) = grid1(ii,jj)
1353  enddo
1354  enddo
1355  endif
1356  ENDIF
1357 
1358  IF ( iget(1000)>0 )THEN
1359 !$omp parallel do private(i,j)
1360  DO j=jsta,jend
1361  DO i=ista,iend
1362  grid1(i,j) = tetran(i,j)
1363  ENDDO
1364  ENDDO
1365  id(1:25) = 0
1366  itprec = nint(tprec)
1367  if (itprec /= 0) then
1368  ifincr = mod(ifhr,itprec)
1369  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
1370  else
1371  ifincr = 0
1372  endif
1373  id(18) = 0
1374  id(19) = ifhr
1375  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1376  id(20) = 4
1377  IF (ifincr==0) THEN
1378  id(18) = ifhr-itprec
1379  ELSE
1380  id(18) = ifhr-ifincr
1381  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1382  ENDIF
1383  IF (id(18)<0) id(18) = 0
1384  if(grib=='grib2') then
1385  cfld=cfld+1
1386  fld_info(cfld)%ifld=iavblfld(iget(1000))
1387  fld_info(cfld)%ntrange=1
1388  fld_info(cfld)%tinvstat=ifhr-id(18)
1389 !$omp parallel do private(i,j,ii,jj)
1390  do j=1,jend-jsta+1
1391  jj = jsta+j-1
1392  do i=1,iend-ista+1
1393  ii = ista+i-1
1394  datapd(i,j,cfld) = grid1(ii,jj)
1395  enddo
1396  enddo
1397  endif
1398  ENDIF
1399 !
1400  IF ( iget(1001)>0 )THEN
1401 !$omp parallel do private(i,j)
1402  DO j=jsta,jend
1403  DO i=ista,iend
1404  grid1(i,j) = tedir(i,j)
1405  ENDDO
1406  ENDDO
1407  id(1:25) = 0
1408  itprec = nint(tprec)
1409  if (itprec /= 0) then
1410  ifincr = mod(ifhr,itprec)
1411  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
1412  else
1413  ifincr = 0
1414  endif
1415  id(18) = 0
1416  id(19) = ifhr
1417  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1418  id(20) = 4
1419  IF (ifincr==0) THEN
1420  id(18) = ifhr-itprec
1421  ELSE
1422  id(18) = ifhr-ifincr
1423  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1424  ENDIF
1425  IF (id(18)<0) id(18) = 0
1426  if(grib=='grib2') then
1427  cfld=cfld+1
1428  fld_info(cfld)%ifld=iavblfld(iget(1001))
1429  fld_info(cfld)%ntrange=1
1430  fld_info(cfld)%tinvstat=ifhr-id(18)
1431 !$omp parallel do private(i,j,ii,jj)
1432  do j=1,jend-jsta+1
1433  jj = jsta+j-1
1434  do i=1,iend-ista+1
1435  ii = ista+i-1
1436  datapd(i,j,cfld) = grid1(ii,jj)
1437  enddo
1438  enddo
1439  endif
1440  ENDIF
1441 !
1442 
1443  IF (iget(1002)>0) THEN
1444  IF(asrfc>0.)THEN
1445  rrnum=1./asrfc
1446  ELSE
1447  rrnum=0.
1448  ENDIF
1449  DO j=jsta,jend
1450  DO i=ista,iend
1451  IF(paha(i,j)/=spval)THEN
1452  grid1(i,j)=-1.*paha(i,j)*rrnum !change the sign to conform with Grib
1453  ELSE
1454  grid1(i,j)=paha(i,j)
1455  END IF
1456  ENDDO
1457  ENDDO
1458  id(1:25) = 0
1459  itsrfc = nint(tsrfc)
1460  IF(itsrfc /= 0) then
1461  ifincr = mod(ifhr,itsrfc)
1462  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
1463  ELSE
1464  ifincr = 0
1465  endif
1466  id(19) = ifhr
1467  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1468  id(20) = 3
1469  IF (ifincr==0) THEN
1470  id(18) = ifhr-itsrfc
1471  ELSE
1472  id(18) = ifhr-ifincr
1473  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1474  ENDIF
1475  IF (id(18)<0) id(18) = 0
1476  if(grib=='grib2') then
1477  cfld=cfld+1
1478  fld_info(cfld)%ifld=iavblfld(iget(1002))
1479  if(itsrfc>0) then
1480  fld_info(cfld)%ntrange=1
1481  else
1482  fld_info(cfld)%ntrange=0
1483  endif
1484  fld_info(cfld)%tinvstat=ifhr-id(18)
1485  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1486  endif
1487  ENDIF
1488 !
1489 !
1490 !
1491 !*** BLOCK 2. SHELTER (2M) LEVEL FIELDS.
1492 !
1493 ! COMPUTE/POST SHELTER LEVEL FIELDS.
1494 !
1495  IF ( (iget(106)>0).OR.(iget(112)>0).OR. &
1496  (iget(113)>0).OR.(iget(114)>0).OR. &
1497  (iget(138)>0).OR.(iget(414)>0).OR. &
1498  (iget(546)>0).OR.(iget(547)>0).OR. &
1499  (iget(548)>0).OR.(iget(739)>0).OR. &
1500  (iget(771)>0)) THEN
1501 
1502  if (.not. allocated(psfc)) allocate(psfc(ista:iend,jsta:jend))
1503 !
1504 !HC COMPUTE SHELTER PRESSURE BECAUSE IT WAS NOT OUTPUT FROM WRF
1505  IF(modelname == 'NCAR' .OR. modelname=='RSM'.OR. modelname=='RAPR')THEN
1506  DO j=jsta,jend
1507  DO i=ista,iend
1508  tlow = t(i,j,nint(lmh(i,j)))
1509  psfc(i,j) = pint(i,j,nint(lmh(i,j))+1) !May not have been set above
1510  pshltr(i,j) = psfc(i,j)*exp(-0.068283/tlow)
1511  ENDDO
1512  ENDDO
1513  ENDIF
1514 !
1515 ! print *,'in, surfc,pshltr=',maxval(PSHLTR(1:im,jsta:jend)), &
1516 ! minval(PSHLTR(1:im,jsta:jend)),PSHLTR(1:3,jsta),'capa=',capa, &
1517 ! 'tshlter=',tshltr(1:3,jsta:jsta+2),'psfc=',psfc(1:3,jsta:jsta+2), &
1518 ! 'th10=',th10(1:3,jsta:jsta+2),'thz0=',thz0(1:3,jsta:jsta+2)
1519 !
1520 ! SHELTER LEVEL TEMPERATURE
1521  IF (iget(106)>0) THEN
1522  grid1=spval
1523  DO j=jsta,jend
1524  DO i=ista,iend
1525 ! GRID1(I,J)=TSHLTR(I,J)
1526 !HC CONVERT FROM THETA TO T
1527  if(tshltr(i,j)/=spval)grid1(i,j)=tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
1528  IF(grid1(i,j)<200)print*,'ABNORMAL 2MT ',i,j, &
1529  tshltr(i,j),pshltr(i,j)
1530 ! TSHLTR(I,J)=GRID1(I,J)
1531  ENDDO
1532  ENDDO
1533 ! print *,'2m tmp=',maxval(TSHLTR(ista:iend,jsta:jend)), &
1534 ! minval(TSHLTR(ista:iend,jsta:jend)),TSHLTR(1:3,jsta),'grd=',grid1(1:3,jsta)
1535  if(grib=='grib2') then
1536  cfld=cfld+1
1537  fld_info(cfld)%ifld=iavblfld(iget(106))
1538  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
1539  endif
1540  ENDIF
1541 !
1542 ! SHELTER LEVEL POT TEMP
1543  IF (iget(546)>0) THEN
1544 ! GRID1=spval
1545 ! DO J=JSTA,JEND
1546 ! DO I=ISTA,IEND
1547 ! GRID1(I,J)=TSHLTR(I,J)
1548 ! ENDDO
1549 ! ENDDO
1550  if(grib=='grib2') then
1551  cfld=cfld+1
1552  fld_info(cfld)%ifld=iavblfld(iget(546))
1553  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = tshltr(ista:iend,jsta:jend)
1554  endif
1555  ENDIF
1556 !
1557 ! SHELTER LEVEL SPECIFIC HUMIDITY.
1558  IF (iget(112)>0) THEN
1559  DO j=jsta,jend
1560  DO i=ista,iend
1561  grid1(i,j) = qshltr(i,j)
1562  ENDDO
1563  ENDDO
1564  CALL bound(grid1,h1m12,h99999)
1565  if(grib=='grib2') then
1566  cfld=cfld+1
1567  fld_info(cfld)%ifld=iavblfld(iget(112))
1568  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
1569  endif
1570  ENDIF
1571 ! GRID1
1572 ! SHELTER MIXING RATIO.
1573  IF (iget(414)>0) THEN
1574  DO j=jsta,jend
1575  DO i=ista,iend
1576  grid1(i,j) = mrshltr(i,j)
1577  ENDDO
1578  ENDDO
1579  if(grib=='grib2') then
1580  cfld=cfld+1
1581  fld_info(cfld)%ifld=iavblfld(iget(414))
1582  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1583  endif
1584  ENDIF
1585 !
1586 ! SHELTER LEVEL DEWPOINT, DEWPOINT DEPRESSION AND SFC EQUIV POT TEMP.
1587  allocate(p1d(ista:iend,jsta:jend), t1d(ista:iend,jsta:jend))
1588  IF ((iget(113)>0) .OR.(iget(547)>0).OR.(iget(548)>0)) THEN
1589 
1590  DO j=jsta,jend
1591  DO i=ista,iend
1592 
1593 !tgs The next 4 lines are GSD algorithm for Dew Point computation
1594 !tgs Results are very close to dew point computed in DEWPOINT subroutine
1595  qv = max(1.e-5,(qshltr(i,j)/(1.-qshltr(i,j))))
1596  e = pshltr(i,j)/100.*qv/(0.62197+qv)
1597  dwpt = (243.5*log(e)-440.8)/(19.48-log(e))+273.15
1598 
1599 ! if(i==335.and.j==295)print*,'Debug: RUC-type DEWPT,i,j' &
1600 ! if(i==ii.and.j==jj)print*,'Debug: RUC-type DEWPT,i,j'
1601 ! , DWPT,i,j,qv,pshltr(i,j),qshltr(i,j)
1602 
1603 ! EGRID1(I,J) = DWPT
1604 
1605  IF(qshltr(i,j)<spval.and.pshltr(i,j)<spval)THEN
1606  evp(i,j) = pshltr(i,j)*qshltr(i,j)/(eps+oneps*qshltr(i,j))
1607  evp(i,j) = evp(i,j)*d001
1608  ELSE
1609  evp(i,j) = spval
1610  ENDIF
1611  ENDDO
1612  ENDDO
1613  CALL dewpoint(evp,egrid1(ista:iend,jsta:jend))
1614 ! print *,' MAX DEWPOINT',maxval(egrid1)
1615 ! DEWPOINT
1616  IF (iget(113)>0) THEN
1617  grid1=spval
1618  if(modelname=='RAPR')THEN
1619  DO j=jsta,jend
1620  DO i=ista,iend
1621 ! DEWPOINT can't be higher than T2
1622  t2=tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
1623  if(qshltr(i,j)/=spval)grid1(i,j)=min(egrid1(i,j),t2)
1624  ENDDO
1625  ENDDO
1626  else
1627  DO j=jsta,jend
1628  DO i=ista,iend
1629  if(qshltr(i,j)/=spval) grid1(i,j) = egrid1(i,j)
1630  ENDDO
1631  ENDDO
1632  endif
1633  if(grib=='grib2') then
1634  cfld=cfld+1
1635  fld_info(cfld)%ifld=iavblfld(iget(113))
1636  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1637  endif
1638  ENDIF
1639 
1640 
1641 !-------------------------------------------------------------------------
1642 ! DEWPOINT at level 1 ------ p1d and t1d are undefined !! -- Moorthi
1643  IF (iget(771)>0) THEN
1644  DO j=jsta,jend
1645  DO i=ista,iend
1646  evp(i,j)=p1d(i,j)*qvl1(i,j)/(eps+oneps*qvl1(i,j))
1647  evp(i,j)=evp(i,j)*d001
1648  ENDDO
1649  ENDDO
1650  CALL dewpoint(evp,egrid1(ista:iend,jsta:jend))
1651 ! print *,' MAX DEWPOINT at level 1',maxval(egrid1)
1652  grid1=spval
1653  DO j=jsta,jend
1654  DO i=ista,iend
1655 !tgs 30 dec 2013 - 1st leel dewpoint can't be higher than 1-st level temperature
1656  if(qvl1(i,j)/=spval)grid1(i,j) = min(egrid1(i,j),t1d(i,j))
1657  ENDDO
1658  ENDDO
1659  if(grib=='grib2') then
1660  cfld=cfld+1
1661  fld_info(cfld)%ifld=iavblfld(iget(771))
1662  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1663  endif
1664  ENDIF
1665 !-------------------------------------------------------------------------
1666 
1667 !
1668  IF ((iget(547)>0).OR.(iget(548)>0)) THEN
1669  grid1=spval
1670  grid2=spval
1671  DO j=jsta,jend
1672  DO i=ista,iend
1673  if(tshltr(i,j)/=spval.and.pshltr(i,j)/=spval.and.qshltr(i,j)/=spval) then
1674 ! DEWPOINT DEPRESSION in GRID1
1675  grid1(i,j)=max(0.,tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa-egrid1(i,j))
1676 
1677 ! SURFACE EQIV POT TEMP in GRID2
1678  ape=(h10e5/pshltr(i,j))**capa
1679  grid2(i,j)=tshltr(i,j)*exp(elocp*qshltr(i,j)*ape/tshltr(i,j))
1680  endif
1681  ENDDO
1682  ENDDO
1683 ! print *,' MAX/MIN --> DEWPOINT DEPRESSION',maxval(grid1(1:im,jsta:jend)),&
1684 ! minval(grid1(1:im,jsta:jend))
1685 ! print *,' MAX/MIN --> SFC EQUIV POT TEMP',maxval(grid2(1:im,jsta:jend)),&
1686 ! minval(grid2(1:im,jsta:jend))
1687 
1688  IF (iget(547)>0) THEN
1689  if(grib=='grib2') then
1690  cfld=cfld+1
1691  fld_info(cfld)%ifld=iavblfld(iget(547))
1692  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1693  endif
1694 
1695  ENDIF
1696  IF (iget(548)>0) THEN
1697  if(grib=='grib2') then
1698  cfld=cfld+1
1699  fld_info(cfld)%ifld=iavblfld(iget(548))
1700  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid2(ista:iend,jsta:jend)
1701  endif
1702  ENDIF
1703  ENDIF
1704 
1705 
1706  ENDIF
1707 !
1708 ! SHELTER LEVEL RELATIVE HUMIDITY AND APPARENT TEMPERATURE
1709  IF (iget(114) > 0 .OR. iget(808) > 0) THEN
1710  allocate(q1d(ista:iend,jsta:jend))
1711 !$omp parallel do private(i,j,llmh)
1712  DO j=jsta,jend
1713  DO i=ista,iend
1714  IF(modelname=='RAPR')THEN
1715  llmh = nint(lmh(i,j))
1716 ! P1D(I,J)=PINT(I,J,LLMH+1)
1717  p1d(i,j) = pmid(i,j,llmh)
1718  t1d(i,j) = t(i,j,llmh)
1719  ELSE
1720  p1d(i,j) = pshltr(i,j)
1721  t1d(i,j) = tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
1722  ENDIF
1723  q1d(i,j) = qshltr(i,j)
1724  ENDDO
1725  ENDDO
1726 
1727  CALL calrh(p1d,t1d,q1d,egrid1(ista:iend,jsta:jend))
1728 
1729  if (allocated(q1d)) deallocate(q1d)
1730 !$omp parallel do private(i,j)
1731  DO j=jsta,jend
1732  DO i=ista,iend
1733  if(qshltr(i,j) /= spval)then
1734  grid1(i,j) = egrid1(i,j)*100.
1735  else
1736  grid1(i,j) = spval
1737  end if
1738  ENDDO
1739  ENDDO
1740  CALL bound(grid1,h1,h100)
1741  IF (iget(114) > 0) THEN
1742  if(grib == 'grib2') then
1743  cfld = cfld+1
1744  fld_info(cfld)%ifld = iavblfld(iget(114))
1745 !$omp parallel do private(i,j,ii,jj)
1746  do j=1,jend-jsta+1
1747  jj = jsta+j-1
1748  do i=1,iend-ista+1
1749  ii = ista+i-1
1750  datapd(i,j,cfld) = grid1(ii,jj)
1751  enddo
1752  enddo
1753  endif
1754  ENDIF
1755 
1756  IF(iget(808)>0)THEN
1757  grid2=spval
1758 !$omp parallel do private(i,j,dum1,dum2,dum3,dum216,dum1s,dum3s)
1759  DO j=jsta,jend
1760  DO i=ista,iend
1761  if(t1d(i,j)/=spval.and.u10h(i,j)/=spval.and.v10h(i,j)<spval) then
1762  dum1 = (t1d(i,j)-tfrz)*1.8+32.
1763  dum2 = sqrt(u10h(i,j)**2.0+v10h(i,j)**2.0)/0.44704
1764  dum3 = egrid1(i,j) * 100.0
1765 ! if(abs(gdlon(i,j)-120.)<1. .and. abs(gdlat(i,j))<1.) &
1766 ! print*,'Debug AT: INPUT', T1D(i,j),dum1,dum2,dum3
1767  IF(dum1 <= 50.) THEN
1768  dum216 = dum2**0.16
1769  grid2(i,j) = 35.74 + 0.6215*dum1 &
1770  - 35.75*dum216 + 0.4275*dum1*dum216
1771  grid2(i,j) =(grid2(i,j)-32.)/1.8+tfrz
1772  ELSE IF(dum1 > 80.) THEN
1773  dum1s = dum1*dum1
1774  dum3s = dum3*dum3
1775  grid2(i,j) = -42.379 + 2.04901523*dum1 &
1776  + 10.14333127*dum3 &
1777  - 0.22475541*dum1*dum3 &
1778  - 0.00683783*dum1s &
1779  - 0.05481717*dum3s &
1780  + 0.00122874*dum1s*dum3 &
1781  + 0.00085282*dum1*dum3s &
1782  - 0.00000199*dum1s*dum3s
1783  grid2(i,j) = (grid2(i,j)-32.)/1.8 + tfrz
1784  ELSE
1785  grid2(i,j) = t1d(i,j)
1786  END IF
1787 ! if(abs(gdlon(i,j)-120.)<1. .and. abs(gdlat(i,j))<1.) &
1788 ! print*,'Debug AT: OUTPUT',Grid2(i,j)
1789  endif
1790  ENDDO
1791  ENDDO
1792 
1793  if(grib == 'grib2') then
1794  cfld = cfld+1
1795  fld_info(cfld)%ifld = iavblfld(iget(808))
1796 !$omp parallel do private(i,j,ii,jj)
1797  do j=1,jend-jsta+1
1798  jj = jsta+j-1
1799  do i=1,iend-ista+1
1800  ii = ista+i-1
1801  datapd(i,j,cfld) = grid2(ii,jj)
1802  enddo
1803  enddo
1804  endif
1805 
1806  ENDIF !for 808
1807 
1808  ENDIF ! ENDIF for shleter RH or apparent T
1809 
1810  if (allocated(p1d)) deallocate (p1d)
1811  if (allocated(t1d)) deallocate (t1d)
1812 !
1813 ! SHELTER LEVEL PRESSURE.
1814  IF (iget(138)>0) THEN
1815 ! DO J=JSTA,JEND
1816 ! DO I=ISTA,IEND
1817 ! GRID1(I,J)=PSHLTR(I,J)
1818 ! ENDDO
1819 ! ENDDO
1820  if(grib=='grib2') then
1821  cfld=cfld+1
1822  fld_info(cfld)%ifld=iavblfld(iget(138))
1823 !$omp parallel do private(i,j,ii,jj)
1824  do j=1,jend-jsta+1
1825  jj = jsta+j-1
1826  do i=1,iend-ista+1
1827  ii = ista+i-1
1828  datapd(i,j,cfld) = pshltr(ii,jj)
1829  enddo
1830  enddo
1831  endif
1832  ENDIF
1833 !
1834  ENDIF
1835 !
1836 ! SHELTER LEVEL MAX TEMPERATURE.
1837  IF (iget(345)>0) THEN
1838 ! DO J=JSTA,JEND
1839 ! DO I=ISTA,IEND
1840 ! GRID1(I,J)=MAXTSHLTR(I,J)
1841 ! ENDDO
1842 ! ENDDO
1843 !mp
1844  tmaxmin = max(tmaxmin,1.)
1845 !mp
1846  itmaxmin = int(tmaxmin)
1847  IF(itmaxmin /= 0) then
1848  ifincr = mod(ifhr,itmaxmin)
1849  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
1850  ELSE
1851  ifincr = 0
1852  endif
1853  id(19) = ifhr
1854  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1855  id(20) = 2
1856  IF (ifincr==0) THEN
1857  id(18) = ifhr-itmaxmin
1858  ELSE
1859  id(18) = ifhr-ifincr
1860  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1861  ENDIF
1862  IF (id(18)<0) id(18) = 0
1863  if(grib=='grib2') then
1864  cfld=cfld+1
1865  fld_info(cfld)%ifld=iavblfld(iget(345))
1866  if(itmaxmin==0) then
1867  fld_info(cfld)%ntrange=0
1868  else
1869  fld_info(cfld)%ntrange=1
1870  endif
1871  fld_info(cfld)%tinvstat=ifhr-id(18)
1872  if(ifhr==0) fld_info(cfld)%tinvstat=0
1873 !$omp parallel do private(i,j,ii,jj)
1874  do j=1,jend-jsta+1
1875  jj = jsta+j-1
1876  do i=1,iend-ista+1
1877  ii = ista+i-1
1878  datapd(i,j,cfld) = maxtshltr(ii,jj)
1879  enddo
1880  enddo
1881  endif
1882  ENDIF
1883 !
1884 ! SHELTER LEVEL MIN TEMPERATURE.
1885  IF (iget(346)>0) THEN
1886 !!$omp parallel do private(i,j)
1887 ! DO J=JSTA,JEND
1888 ! DO I=ISTA,IEND
1889 ! GRID1(I,J) = MINTSHLTR(I,J)
1890 ! ENDDO
1891 ! ENDDO
1892  id(1:25) = 0
1893  itmaxmin = int(tmaxmin)
1894  IF(itmaxmin /= 0) then
1895  ifincr = mod(ifhr,itmaxmin)
1896  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
1897  ELSE
1898  ifincr = 0
1899  endif
1900  id(19) = ifhr
1901  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1902  id(20) = 2
1903  IF (ifincr==0) THEN
1904  id(18) = ifhr-itmaxmin
1905  ELSE
1906  id(18) = ifhr-ifincr
1907  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1908  ENDIF
1909  IF (id(18)<0) id(18) = 0
1910  if(grib=='grib2') then
1911  cfld=cfld+1
1912  fld_info(cfld)%ifld=iavblfld(iget(346))
1913  if(itmaxmin==0) then
1914  fld_info(cfld)%ntrange=0
1915  else
1916  fld_info(cfld)%ntrange=1
1917  endif
1918  fld_info(cfld)%tinvstat=ifhr-id(18)
1919  if(ifhr==0) fld_info(cfld)%tinvstat=0
1920 !$omp parallel do private(i,j,ii,jj)
1921  do j=1,jend-jsta+1
1922  jj = jsta+j-1
1923  do i=1,iend-ista+1
1924  ii = ista+i-1
1925  datapd(i,j,cfld) = mintshltr(ii,jj)
1926  enddo
1927  enddo
1928  endif
1929  ENDIF
1930 !
1931 ! SHELTER LEVEL MAX RH.
1932  IF (iget(347)>0) THEN
1933  grid1=spval
1934  DO j=jsta,jend
1935  DO i=ista,iend
1936  if(maxrhshltr(i,j)/=spval) grid1(i,j)=maxrhshltr(i,j)*100.
1937  ENDDO
1938  ENDDO
1939  id(1:25) = 0
1940  id(02)=129
1941  itmaxmin = int(tmaxmin)
1942  IF(itmaxmin /= 0) then
1943  ifincr = mod(ifhr,itmaxmin)
1944  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
1945  ELSE
1946  ifincr = 0
1947  endif
1948  id(19) = ifhr
1949  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1950  id(20) = 2
1951  IF (ifincr==0) THEN
1952  id(18) = ifhr-itmaxmin
1953  ELSE
1954  id(18) = ifhr-ifincr
1955  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1956  ENDIF
1957  IF (id(18)<0) id(18) = 0
1958  if(grib=='grib2') then
1959  cfld=cfld+1
1960  fld_info(cfld)%ifld=iavblfld(iget(347))
1961  if(itmaxmin==0) then
1962  fld_info(cfld)%ntrange=0
1963  else
1964 !Meng 03/2019
1965 ! fld_info(cfld)%ntrange=(IFHR-ID(18))/ITMAXMIN
1966  fld_info(cfld)%ntrange=1
1967  endif
1968 ! fld_info(cfld)%tinvstat=ITMAXMIN
1969  fld_info(cfld)%tinvstat=ifhr-id(18)
1970  if(ifhr==0) fld_info(cfld)%tinvstat=0
1971 ! print*,'id(18),tinvstat,IFHR,ITMAXMIN in rhmax= ',ID(18),fld_info(cfld)%tinvstat, &
1972 ! IFHR, ITMAXMIN
1973 !$omp parallel do private(i,j,ii,jj)
1974  do j=1,jend-jsta+1
1975  jj = jsta+j-1
1976  do i=1,iend-ista+1
1977  ii = ista+i-1
1978  datapd(i,j,cfld) = grid1(ii,jj)
1979  enddo
1980  enddo
1981  endif
1982  ENDIF
1983 !
1984 ! SHELTER LEVEL MIN RH.
1985  IF (iget(348)>0) THEN
1986  grid1=spval
1987  DO j=jsta,jend
1988  DO i=ista,iend
1989  if(minrhshltr(i,j)/=spval) grid1(i,j)=minrhshltr(i,j)*100.
1990  ENDDO
1991  ENDDO
1992  id(1:25) = 0
1993  id(02)=129
1994  itmaxmin = int(tmaxmin)
1995  IF(itmaxmin /= 0) then
1996  ifincr = mod(ifhr,itmaxmin)
1997  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
1998  ELSE
1999  ifincr = 0
2000  endif
2001  id(19) = ifhr
2002  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2003  id(20) = 2
2004  IF (ifincr==0) THEN
2005  id(18) = ifhr-itmaxmin
2006  ELSE
2007  id(18) = ifhr-ifincr
2008  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2009  ENDIF
2010  IF (id(18)<0) id(18) = 0
2011  if(grib=='grib2') then
2012  cfld=cfld+1
2013  fld_info(cfld)%ifld=iavblfld(iget(348))
2014  if(itmaxmin==0) then
2015  fld_info(cfld)%ntrange=0
2016  else
2017 !Meng 03/2019
2018 ! fld_info(cfld)%ntrange=(IFHR-ID(18))/ITMAXMIN
2019  fld_info(cfld)%ntrange=1
2020  endif
2021 ! fld_info(cfld)%tinvstat=ITMAXMIN
2022  fld_info(cfld)%tinvstat=ifhr-id(18)
2023  if(ifhr==0) fld_info(cfld)%tinvstat=0
2024 !$omp parallel do private(i,j,ii,jj)
2025  do j=1,jend-jsta+1
2026  jj = jsta+j-1
2027  do i=1,iend-ista+1
2028  ii = ista+i-1
2029  datapd(i,j,cfld) = grid1(ii,jj)
2030  enddo
2031  enddo
2032  endif
2033  ENDIF
2034 
2035 !
2036 ! SHELTER LEVEL MAX SPFH
2037  IF (iget(510)>0) THEN
2038  id(1:25) = 0
2039  itmaxmin = int(tmaxmin)
2040  IF(itmaxmin /= 0) then
2041  ifincr = mod(ifhr,itmaxmin)
2042  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
2043  ELSE
2044  ifincr = 0
2045  endif
2046  id(19) = ifhr
2047  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2048  id(20) = 2
2049  IF (ifincr==0) THEN
2050  id(18) = ifhr-itmaxmin
2051  ELSE
2052  id(18) = ifhr-ifincr
2053  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2054  ENDIF
2055  IF (id(18)<0) id(18) = 0
2056  if(grib=='grib2') then
2057  cfld=cfld+1
2058  fld_info(cfld)%ifld=iavblfld(iget(510))
2059  if(itmaxmin==0) then
2060  fld_info(cfld)%ntrange=0
2061  else
2062  fld_info(cfld)%ntrange=1
2063  endif
2064  fld_info(cfld)%tinvstat=ifhr-id(18)
2065 !$omp parallel do private(i,j,ii,jj)
2066  do j=1,jend-jsta+1
2067  jj = jsta+j-1
2068  do i=1,iend-ista+1
2069  ii = ista+i-1
2070  datapd(i,j,cfld) = maxqshltr(ii,jj)
2071  enddo
2072  enddo
2073  endif
2074  ENDIF
2075 !
2076 ! SHELTER LEVEL MIN SPFH
2077  IF (iget(511)>0) THEN
2078  id(1:25) = 0
2079  itmaxmin = int(tmaxmin)
2080  IF(itmaxmin /= 0) then
2081  ifincr = mod(ifhr,itmaxmin)
2082  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itmaxmin*60)
2083  ELSE
2084  ifincr = 0
2085  endif
2086  id(19) = ifhr
2087  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2088  id(20) = 2
2089  IF (ifincr==0) THEN
2090  id(18) = ifhr-itmaxmin
2091  ELSE
2092  id(18) = ifhr-ifincr
2093  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2094  ENDIF
2095  IF (id(18)<0) id(18) = 0
2096  if(grib=='grib2') then
2097  cfld=cfld+1
2098  fld_info(cfld)%ifld=iavblfld(iget(511))
2099  if(itmaxmin==0) then
2100  fld_info(cfld)%ntrange=0
2101  else
2102  fld_info(cfld)%ntrange=1
2103  endif
2104  fld_info(cfld)%tinvstat=ifhr-id(18)
2105 !$omp parallel do private(i,j,ii,jj)
2106  do j=1,jend-jsta+1
2107  jj = jsta+j-1
2108  do i=1,iend-ista+1
2109  ii = ista+i-1
2110  datapd(i,j,cfld) = minqshltr(ii,jj)
2111  enddo
2112  enddo
2113  endif
2114  ENDIF
2115 !
2116 ! E. James - 12 Sep 2018: SMOKE from WRF-CHEM on lowest model level
2117 !
2118  IF (iget(739)>0) THEN
2119  grid1=spval
2120  DO j=jsta,jend
2121  DO i=ista,iend
2122  if(t(i,j,lm)/=spval.and.pmid(i,j,lm)/=spval.and.smoke(i,j,lm,1)/=spval)&
2123  grid1(i,j) = (1./rd)*(pmid(i,j,lm)/t(i,j,lm))*smoke(i,j,lm,1)
2124  ENDDO
2125  ENDDO
2126  if(grib=='grib2') then
2127  cfld=cfld+1
2128  fld_info(cfld)%ifld=iavblfld(iget(739))
2129  datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
2130  endif
2131  ENDIF
2132 !
2133 ! BLOCK 3. ANEMOMETER LEVEL (10M) WINDS, THETA, AND Q.
2134 !
2135  IF ( (iget(064)>0).OR.(iget(065)>0).OR. &
2136  (iget(506)>0).OR.(iget(507)>0) ) THEN
2137 !
2138 ! ANEMOMETER LEVEL U WIND AND/OR V WIND.
2139  IF ((iget(064)>0).OR.(iget(065)>0)) THEN
2140 !$omp parallel do private(i,j)
2141  DO j=jsta,jend
2142  DO i=ista,iend
2143  grid1(i,j) = u10(i,j)
2144  grid2(i,j) = v10(i,j)
2145  ENDDO
2146  ENDDO
2147  if(grib=='grib2') then
2148  cfld=cfld+1
2149  fld_info(cfld)%ifld=iavblfld(iget(064))
2150 !$omp parallel do private(i,j,ii,jj)
2151  do j=1,jend-jsta+1
2152  jj = jsta+j-1
2153  do i=1,iend-ista+1
2154  ii = ista+i-1
2155  datapd(i,j,cfld) = grid1(ii,jj)
2156  enddo
2157  enddo
2158  cfld=cfld+1
2159  fld_info(cfld)%ifld=iavblfld(iget(065))
2160 !$omp parallel do private(i,j,ii,jj)
2161  do j=1,jend-jsta+1
2162  jj = jsta+j-1
2163  do i=1,iend-ista+1
2164  ii = ista+i-1
2165  datapd(i,j,cfld) = grid2(ii,jj)
2166  enddo
2167  enddo
2168  endif
2169  ENDIF
2170 ! GSD - Time-averaged wind speed (forecast time labels will all be in minutes)
2171  IF (iget(730)>0) THEN
2172  ifincr = 5
2173  DO j=jsta,jend
2174  DO i=ista,iend
2175  grid1(i,j)=spduv10mean(i,j)
2176  ENDDO
2177  ENDDO
2178  if(grib=='grib2') then
2179 ! print*,'Outputting time-averaged winds'
2180  cfld=cfld+1
2181  fld_info(cfld)%ifld=iavblfld(iget(730))
2182  if(fld_info(cfld)%ntrange==0) then
2183  if (ifhr==0 .and. ifmin==0) then
2184  fld_info(cfld)%tinvstat=0
2185  else
2186  fld_info(cfld)%tinvstat=ifincr
2187  endif
2188  fld_info(cfld)%ntrange=1
2189  end if
2190  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2191  endif
2192  ENDIF
2193 !---
2194 ! GSD - Time-averaged U wind speed (forecast time labels will all be in minutes)
2195  IF (iget(731)>0) THEN
2196  ifincr = 5
2197  DO j=jsta,jend
2198  DO i=ista,iend
2199  grid1(i,j)=u10mean(i,j)
2200  ENDDO
2201  ENDDO
2202  if(grib=='grib2') then
2203  cfld=cfld+1
2204  fld_info(cfld)%ifld=iavblfld(iget(731))
2205  if(fld_info(cfld)%ntrange==0) then
2206  if (ifhr==0 .and. ifmin==0) then
2207  fld_info(cfld)%tinvstat=0
2208  else
2209  fld_info(cfld)%tinvstat=ifincr
2210  endif
2211  fld_info(cfld)%ntrange=1
2212  end if
2213  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2214  endif
2215  ENDIF
2216 ! GSD - Time-averaged V wind speed (forecast time labels will all be in minutes)
2217  IF (iget(732)>0) THEN
2218  ifincr = 5
2219  DO j=jsta,jend
2220  DO i=ista,iend
2221  grid1(i,j)=v10mean(i,j)
2222  ENDDO
2223  ENDDO
2224  if(grib=='grib2') then
2225  cfld=cfld+1
2226  fld_info(cfld)%ifld=iavblfld(iget(732))
2227  if(fld_info(cfld)%ntrange==0) then
2228  if (ifhr==0 .and. ifmin==0) then
2229  fld_info(cfld)%tinvstat=0
2230  else
2231  fld_info(cfld)%tinvstat=ifincr
2232  endif
2233  fld_info(cfld)%ntrange=1
2234  end if
2235  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2236  endif
2237  ENDIF
2238 ! Time-averaged SWDOWN (forecast time labels will all be in minutes)
2239  IF (iget(733)>0 )THEN
2240  ifincr = 15
2241  DO j=jsta,jend
2242  DO i=ista,iend
2243  grid1(i,j) = swradmean(i,j)
2244  ENDDO
2245  ENDDO
2246  if(grib=='grib2') then
2247  cfld=cfld+1
2248  fld_info(cfld)%ifld=iavblfld(iget(733))
2249  if(fld_info(cfld)%ntrange==0) then
2250  if (ifhr==0 .and. ifmin==0) then
2251  fld_info(cfld)%tinvstat=0
2252  else
2253  fld_info(cfld)%tinvstat=ifincr
2254  endif
2255  fld_info(cfld)%ntrange=1
2256  end if
2257  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2258  endif
2259  ENDIF
2260 ! Time-averaged SWNORM (forecast time labels will all be in minutes)
2261  IF (iget(734)>0 )THEN
2262  ifincr = 15
2263  DO j=jsta,jend
2264  DO i=ista,iend
2265  grid1(i,j) = swnormmean(i,j)
2266  ENDDO
2267  ENDDO
2268  if(grib=='grib2') then
2269  cfld=cfld+1
2270  fld_info(cfld)%ifld=iavblfld(iget(734))
2271  if(fld_info(cfld)%ntrange==0) then
2272  if (ifhr==0 .and. ifmin==0) then
2273  fld_info(cfld)%tinvstat=0
2274  else
2275  fld_info(cfld)%tinvstat=ifincr
2276  endif
2277  fld_info(cfld)%ntrange=1
2278  endif
2279  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2280  endif
2281  ENDIF
2282 !
2283  IF ((iget(506)>0).OR.(iget(507)>0)) THEN
2284  id(02)=129
2285  id(20) = 2
2286  id(19) = ifhr
2287  IF (ifhr==0) THEN
2288  id(18) = 0
2289  ELSE
2290  id(18) = ifhr - 1
2291  ENDIF
2292 !$omp parallel do private(i,j)
2293  DO j=jsta,jend
2294  DO i=ista,iend
2295  grid1(i,j) = u10max(i,j)
2296  grid2(i,j) = v10max(i,j)
2297  ENDDO
2298  ENDDO
2299  if(grib=='grib2') then
2300  cfld=cfld+1
2301  fld_info(cfld)%ifld=iavblfld(iget(506))
2302  fld_info(cfld)%ntrange=ifhr-id(18)
2303  fld_info(cfld)%tinvstat=1
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  cfld=cfld+1
2313  fld_info(cfld)%ifld=iavblfld(iget(507))
2314  fld_info(cfld)%ntrange=ifhr-id(18)
2315  fld_info(cfld)%tinvstat=1
2316 !$omp parallel do private(i,j,ii,jj)
2317  do j=1,jend-jsta+1
2318  jj = jsta+j-1
2319  do i=1,iend-ista+1
2320  ii = ista+i-1
2321  datapd(i,j,cfld) = grid2(ii,jj)
2322  enddo
2323  enddo
2324  endif
2325  ENDIF
2326 
2327  ENDIF
2328 !
2329 ! ANEMOMETER LEVEL (10 M) POTENTIAL TEMPERATURE.
2330 ! NOT A OUTPUT FROM WRF
2331  IF (iget(158)>0) THEN
2332 !$omp parallel do private(i,j)
2333  DO j=jsta,jend
2334  DO i=ista,iend
2335  grid1(i,j)=th10(i,j)
2336  ENDDO
2337  ENDDO
2338  if(grib=='grib2') then
2339  cfld=cfld+1
2340  fld_info(cfld)%ifld=iavblfld(iget(158))
2341 !$omp parallel do private(i,j,ii,jj)
2342  do j=1,jend-jsta+1
2343  jj = jsta+j-1
2344  do i=1,iend-ista+1
2345  ii = ista+i-1
2346  datapd(i,j,cfld) = grid1(ii,jj)
2347  enddo
2348  enddo
2349  endif
2350  ENDIF
2351 
2352 ! ANEMOMETER LEVEL (10 M) SENSIBLE TEMPERATURE.
2353 ! NOT A OUTPUT FROM WRF
2354  IF (iget(505)>0) THEN
2355 !$omp parallel do private(i,j)
2356  DO j=jsta,jend
2357  DO i=ista,iend
2358  grid1(i,j)=t10m(i,j)
2359  ENDDO
2360  ENDDO
2361  if(grib=='grib2') then
2362  cfld=cfld+1
2363  fld_info(cfld)%ifld=iavblfld(iget(505))
2364 !$omp parallel do private(i,j,ii,jj)
2365  do j=1,jend-jsta+1
2366  jj = jsta+j-1
2367  do i=1,iend-ista+1
2368  ii = ista+i-1
2369  datapd(i,j,cfld) = grid1(ii,jj)
2370  enddo
2371  enddo
2372  endif
2373  ENDIF
2374 !
2375 ! ANEMOMETER LEVEL (10 M) SPECIFIC HUMIDITY.
2376 !
2377  IF (iget(159)>0) THEN
2378 !$omp parallel do private(i,j)
2379  DO j=jsta,jend
2380  DO i=ista,iend
2381  grid1(i,j) = q10(i,j)
2382  ENDDO
2383  ENDDO
2384  if(grib=='grib2') then
2385  cfld=cfld+1
2386  fld_info(cfld)%ifld=iavblfld(iget(159))
2387 !$omp parallel do private(i,j,ii,jj)
2388  do j=1,jend-jsta+1
2389  jj = jsta+j-1
2390  do i=1,iend-ista+1
2391  ii = ista+i-1
2392  datapd(i,j,cfld) = grid1(ii,jj)
2393  enddo
2394  enddo
2395  endif
2396  ENDIF
2397 !
2398 ! SRD
2399 !
2400 ! ANEMOMETER LEVEL (10 M) MAX WIND SPEED.
2401 !
2402  IF (iget(422)>0) THEN
2403 !$omp parallel do private(i,j)
2404  DO j=jsta,jend
2405  DO i=ista,iend
2406  grid1(i,j) = wspd10max(i,j)
2407  ENDDO
2408  ENDDO
2409  if(grib=='grib2') then
2410  cfld=cfld+1
2411  fld_info(cfld)%ifld=iavblfld(iget(422))
2412  if (ifhr==0) then
2413  fld_info(cfld)%tinvstat=0
2414  else
2415  fld_info(cfld)%tinvstat=1
2416  endif
2417  fld_info(cfld)%ntrange=1
2418 !$omp parallel do private(i,j,ii,jj)
2419  do j=1,jend-jsta+1
2420  jj = jsta+j-1
2421  do i=1,iend-ista+1
2422  ii = ista+i-1
2423  datapd(i,j,cfld) = grid1(ii,jj)
2424  enddo
2425  enddo
2426  endif
2427  ENDIF
2428 
2429 ! ANEMOMETER LEVEL (10 M) MAX WIND SPEED U COMPONENT.
2430 !
2431  IF (iget(783)>0) THEN
2432 !$omp parallel do private(i,j)
2433  DO j=jsta,jend
2434  DO i=ista,iend
2435  grid1(i,j) = wspd10umax(i,j)
2436  ENDDO
2437  ENDDO
2438  if(grib=='grib2') then
2439  cfld=cfld+1
2440  fld_info(cfld)%ifld=iavblfld(iget(783))
2441  if (ifhr==0) then
2442  fld_info(cfld)%tinvstat=0
2443  else
2444  fld_info(cfld)%tinvstat=1
2445  endif
2446  fld_info(cfld)%ntrange=1
2447 !$omp parallel do private(i,j,ii,jj)
2448  do j=1,jend-jsta+1
2449  jj = jsta+j-1
2450  do i=1,iend-ista+1
2451  ii = ista+i-1
2452  datapd(i,j,cfld) = grid1(ii,jj)
2453  enddo
2454  enddo
2455  endif
2456  ENDIF
2457 
2458 ! ANEMOMETER LEVEL (10 M) MAX WIND SPEED V COMPONENT.
2459 !
2460  IF (iget(784)>0) THEN
2461 !$omp parallel do private(i,j)
2462  DO j=jsta,jend
2463  DO i=ista,iend
2464  grid1(i,j) = wspd10vmax(i,j)
2465  ENDDO
2466  ENDDO
2467  if(grib=='grib2') then
2468  cfld=cfld+1
2469  fld_info(cfld)%ifld=iavblfld(iget(784))
2470  if (ifhr==0) then
2471  fld_info(cfld)%tinvstat=0
2472  else
2473  fld_info(cfld)%tinvstat=1
2474  endif
2475  fld_info(cfld)%ntrange=1
2476 !$omp parallel do private(i,j,ii,jj)
2477  do j=1,jend-jsta+1
2478  jj = jsta+j-1
2479  do i=1,iend-ista+1
2480  ii = ista+i-1
2481  datapd(i,j,cfld) = grid1(i,jj)
2482  enddo
2483  enddo
2484  endif
2485  ENDIF
2486 
2487 !
2488 ! SRD
2489 !
2490 
2491 ! Ice Growth Rate
2492 !
2493  IF (iget(588)>0) THEN
2494 
2495  CALL calvessel(iceg(ista:iend,jsta:jend))
2496 
2497  DO j=jsta,jend
2498  DO i=ista,iend
2499  grid1(i,j) = iceg(i,j)
2500  ENDDO
2501  ENDDO
2502 
2503  if(grib=='grib2') then
2504  cfld=cfld+1
2505  fld_info(cfld)%ifld=iavblfld(iget(588))
2506  if (ifhr==0) then
2507  fld_info(cfld)%tinvstat=0
2508  else
2509  fld_info(cfld)%tinvstat=1
2510  endif
2511  fld_info(cfld)%ntrange=1
2512 
2513 !$omp parallel do private(i,j,ii,jj)
2514  do j=1,jend-jsta+1
2515  jj = jsta+j-1
2516  do i=1,iend-ista+1
2517  ii = ista+i-1
2518  datapd(i,j,cfld) = grid1(ii,jj)
2519  enddo
2520  enddo
2521  endif
2522 
2523  ENDIF
2524 
2525 !
2526 !*** BLOCK 4. PRECIPITATION RELATED FIELDS.
2527 !MEB 6/17/02 ASSUMING THAT ALL ACCUMULATED FIELDS NEVER EMPTY
2528 ! THEIR BUCKETS. THIS IS THE EASIEST WAY TO DEAL WITH
2529 ! ACCUMULATED FIELDS. SHORTER TIME ACCUMULATIONS CAN
2530 ! BE COMPUTED AFTER THE FACT IN A SEPARATE CODE ONCE
2531 ! THE POST HAS FINISHED. I HAVE LEFT IN THE OLD
2532 ! ETAPOST CODE FOR COMPUTING THE BEGINNING TIME OF
2533 ! THE ACCUMULATION PERIOD IF THIS IS CHANGED BACK
2534 ! TO A 12H OR 3H BUCKET. I AM NOT SURE WHAT
2535 ! TO DO WITH THE TIME AVERAGED FIELDS, SO
2536 ! LEAVING THAT UNCHANGED.
2537 !
2538 ! SNOW FRACTION FROM EXPLICIT CLOUD SCHEME. LABELLED AS
2539 ! 'PROB OF FROZEN PRECIP' IN GRIB,
2540 ! DIDN'T KNOW WHAT ELSE TO CALL IT
2541  IF (iget(172)>0) THEN
2542 !$omp parallel do private(i,j)
2543  DO j=jsta,jend
2544  DO i=ista,iend
2545  IF (prec(i,j) <= pthresh .OR. sr(i,j)==spval) THEN
2546  grid1(i,j) = -50.
2547  ELSE
2548  grid1(i,j) = sr(i,j)*100.
2549  ENDIF
2550  ENDDO
2551  ENDDO
2552  if(grib=='grib2') then
2553  cfld=cfld+1
2554  fld_info(cfld)%ifld=iavblfld(iget(172))
2555 !$omp parallel do private(i,j,ii,jj)
2556  do j=1,jend-jsta+1
2557  jj = jsta+j-1
2558  do i=1,iend-ista+1
2559  ii = ista+i-1
2560  datapd(i,j,cfld) = grid1(ii,jj)
2561  enddo
2562  enddo
2563  endif
2564  ENDIF
2565 !
2566 ! INSTANTANEOUS CONVECTIVE PRECIPITATION RATE.
2567 ! SUBSTITUTE WITH CUPPT IN WRF FOR NOW
2568  IF (iget(249)>0) THEN
2569  rdtphs=1000./dtq2 !--- 1000 kg/m**3, density of liquid water
2570 ! RDTPHS=1000./(TRDLW*3600.)
2571  grid1=spval
2572 !$omp parallel do private(i,j)
2573  DO j=jsta,jend
2574  DO i=ista,iend
2575  if(cprate(i,j)/=spval) grid1(i,j) = cprate(i,j)*rdtphs
2576 ! GRID1(I,J) = CUPPT(I,J)*RDTPHS
2577  ENDDO
2578  ENDDO
2579  if(grib=='grib2') then
2580  cfld=cfld+1
2581  fld_info(cfld)%ifld=iavblfld(iget(249))
2582 !$omp parallel do private(i,j,ii,jj)
2583  do j=1,jend-jsta+1
2584  jj = jsta+j-1
2585  do i=1,iend-ista+1
2586  ii = ista+i-1
2587  datapd(i,j,cfld) = grid1(ii,jj)
2588  enddo
2589  enddo
2590  endif
2591  ENDIF
2592 !
2593 ! INSTANTANEOUS PRECIPITATION RATE.
2594  IF (iget(167)>0) THEN
2595 !MEB need to get physics DT
2596  rdtphs=1./(dtq2)
2597 !MEB need to get physics DT
2598  grid1=spval
2599 !$omp parallel do private(i,j)
2600  DO j=jsta,jend
2601  DO i=ista,iend
2602  if(prec(i,j)/=spval) then
2603  IF(modelname /= 'RSM') THEN
2604  grid1(i,j) = prec(i,j)*rdtphs*1000.
2605  ELSE !Add by Binbin
2606  grid1(i,j) = prec(i,j)
2607  END IF
2608  endif
2609  ENDDO
2610  ENDDO
2611  if(grib=='grib2') then
2612  cfld=cfld+1
2613  fld_info(cfld)%ifld=iavblfld(iget(167))
2614 !$omp parallel do private(i,j,ii,jj)
2615  do j=1,jend-jsta+1
2616  jj = jsta+j-1
2617  do i=1,iend-ista+1
2618  ii = ista+i-1
2619  datapd(i,j,cfld) = grid1(ii,jj)
2620  enddo
2621  enddo
2622  endif
2623  ENDIF
2624 !
2625 ! MAXIMUM INSTANTANEOUS PRECIPITATION RATE.
2626  IF (iget(508)>0) THEN
2627 !-- PRATE_MAX in units of mm/h from NMMB history files
2628  grid1=spval
2629  DO j=jsta,jend
2630  DO i=ista,iend
2631  if(prate_max(i,j)/=spval) grid1(i,j)=prate_max(i,j)*sec2hr
2632  ENDDO
2633  ENDDO
2634  if(grib=='grib2') then
2635  cfld=cfld+1
2636  fld_info(cfld)%ifld=iavblfld(iget(508))
2637  fld_info(cfld)%lvl=lvlsxml(1,iget(508))
2638  fld_info(cfld)%tinvstat=1
2639  if (ifhr > 0) then
2640  fld_info(cfld)%ntrange=1
2641  else
2642  fld_info(cfld)%ntrange=0
2643  endif
2644 !$omp parallel do private(i,j,ii,jj)
2645  do j=1,jend-jsta+1
2646  jj = jsta+j-1
2647  do i=1,iend-ista+1
2648  ii = ista+i-1
2649  datapd(i,j,cfld) = grid1(ii,jj)
2650  enddo
2651  enddo
2652  endif
2653  ENDIF
2654 !
2655 ! MAXIMUM INSTANTANEOUS *FROZEN* PRECIPITATION RATE.
2656  IF (iget(509)>0) THEN
2657 !-- FPRATE_MAX in units of mm/h from NMMB history files
2658  grid1=spval
2659  DO j=jsta,jend
2660  DO i=ista,iend
2661  if(fprate_max(i,j)/=spval) grid1(i,j)=fprate_max(i,j)*sec2hr
2662  ENDDO
2663  ENDDO
2664  if(grib=='grib2') then
2665  cfld=cfld+1
2666  fld_info(cfld)%ifld=iavblfld(iget(509))
2667  fld_info(cfld)%lvl=lvlsxml(1,iget(509))
2668  fld_info(cfld)%tinvstat=1
2669  if (ifhr > 0) then
2670  fld_info(cfld)%ntrange=1
2671  else
2672  fld_info(cfld)%ntrange=0
2673  endif
2674 !$omp parallel do private(i,j,ii,jj)
2675  do j=1,jend-jsta+1
2676  jj = jsta+j-1
2677  do i=1,iend-ista+1
2678  ii = ista+i-1
2679  datapd(i,j,cfld) = grid1(ii,jj)
2680  enddo
2681  enddo
2682  endif
2683  ENDIF
2684 !
2685 ! TIME-AVERAGED CONVECTIVE PRECIPITATION RATE.
2686  IF (iget(272)>0) THEN
2687  rdtphs=1000./dtq2 !--- 1000 kg/m**3, density of liquid water
2688  id(1:25) = 0
2689  itprec = nint(tprec)
2690 !mp
2691  if (itprec /= 0) then
2692  ifincr = mod(ifhr,itprec)
2693  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
2694  else
2695  ifincr = 0
2696  endif
2697 !mp
2698  id(18) = 0
2699  id(19) = ifhr
2700  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2701  id(20) = 3
2702  IF (ifincr==0) THEN
2703  id(18) = ifhr-itprec
2704  ELSE
2705  id(18) = ifhr-ifincr
2706  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2707  ENDIF
2708  IF (id(18)<0) id(18) = 0
2709  grid1=spval
2710 !$omp parallel do private(i,j)
2711  DO j=jsta,jend
2712  DO i=ista,iend
2713  if(avgcprate(i,j)/=spval) grid1(i,j) = avgcprate(i,j)*rdtphs
2714  ENDDO
2715  ENDDO
2716 
2717 ! print *,'in surf,iget(272)=',iget(272),'RDTPHS=',RDTPHS, &
2718 ! 'AVGCPRATE=',maxval(AVGCPRATE(1:im,jsta:jend)),minval(AVGCPRATE(1:im,jsta:jend)), &
2719 ! 'grid1=',maxval(grid1(1:im,jsta:jend)),minval(grid1(1:im,jsta:jend))
2720  if(grib=='grib2') then
2721  cfld=cfld+1
2722  fld_info(cfld)%ifld=iavblfld(iget(272))
2723 
2724  if(itprec==0) then
2725  fld_info(cfld)%ntrange=0
2726  else
2727  fld_info(cfld)%ntrange=1
2728  endif
2729  fld_info(cfld)%tinvstat=ifhr-id(18)
2730 
2731 !$omp parallel do private(i,j,ii,jj)
2732  do j=1,jend-jsta+1
2733  jj = jsta+j-1
2734  do i=1,iend-ista+1
2735  ii = ista+i-1
2736  datapd(i,j,cfld) = grid1(ii,jj)
2737  enddo
2738  enddo
2739  endif
2740  ENDIF
2741 !
2742 ! TIME-AVERAGED PRECIPITATION RATE.
2743  IF (iget(271)>0) THEN
2744  rdtphs=1000./dtq2 !--- 1000 kg/m**3, density of liquid water
2745 ! RDTPHS=1000./(TRDLW*3600.)
2746  id(1:25) = 0
2747  itprec = nint(tprec)
2748 !mp
2749  if (itprec /= 0) then
2750  ifincr = mod(ifhr,itprec)
2751  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
2752  else
2753  ifincr = 0
2754  endif
2755 !mp
2756  id(18) = 0
2757  id(19) = ifhr
2758  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2759  id(20) = 3
2760  IF (ifincr==0) THEN
2761  id(18) = ifhr-itprec
2762  ELSE
2763  id(18) = ifhr-ifincr
2764  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2765  ENDIF
2766  IF (id(18)<0) id(18) = 0
2767  grid1=spval
2768 !$omp parallel do private(i,j)
2769  DO j=jsta,jend
2770  DO i=ista,iend
2771  if(avgprec(i,j)/=spval) grid1(i,j) = avgprec(i,j)*rdtphs
2772  ENDDO
2773  ENDDO
2774 
2775  if(grib=='grib2') then
2776  cfld=cfld+1
2777  fld_info(cfld)%ifld=iavblfld(iget(271))
2778 
2779  if(itprec==0) then
2780  fld_info(cfld)%ntrange=0
2781  else
2782  fld_info(cfld)%ntrange=1
2783  endif
2784  fld_info(cfld)%tinvstat=ifhr-id(18)
2785 
2786 !$omp parallel do private(i,j,ii,jj)
2787  do j=1,jend-jsta+1
2788  jj = jsta+j-1
2789  do i=1,iend-ista+1
2790  ii = ista+i-1
2791  datapd(i,j,cfld) = grid1(ii,jj)
2792  enddo
2793  enddo
2794  endif
2795  ENDIF
2796 !
2797 ! ACCUMULATED TOTAL PRECIPITATION.
2798  IF (iget(087)>0) THEN
2799  id(1:25) = 0
2800  itprec = nint(tprec)
2801 !mp
2802  if (itprec /= 0) then
2803  ifincr = mod(ifhr,itprec)
2804  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
2805  else
2806  ifincr = 0
2807  endif
2808 !mp
2809  id(18) = 0
2810  id(19) = ifhr
2811  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2812  id(20) = 4
2813  IF (ifincr==0) THEN
2814  id(18) = ifhr-itprec
2815  ELSE
2816  id(18) = ifhr-ifincr
2817  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2818  ENDIF
2819  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
2820 !$omp parallel do private(i,j)
2821  DO j=jsta,jend
2822  DO i=ista,iend
2823  IF(avgprec(i,j) < spval)THEN
2824  grid1(i,j) = avgprec(i,j)*float(id(19)-id(18))*3600.*1000./dtq2
2825  ELSE
2826  grid1(i,j) = spval
2827  END IF
2828  ENDDO
2829  ENDDO
2830 !! Chuang 3/29/2018: add continuous bucket
2831 ! DO J=JSTA,JEND
2832 ! DO I=ISTA,IEND
2833 ! IF(AVGPREC_CONT(I,J) < SPVAL)THEN
2834 ! GRID2(I,J) = AVGPREC_CONT(I,J)*FLOAT(IFHR)*3600.*1000./DTQ2
2835 ! ELSE
2836 ! GRID2(I,J) = SPVAL
2837 ! END IF
2838 ! ENDDO
2839 ! ENDDO
2840  ELSE
2841 !$omp parallel do private(i,j)
2842  DO j=jsta,jend
2843  DO i=ista,iend
2844  IF(acprec(i,j) < spval)THEN
2845  grid1(i,j) = acprec(i,j)*1000.
2846  ELSE
2847  grid1(i,j) = spval
2848  ENDIF
2849  ENDDO
2850  ENDDO
2851  END IF
2852 ! IF(IFMIN >= 1 .AND. ID(19) > 256)THEN
2853 ! IF(ITPREC==3)ID(17)=10
2854 ! IF(ITPREC==6)ID(17)=11
2855 ! IF(ITPREC==12)ID(17)=12
2856 ! END IF
2857  IF (id(18)<0) id(18) = 0
2858 ! write(6,*) 'call gribit...total precip'
2859  if(grib=='grib2') then
2860  cfld=cfld+1
2861  fld_info(cfld)%ifld=iavblfld(iget(087))
2862  fld_info(cfld)%ntrange=1
2863  fld_info(cfld)%tinvstat=ifhr-id(18)
2864 ! print*,'id(18),tinvstat in apcp= ',ID(18),fld_info(cfld)%tinvstat
2865 !$omp parallel do private(i,j,ii,jj)
2866  do j=1,jend-jsta+1
2867  jj = jsta+j-1
2868  do i=1,iend-ista+1
2869  ii = ista+i-1
2870  datapd(i,j,cfld) = grid1(ii,jj)
2871  enddo
2872  enddo
2873 !! add continuous bucket
2874 ! if(MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R') then
2875 ! cfld=cfld+1
2876 ! fld_info(cfld)%ifld=IAVBLFLD(IGET(087))
2877 ! fld_info(cfld)%ntrange=1
2878 ! fld_info(cfld)%tinvstat=IFHR
2879 ! print*,'tinvstat in cont bucket= ',fld_info(cfld)%tinvstat
2880 ! do j=1,jend-jsta+1
2881 ! jj = jsta+j-1
2882 ! do i=1,im
2883 ! datapd(i,j,cfld) = GRID2(i,jj)
2884 ! enddo
2885 ! enddo
2886 ! endif
2887  endif
2888  ENDIF
2889 
2890 !
2891 ! CONTINOUS ACCUMULATED TOTAL PRECIPITATION.
2892  IF (iget(417)>0) THEN
2893  id(1:25) = 0
2894  itprec = nint(tprec)
2895 !mp
2896  if (itprec /= 0) then
2897  ifincr = mod(ifhr,itprec)
2898  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
2899  else
2900  ifincr = 0
2901  endif
2902 !mp
2903  id(18) = 0
2904  id(19) = ifhr
2905  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2906  id(20) = 4
2907  IF (ifincr==0) THEN
2908  id(18) = ifhr-itprec
2909  ELSE
2910  id(18) = ifhr-ifincr
2911  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2912  ENDIF
2913  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
2914 ! Chuang 3/29/2018: add continuous bucket
2915 !$omp parallel do private(i,j)
2916  DO j=jsta,jend
2917  DO i=ista,iend
2918  IF(avgprec_cont(i,j) < spval)THEN
2919  grid2(i,j) = avgprec_cont(i,j)*float(ifhr)*3600.*1000./dtq2
2920  ELSE
2921  grid2(i,j) = spval
2922  END IF
2923  ENDDO
2924  ENDDO
2925  ENDIF
2926  IF (id(18)<0) id(18) = 0
2927  if(grib=='grib2') then
2928 ! add continuous bucket
2929  if(modelname == 'GFS' .OR. modelname == 'FV3R') then
2930  cfld=cfld+1
2931  fld_info(cfld)%ifld=iavblfld(iget(417))
2932  fld_info(cfld)%ntrange=1
2933  fld_info(cfld)%tinvstat=ifhr
2934 ! print*,'tinvstat in cont bucket= ',fld_info(cfld)%tinvstat
2935 !$omp parallel do private(i,j,ii,jj)
2936  do j=1,jend-jsta+1
2937  jj = jsta+j-1
2938  do i=1,iend-ista+1
2939  ii = ista+i-1
2940  datapd(i,j,cfld) = grid2(ii,jj)
2941  enddo
2942  enddo
2943  endif
2944  endif
2945  ENDIF
2946 !
2947 ! ACCUMULATED CONVECTIVE PRECIPITATION.
2948  IF (iget(033)>0) THEN
2949  id(1:25) = 0
2950  itprec = nint(tprec)
2951 !mp
2952  if (itprec /= 0) then
2953  ifincr = mod(ifhr,itprec)
2954  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
2955  else
2956  ifincr = 0
2957  endif
2958 !mp
2959  id(18) = 0
2960  id(19) = ifhr
2961  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2962  id(20) = 4
2963  IF (ifincr==0) THEN
2964  id(18) = ifhr-itprec
2965  ELSE
2966  id(18) = ifhr-ifincr
2967  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2968  ENDIF
2969  IF (id(18)<0) id(18) = 0
2970  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
2971 !$omp parallel do private(i,j)
2972  DO j=jsta,jend
2973  DO i=ista,iend
2974  IF(avgcprate(i,j) < spval)THEN
2975  grid1(i,j) = avgcprate(i,j)* &
2976  float(id(19)-id(18))*3600.*1000./dtq2
2977  ELSE
2978  grid1(i,j) = spval
2979  END IF
2980  ENDDO
2981  ENDDO
2982 !! Chuang 3/29/2018: add continuous bucket
2983 ! DO J=JSTA,JEND
2984 ! DO I=ISTA,IEND
2985 ! IF(AVGCPRATE_CONT(I,J) < SPVAL)THEN
2986 ! GRID2(I,J) = AVGCPRATE_CONT(I,J)*FLOAT(IFHR)*3600.*1000./DTQ2
2987 ! ELSE
2988 ! GRID2(I,J) = SPVAL
2989 ! END IF
2990 ! ENDDO
2991 ! ENDDO
2992  ELSE
2993 !$omp parallel do private(i,j)
2994  DO j=jsta,jend
2995  DO i=ista,iend
2996  IF(cuprec(i,j) < spval)THEN
2997  grid1(i,j) = cuprec(i,j)*1000.
2998  ELSE
2999  grid1(i,j) = spval
3000  ENDIF
3001  ENDDO
3002  ENDDO
3003  END IF
3004 ! write(6,*) 'call gribit...convective precip'
3005  if(grib=='grib2') then
3006  cfld=cfld+1
3007  fld_info(cfld)%ifld=iavblfld(iget(033))
3008  fld_info(cfld)%ntrange=1
3009  fld_info(cfld)%tinvstat=ifhr-id(18)
3010 !$omp parallel do private(i,j,ii,jj)
3011  do j=1,jend-jsta+1
3012  jj = jsta+j-1
3013  do i=1,iend-ista+1
3014  ii = ista+i-1
3015  datapd(i,j,cfld) = grid1(ii,jj)
3016  enddo
3017  enddo
3018 !! add continuous bucket
3019 ! if(MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R') then
3020 ! cfld=cfld+1
3021 ! fld_info(cfld)%ifld=IAVBLFLD(IGET(033))
3022 ! fld_info(cfld)%ntrange=1
3023 ! fld_info(cfld)%tinvstat=IFHR
3024 ! do j=1,jend-jsta+1
3025 ! jj = jsta+j-1
3026 ! do i=1,im
3027 ! datapd(i,j,cfld) = GRID2(i,jj)
3028 ! enddo
3029 ! enddo
3030 ! endif
3031  endif
3032  ENDIF
3033 
3034 ! CONTINOUS ACCUMULATED CONVECTIVE PRECIPITATION.
3035  IF (iget(418)>0) THEN
3036  id(1:25) = 0
3037  itprec = nint(tprec)
3038 !mp
3039  if (itprec /= 0) then
3040  ifincr = mod(ifhr,itprec)
3041  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3042  else
3043  ifincr = 0
3044  endif
3045 !mp
3046  id(18) = 0
3047  id(19) = ifhr
3048  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3049  id(20) = 4
3050  IF (ifincr==0) THEN
3051  id(18) = ifhr-itprec
3052  ELSE
3053  id(18) = ifhr-ifincr
3054  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3055  ENDIF
3056  IF (id(18)<0) id(18) = 0
3057  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
3058 ! Chuang 3/29/2018: add continuous bucket
3059 !$omp parallel do private(i,j)
3060  DO j=jsta,jend
3061  DO i=ista,iend
3062  IF(avgcprate_cont(i,j) < spval)THEN
3063  grid2(i,j) = avgcprate_cont(i,j)*float(ifhr)*3600.*1000./dtq2
3064  ELSE
3065  grid2(i,j) = spval
3066  END IF
3067  ENDDO
3068  ENDDO
3069  ENDIF
3070 ! write(6,*) 'call gribit...convective precip'
3071  if(grib=='grib2') then
3072 ! add continuous bucket
3073  if(modelname == 'GFS' .OR. modelname == 'FV3R') then
3074  cfld=cfld+1
3075  fld_info(cfld)%ifld=iavblfld(iget(418))
3076  fld_info(cfld)%ntrange=1
3077  fld_info(cfld)%tinvstat=ifhr
3078 !$omp parallel do private(i,j,ii,jj)
3079  do j=1,jend-jsta+1
3080  jj = jsta+j-1
3081  do i=1,iend-ista+1
3082  ii = ista+i-1
3083  datapd(i,j,cfld) = grid2(ii,jj)
3084  enddo
3085  enddo
3086  endif
3087  endif
3088  ENDIF
3089 !
3090 ! ACCUMULATED GRID-SCALE PRECIPITATION.
3091  IF (iget(034)>0) THEN
3092 
3093  id(1:25) = 0
3094  itprec = nint(tprec)
3095 !mp
3096  if (itprec /= 0) then
3097  ifincr = mod(ifhr,itprec)
3098  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3099  else
3100  ifincr = 0
3101  endif
3102 !mp
3103  id(18) = 0
3104  id(19) = ifhr
3105  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3106  id(20) = 4
3107  IF (ifincr==0) THEN
3108  id(18) = ifhr-itprec
3109  ELSE
3110  id(18) = ifhr-ifincr
3111  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3112  ENDIF
3113  IF (id(18)<0) id(18) = 0
3114  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
3115 !$omp parallel do private(i,j)
3116  DO j=jsta,jend
3117  DO i=ista,iend
3118  IF(avgcprate(i,j) < spval .AND. avgprec(i,j) < spval) then
3119  grid1(i,j) = ( avgprec(i,j) - avgcprate(i,j) ) * &
3120  float(id(19)-id(18))*3600.*1000./dtq2
3121  ELSE
3122  grid1(i,j) = spval
3123  END IF
3124  ENDDO
3125  ENDDO
3126 !! Chuang 3/29/2018: add continuous bucket
3127 ! DO J=JSTA,JEND
3128 ! DO I=ISTA,IEND
3129 ! IF(AVGCPRATE_CONT(I,J) < SPVAL .AND. AVGPREC_CONT(I,J) < SPVAL)THEN
3130 ! GRID2(I,J) = (AVGPREC_CONT(I,J) - AVGCPRATE_CONT(I,J)) &
3131 ! *FLOAT(IFHR)*3600.*1000./DTQ2
3132 ! ELSE
3133 ! GRID2(I,J) = SPVAL
3134 ! END IF
3135 ! ENDDO
3136 ! ENDDO
3137  ELSE
3138 !$omp parallel do private(i,j)
3139  DO j=jsta,jend
3140  DO i=ista,iend
3141  grid1(i,j) = ancprc(i,j)*1000.
3142  ENDDO
3143  ENDDO
3144  END IF
3145 ! write(6,*) 'call gribit...grid-scale precip'
3146  if(grib=='grib2') then
3147  cfld=cfld+1
3148  fld_info(cfld)%ifld=iavblfld(iget(034))
3149  fld_info(cfld)%ntrange=1
3150  fld_info(cfld)%tinvstat=ifhr-id(18)
3151 !$omp parallel do private(i,j,ii,jj)
3152  do j=1,jend-jsta+1
3153  jj = jsta+j-1
3154  do i=1,iend-ista+1
3155  ii = ista+i-1
3156  datapd(i,j,cfld) = grid1(ii,jj)
3157  enddo
3158  enddo
3159 !! add continuous bucket
3160 ! if(MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R') then
3161 ! cfld=cfld+1
3162 ! fld_info(cfld)%ifld=IAVBLFLD(IGET(034))
3163 ! fld_info(cfld)%ntrange=1
3164 ! fld_info(cfld)%tinvstat=IFHR
3165 ! do j=1,jend-jsta+1
3166 ! jj = jsta+j-1
3167 ! do i=1,iend-ista+1
3168 ! ii = ista+1-1
3169 ! datapd(i,j,cfld) = GRID2(ii,jj)
3170 ! enddo
3171 ! enddo
3172 ! endif
3173  endif
3174  ENDIF
3175 
3176 ! CONTINOUS ACCUMULATED GRID-SCALE PRECIPITATION.
3177  IF (iget(419)>0) THEN
3178  id(1:25) = 0
3179  itprec = nint(tprec)
3180 !mp
3181  if (itprec /= 0) then
3182  ifincr = mod(ifhr,itprec)
3183  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3184  else
3185  ifincr = 0
3186  endif
3187 !mp
3188  id(18) = 0
3189  id(19) = ifhr
3190  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3191  id(20) = 4
3192  IF (ifincr==0) THEN
3193  id(18) = ifhr-itprec
3194  ELSE
3195  id(18) = ifhr-ifincr
3196  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3197  ENDIF
3198  IF (id(18)<0) id(18) = 0
3199  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
3200 ! Chuang 3/29/2018: add continuous bucket
3201 !$omp parallel do private(i,j)
3202  DO j=jsta,jend
3203  DO i=ista,iend
3204  IF(avgcprate_cont(i,j) < spval .AND. avgprec_cont(i,j) < spval)THEN
3205  grid2(i,j) = (avgprec_cont(i,j) - avgcprate_cont(i,j)) &
3206  *float(ifhr)*3600.*1000./dtq2
3207  ELSE
3208  grid2(i,j) = spval
3209  END IF
3210  ENDDO
3211  ENDDO
3212  ENDIF
3213 ! write(6,*) 'call gribit...grid-scale precip'
3214  if(grib=='grib2') then
3215 ! add continuous bucket
3216  if(modelname == 'GFS' .OR. modelname == 'FV3R') then
3217  cfld=cfld+1
3218  fld_info(cfld)%ifld=iavblfld(iget(419))
3219  fld_info(cfld)%ntrange=1
3220  fld_info(cfld)%tinvstat=ifhr
3221 !$omp parallel do private(i,j,ii,jj)
3222  do j=1,jend-jsta+1
3223  jj = jsta+j-1
3224  do i=1,iend-ista+1
3225  ii = ista+i-1
3226  datapd(i,j,cfld) = grid2(ii,jj)
3227  enddo
3228  enddo
3229  endif
3230  endif
3231  ENDIF
3232 !
3233 ! ACCUMULATED LAND SURFACE PRECIPITATION.
3234  IF (iget(256)>0) THEN
3235  grid1=spval
3236 !$omp parallel do private(i,j)
3237  DO j=jsta,jend
3238  DO i=ista,iend
3239  IF(lspa(i,j)<=-1.0e-6)THEN
3240  if(acprec(i,j)/=spval) grid1(i,j) = acprec(i,j)*1000
3241  ELSE
3242  if(lspa(i,j)/=spval) grid1(i,j) = lspa(i,j)*1000.
3243  END IF
3244  ENDDO
3245  ENDDO
3246  id(1:25) = 0
3247  itprec = nint(tprec)
3248 !mp
3249  if (itprec /= 0) then
3250  ifincr = mod(ifhr,itprec)
3251  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3252  else
3253  ifincr = 0
3254  endif
3255 !mp
3256  id(18) = 0
3257  id(19) = ifhr
3258  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3259  id(20) = 4
3260  IF (ifincr==0) THEN
3261  id(18) = ifhr-itprec
3262  ELSE
3263  id(18) = ifhr-ifincr
3264  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3265  ENDIF
3266  IF (id(18)<0) id(18) = 0
3267  id(02)= 130
3268  if(grib=='grib2') then
3269  cfld=cfld+1
3270  fld_info(cfld)%ifld=iavblfld(iget(256))
3271  fld_info(cfld)%ntrange=1
3272  fld_info(cfld)%tinvstat=ifhr-id(18)
3273 !$omp parallel do private(i,j,ii,jj)
3274  do j=1,jend-jsta+1
3275  jj = jsta+j-1
3276  do i=1,iend-ista+1
3277  ii = ista+i-1
3278  datapd(i,j,cfld) = grid1(ii,jj)
3279  enddo
3280  enddo
3281  endif
3282  ENDIF
3283 !
3284 ! ACCUMULATED SNOWFALL.
3285  IF (iget(035)>0) THEN
3286 !$omp parallel do private(i,j)
3287  DO j=jsta,jend
3288  DO i=ista,iend
3289 ! GRID1(I,J) = ACSNOW(I,J)*1000.
3290  grid1(i,j) = acsnow(i,j)
3291  ENDDO
3292  ENDDO
3293  id(1:25) = 0
3294  itprec = nint(tprec)
3295 !mp
3296  if (itprec /= 0) then
3297  ifincr = mod(ifhr,itprec)
3298  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3299  else
3300  ifincr = 0
3301  endif
3302 !mp
3303  id(18) = 0
3304  id(19) = ifhr
3305  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3306  id(20) = 4
3307  IF (ifincr==0) THEN
3308  id(18) = ifhr-itprec
3309  ELSE
3310  id(18) = ifhr-ifincr
3311  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3312  ENDIF
3313  IF (id(18)<0) id(18) = 0
3314  if(grib=='grib2') then
3315  cfld=cfld+1
3316  fld_info(cfld)%ifld=iavblfld(iget(035))
3317  fld_info(cfld)%ntrange=1
3318  fld_info(cfld)%tinvstat=ifhr-id(18)
3319 !$omp parallel do private(i,j,ii,jj)
3320  do j=1,jend-jsta+1
3321  jj = jsta+j-1
3322  do i=1,iend-ista+1
3323  ii = ista+i-1
3324  datapd(i,j,cfld) = grid1(ii,jj)
3325  enddo
3326  enddo
3327  endif
3328  ENDIF
3329 !
3330 ! ACCUMULATED GRAUPEL.
3331  IF (iget(746)>0) THEN
3332 !$omp parallel do private(i,j)
3333  DO j=jsta,jend
3334  DO i=ista,iend
3335  grid1(i,j) = acgraup(i,j)
3336  ENDDO
3337  ENDDO
3338  id(1:25) = 0
3339  itprec = nint(tprec)
3340 !mp
3341  if (itprec /= 0) then
3342  ifincr = mod(ifhr,itprec)
3343  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3344  else
3345  ifincr = 0
3346  endif
3347 !mp
3348  id(18) = 0
3349  id(19) = ifhr
3350  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3351  id(20) = 4
3352  IF (ifincr==0) THEN
3353  id(18) = ifhr-itprec
3354  ELSE
3355  id(18) = ifhr-ifincr
3356  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3357  ENDIF
3358  IF (id(18)<0) id(18) = 0
3359  if(grib=='grib2') then
3360  cfld=cfld+1
3361  fld_info(cfld)%ifld=iavblfld(iget(746))
3362  fld_info(cfld)%ntrange=1
3363  fld_info(cfld)%tinvstat=ifhr-id(18)
3364 !$omp parallel do private(i,j,ii,jj)
3365  do j=1,jend-jsta+1
3366  jj = jsta+j-1
3367  do i=1,iend-ista+1
3368  ii = ista+i-1
3369  datapd(i,j,cfld) = grid1(ii,jj)
3370  enddo
3371  enddo
3372  endif
3373  ENDIF
3374 !
3375 ! ACCUMULATED FREEZING RAIN.
3376  IF (iget(782)>0) THEN
3377 !$omp parallel do private(i,j)
3378  DO j=jsta,jend
3379  DO i=ista,iend
3380  grid1(i,j) = acfrain(i,j)
3381  ENDDO
3382  ENDDO
3383  id(1:25) = 0
3384  itprec = nint(tprec)
3385 !mp
3386  if (itprec /= 0) then
3387  ifincr = mod(ifhr,itprec)
3388  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3389  else
3390  ifincr = 0
3391  endif
3392 !mp
3393  id(18) = 0
3394  id(19) = ifhr
3395  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3396  id(20) = 4
3397  IF (ifincr==0) THEN
3398  id(18) = ifhr-itprec
3399  ELSE
3400  id(18) = ifhr-ifincr
3401  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3402  ENDIF
3403  IF (id(18)<0) id(18) = 0
3404  if(grib=='grib2') then
3405  cfld=cfld+1
3406  fld_info(cfld)%ifld=iavblfld(iget(782))
3407  fld_info(cfld)%ntrange=1
3408  fld_info(cfld)%tinvstat=ifhr-id(18)
3409 !$omp parallel do private(i,j,ii,jj)
3410  do j=1,jend-jsta+1
3411  jj = jsta+j-1
3412  do i=1,iend-ista+1
3413  ii = ista+i-1
3414  datapd(i,j,cfld) = grid1(ii,jj)
3415  enddo
3416  enddo
3417  endif
3418  ENDIF
3419 !
3420 ! ACCUMULATED SNOW MELT.
3421  IF (iget(121)>0) THEN
3422 !$omp parallel do private(i,j)
3423  DO j=jsta,jend
3424  DO i=ista,iend
3425 ! GRID1(I,J) = ACSNOM(I,J)*1000.
3426  grid1(i,j) = acsnom(i,j)
3427  ENDDO
3428  ENDDO
3429  id(1:25) = 0
3430  itprec = nint(tprec)
3431 !mp
3432  if (itprec /= 0) then
3433  ifincr = mod(ifhr,itprec)
3434  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3435  else
3436  ifincr = 0
3437  endif
3438 !mp
3439  id(18) = 0
3440  id(19) = ifhr
3441  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3442  id(20) = 4
3443  IF (ifincr==0) THEN
3444  id(18) = ifhr-itprec
3445  ELSE
3446  id(18) = ifhr-ifincr
3447  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3448  ENDIF
3449  IF (id(18)<0) id(18) = 0
3450  if(grib=='grib2') then
3451  cfld=cfld+1
3452  fld_info(cfld)%ifld=iavblfld(iget(121))
3453  fld_info(cfld)%ntrange=1
3454  fld_info(cfld)%tinvstat=ifhr-id(18)
3455 !$omp parallel do private(i,j,ii,jj)
3456  do j=1,jend-jsta+1
3457  jj = jsta+j-1
3458  do i=1,iend-ista+1
3459  ii = ista+i-1
3460  datapd(i,j,cfld) = grid1(ii,jj)
3461  enddo
3462  enddo
3463  endif
3464  ENDIF
3465 !
3466 ! ACCUMULATED SNOWFALL RATE
3467  IF (iget(405)>0) THEN
3468 !$omp parallel do private(i,j)
3469  DO j=jsta,jend
3470  DO i=ista,iend
3471  grid1(i,j) = snowfall(i,j)
3472  ENDDO
3473  ENDDO
3474  id(1:25) = 0
3475  itprec = nint(tprec)
3476 !mp
3477  if (itprec /= 0) then
3478  ifincr = mod(ifhr,itprec)
3479  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3480  else
3481  ifincr = 0
3482  endif
3483 !mp
3484  id(18) = 0
3485  id(19) = ifhr
3486  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3487  id(20) = 4
3488  IF (ifincr==0) THEN
3489  id(18) = ifhr-itprec
3490  ELSE
3491  id(18) = ifhr-ifincr
3492  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3493  ENDIF
3494  IF (id(18)<0) id(18) = 0
3495  IF(itprec < 0)id(1:25)=0
3496  if(grib=='grib2') then
3497  cfld=cfld+1
3498  fld_info(cfld)%ifld=iavblfld(iget(405))
3499  fld_info(cfld)%ntrange=1
3500  fld_info(cfld)%tinvstat=ifhr-id(18)
3501 !$omp parallel do private(i,j,ii,jj)
3502  do j=1,jend-jsta+1
3503  jj = jsta+j-1
3504  do i=1,iend-ista+1
3505  ii = ista+i-1
3506  datapd(i,j,cfld) = grid1(ii,jj)
3507  enddo
3508  enddo
3509  endif
3510  ENDIF
3511 !
3512 ! ACCUMULATED STORM SURFACE RUNOFF.
3513  IF (iget(122)>0) THEN
3514 !$omp parallel do private(i,j)
3515  DO j=jsta,jend
3516  DO i=ista,iend
3517 ! GRID1(I,J) = SSROFF(I,J)*1000.
3518  grid1(i,j) = ssroff(i,j)
3519  ENDDO
3520  ENDDO
3521  id(1:25) = 0
3522  itprec = nint(tprec)
3523 !mp
3524  if (itprec /= 0) then
3525  ifincr = mod(ifhr,itprec)
3526  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3527  else
3528  ifincr = 0
3529  endif
3530 !mp
3531  id(18) = 0
3532  id(19) = ifhr
3533  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3534  id(20) = 4
3535  IF (ifincr==0) THEN
3536  id(18) = ifhr-itprec
3537  ELSE
3538  id(18) = ifhr-ifincr
3539  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3540  ENDIF
3541  IF (id(18)<0) id(18) = 0
3542 ! 1-HR RUNOFF ACCUMULATIONS IN RR
3543  IF (modelname=='RAPR') THEN
3544  IF (ifhr > 0) THEN
3545  id(18)=ifhr-1
3546  ELSE
3547  id(18)=0
3548  ENDIF
3549  ENDIF
3550  if(grib=='grib2') then
3551  cfld=cfld+1
3552  fld_info(cfld)%ifld=iavblfld(iget(122))
3553  fld_info(cfld)%ntrange=1
3554  fld_info(cfld)%tinvstat=ifhr-id(18)
3555 !$omp parallel do private(i,j,ii,jj)
3556  do j=1,jend-jsta+1
3557  jj = jsta+j-1
3558  do i=1,iend-ista+1
3559  ii = ista+i-1
3560  datapd(i,j,cfld) = grid1(ii,jj)
3561  enddo
3562  enddo
3563  endif
3564  ENDIF
3565 !
3566 ! ACCUMULATED BASEFLOW-GROUNDWATER RUNOFF.
3567  IF (iget(123)>0) THEN
3568 !$omp parallel do private(i,j)
3569  DO j=jsta,jend
3570  DO i=ista,iend
3571 ! GRID1(I,J) = BGROFF(I,J)*1000.
3572  grid1(i,j) = bgroff(i,j)
3573  ENDDO
3574  ENDDO
3575  id(1:25) = 0
3576  itprec = nint(tprec)
3577 !mp
3578  if (itprec /= 0) then
3579  ifincr = mod(ifhr,itprec)
3580  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3581  else
3582  ifincr = 0
3583  endif
3584 !mp
3585  id(18) = ifhr - 1
3586  id(19) = ifhr
3587  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3588  id(20) = 4
3589  IF (ifincr==0) THEN
3590  id(18) = ifhr-itprec
3591  ELSE
3592  id(18) = ifhr-ifincr
3593  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3594  ENDIF
3595  IF (id(18)<0) id(18) = 0
3596 ! 1-HR RUNOFF ACCUMULATIONS IN RR
3597  IF (modelname=='RAPR') THEN
3598  IF (ifhr > 0) THEN
3599  id(18)=ifhr-1
3600  ELSE
3601  id(18)=0
3602  ENDIF
3603  ENDIF
3604  if(grib=='grib2') then
3605  cfld=cfld+1
3606  fld_info(cfld)%ifld=iavblfld(iget(123))
3607  fld_info(cfld)%ntrange=1
3608  fld_info(cfld)%tinvstat=ifhr-id(18)
3609 !$omp parallel do private(i,j,ii,jj)
3610  do j=1,jend-jsta+1
3611  jj = jsta+j-1
3612  do i=1,iend-ista+1
3613  ii = ista+i-1
3614  datapd(i,j,cfld) = grid1(ii,jj)
3615  enddo
3616  enddo
3617  endif
3618  ENDIF
3619 !
3620 ! ACCUMULATED WATER RUNOFF.
3621  IF (iget(343)>0) THEN
3622 !$omp parallel do private(i,j)
3623  DO j=jsta,jend
3624  DO i=ista,iend
3625  grid1(i,j) = runoff(i,j)
3626  ENDDO
3627  ENDDO
3628  id(1:25) = 0
3629  itprec = nint(tprec)
3630 ! GFS starts to use continuous bucket for precipitation only
3631 ! so have to change water runoff to use different bucket
3632  if(modelname == 'GFS')itprec=nint(tmaxmin)
3633 !mp
3634  if (itprec /= 0) then
3635  ifincr = mod(ifhr,itprec)
3636  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3637  else
3638  ifincr = 0
3639  endif
3640 !mp
3641  id(18) = 0
3642  id(19) = ifhr
3643  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3644  id(20) = 4
3645  IF (ifincr==0) THEN
3646  id(18) = ifhr-itprec
3647  ELSE
3648  id(18) = ifhr-ifincr
3649  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3650  ENDIF
3651  IF (id(18)<0) id(18) = 0
3652  if(grib=='grib2') then
3653  cfld=cfld+1
3654  fld_info(cfld)%ifld=iavblfld(iget(343))
3655  fld_info(cfld)%ntrange=1
3656  fld_info(cfld)%tinvstat=ifhr-id(18)
3657 !$omp parallel do private(i,j,ii,jj)
3658  do j=1,jend-jsta+1
3659  jj = jsta+j-1
3660  do i=1,iend-ista+1
3661  ii = ista+i-1
3662  datapd(i,j,cfld) = grid1(ii,jj)
3663  enddo
3664  enddo
3665  endif
3666  ENDIF
3667 
3668 ! PRECIPITATION BUCKETS - accumulated between output times
3669 ! 'BUCKET TOTAL PRECIP '
3670  IF (iget(434)>0.) THEN
3671 !$omp parallel do private(i,j)
3672  DO j=jsta,jend
3673  DO i=ista,iend
3674  IF (ifhr == 0) THEN
3675  grid1(i,j) = 0.0
3676  ELSE
3677  grid1(i,j) = pcp_bucket(i,j)
3678  ENDIF
3679  ENDDO
3680  ENDDO
3681  id(1:25) = 0
3682  itprec = nint(tprec)
3683 !mp
3684  if (itprec /= 0) then
3685  ifincr = mod(ifhr,itprec)
3686  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3687  else
3688  ifincr = 0
3689  endif
3690 !mp
3691  if(modelname=='NCAR' .OR. modelname=='RAPR') ifincr = nint(prec_acc_dt)/60
3692  id(18) = 0
3693  id(19) = ifhr
3694  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3695  id(20) = 4
3696  IF (ifincr==0) THEN
3697  id(18) = ifhr-itprec
3698  ELSE
3699  id(18) = ifhr-ifincr
3700  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3701  ENDIF
3702  IF (id(18)<0) id(18) = 0
3703  if(grib=='grib2') then
3704  cfld=cfld+1
3705  fld_info(cfld)%ifld=iavblfld(iget(434))
3706  if(itprec>0) then
3707  fld_info(cfld)%ntrange=(ifhr-id(18))/itprec
3708  else
3709  fld_info(cfld)%ntrange=0
3710  endif
3711  fld_info(cfld)%tinvstat=itprec
3712  if(fld_info(cfld)%ntrange==0) then
3713  if (ifhr==0) then
3714  fld_info(cfld)%tinvstat=0
3715  else
3716  fld_info(cfld)%tinvstat=1
3717  endif
3718  fld_info(cfld)%ntrange=1
3719  end if
3720 !$omp parallel do private(i,j,ii,jj)
3721  do j=1,jend-jsta+1
3722  jj = jsta+j-1
3723  do i=1,iend-ista+1
3724  ii = ista+i-1
3725  datapd(i,j,cfld) = grid1(ii,jj)
3726  enddo
3727  enddo
3728  endif
3729  ENDIF
3730 
3731 ! PRECIPITATION BUCKETS - accumulated between output times
3732 ! 'BUCKET CONV PRECIP '
3733  IF (iget(435)>0.) THEN
3734 !$omp parallel do private(i,j)
3735  DO j=jsta,jend
3736  DO i=ista,iend
3737  IF (ifhr == 0) THEN
3738  grid1(i,j) = 0.0
3739  ELSE
3740  grid1(i,j) = rainc_bucket(i,j)
3741  ENDIF
3742  ENDDO
3743  ENDDO
3744  id(1:25) = 0
3745  itprec = nint(tprec)
3746 !mp
3747  if (itprec /= 0) then
3748  ifincr = mod(ifhr,itprec)
3749  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3750  else
3751  ifincr = 0
3752  endif
3753 
3754  if(modelname=='NCAR' .OR. modelname=='RAPR') ifincr = nint(prec_acc_dt)/60
3755 !mp
3756  id(18) = 0
3757  id(19) = ifhr
3758  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3759  id(20) = 4
3760  IF (ifincr==0) THEN
3761  id(18) = ifhr-itprec
3762  ELSE
3763  id(18) = ifhr-ifincr
3764  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3765  ENDIF
3766  IF (id(18)<0) id(18) = 0
3767 
3768 ! print *,'IFMIN,IFHR,ITPREC',IFMIN,IFHR,ITPREC
3769  if(debugprint .and. me==0)then
3770  print *,'PREC_ACC_DT,ID(18),ID(19)',prec_acc_dt,id(18),id(19)
3771  endif
3772 
3773  if(grib=='grib2') then
3774  cfld=cfld+1
3775  fld_info(cfld)%ifld=iavblfld(iget(435))
3776  if(itprec>0) then
3777  fld_info(cfld)%ntrange=(ifhr-id(18))/itprec
3778  else
3779  fld_info(cfld)%ntrange=0
3780  endif
3781  fld_info(cfld)%tinvstat=itprec
3782  if(fld_info(cfld)%ntrange==0) then
3783  if (ifhr==0) then
3784  fld_info(cfld)%tinvstat=0
3785  else
3786  fld_info(cfld)%tinvstat=1
3787  endif
3788  fld_info(cfld)%ntrange=1
3789  end if
3790 !$omp parallel do private(i,j,ii,jj)
3791  do j=1,jend-jsta+1
3792  jj = jsta+j-1
3793  do i=1,iend-ista+1
3794  ii = ista+i-1
3795  datapd(i,j,cfld) = grid1(ii,jj)
3796  enddo
3797  enddo
3798  endif
3799  ENDIF
3800 ! PRECIPITATION BUCKETS - accumulated between output times
3801 ! 'BUCKET GRDSCALE PRCP'
3802  IF (iget(436)>0.) THEN
3803 !$omp parallel do private(i,j)
3804  DO j=jsta,jend
3805  DO i=ista,iend
3806  IF (ifhr == 0) THEN
3807  grid1(i,j) = 0.0
3808  ELSE
3809  grid1(i,j) = rainnc_bucket(i,j)
3810  ENDIF
3811  ENDDO
3812  ENDDO
3813  id(1:25) = 0
3814  itprec = nint(tprec)
3815 !mp
3816  if (itprec /= 0) then
3817  ifincr = mod(ifhr,itprec)
3818  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3819  else
3820  ifincr = 0
3821  endif
3822 !mp
3823  if(modelname=='NCAR' .OR. modelname=='RAPR') ifincr = nint(prec_acc_dt)/60
3824  id(18) = 0
3825  id(19) = ifhr
3826  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3827  id(20) = 4
3828  IF (ifincr==0) THEN
3829  id(18) = ifhr-itprec
3830  ELSE
3831  id(18) = ifhr-ifincr
3832  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3833  ENDIF
3834  IF (id(18)<0) id(18) = 0
3835  if(grib=='grib2') then
3836  cfld=cfld+1
3837  fld_info(cfld)%ifld=iavblfld(iget(436))
3838  if(itprec>0) then
3839  fld_info(cfld)%ntrange=(ifhr-id(18))/itprec
3840  else
3841  fld_info(cfld)%ntrange=0
3842  endif
3843  fld_info(cfld)%tinvstat=itprec
3844  if(fld_info(cfld)%ntrange==0) then
3845  if (ifhr==0) then
3846  fld_info(cfld)%tinvstat=0
3847  else
3848  fld_info(cfld)%tinvstat=1
3849  endif
3850  fld_info(cfld)%ntrange=1
3851  end if
3852 !$omp parallel do private(i,j,ii,jj)
3853  do j=1,jend-jsta+1
3854  jj = jsta+j-1
3855  do i=1,iend-ista+1
3856  ii = ista+i-1
3857  datapd(i,j,cfld) = grid1(ii,jj)
3858  enddo
3859  enddo
3860  endif
3861  ENDIF
3862 ! PRECIPITATION BUCKETS - accumulated between output times
3863 ! 'BUCKET SNOW PRECIP '
3864  IF (iget(437)>0.) THEN
3865 !$omp parallel do private(i,j)
3866  DO j=jsta,jend
3867  DO i=ista,iend
3868  grid1(i,j) = snow_bucket(i,j)
3869  ENDDO
3870  ENDDO
3871  id(1:25) = 0
3872  itprec = nint(tprec)
3873 !mp
3874  if (itprec /= 0) then
3875  ifincr = mod(ifhr,itprec)
3876  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3877  else
3878  ifincr = 0
3879  endif
3880 !mp
3881  if(modelname=='NCAR' .OR. modelname=='RAPR') ifincr = nint(prec_acc_dt)/60
3882  id(18) = 0
3883  id(19) = ifhr
3884  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3885  id(20) = 4
3886  IF (ifincr==0) THEN
3887  id(18) = ifhr-itprec
3888  ELSE
3889  id(18) = ifhr-ifincr
3890  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3891  ENDIF
3892  IF (id(18)<0) id(18) = 0
3893 ! if(me==0)print*,'maxval BUCKET SNOWFALL: ', maxval(GRID1)
3894  if(grib=='grib2') then
3895  cfld=cfld+1
3896  fld_info(cfld)%ifld=iavblfld(iget(437))
3897  if(itprec>0) then
3898  fld_info(cfld)%ntrange=(ifhr-id(18))/itprec
3899  else
3900  fld_info(cfld)%ntrange=0
3901  endif
3902  fld_info(cfld)%tinvstat=itprec
3903  if(fld_info(cfld)%ntrange==0) then
3904  if (ifhr==0) then
3905  fld_info(cfld)%tinvstat=0
3906  else
3907  fld_info(cfld)%tinvstat=1
3908  endif
3909  fld_info(cfld)%ntrange=1
3910  end if
3911 !$omp parallel do private(i,j,ii,jj)
3912  do j=1,jend-jsta+1
3913  jj = jsta+j-1
3914  do i=1,iend-ista+1
3915  ii = ista+i-1
3916  datapd(i,j,cfld) = grid1(ii,jj)
3917  enddo
3918  enddo
3919  endif
3920  ENDIF
3921 ! PRECIPITATION BUCKETS - accumulated between output times
3922 ! 'BUCKET GRAUPEL PRECIP '
3923  IF (iget(775)>0.) THEN
3924 !$omp parallel do private(i,j)
3925  DO j=jsta,jend
3926  DO i=ista,iend
3927  grid1(i,j) = graup_bucket(i,j)
3928  ENDDO
3929  ENDDO
3930  id(1:25) = 0
3931  itprec = nint(tprec)
3932 !mp
3933  if (itprec /= 0) then
3934  ifincr = mod(ifhr,itprec)
3935  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
3936  else
3937  ifincr = 0
3938  endif
3939 !mp
3940  if(modelname=='NCAR' .OR. modelname=='RAPR') ifincr = nint(prec_acc_dt)/60
3941  id(18) = 0
3942  id(19) = ifhr
3943  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3944  id(20) = 4
3945  IF (ifincr==0) THEN
3946  id(18) = ifhr-itprec
3947  ELSE
3948  id(18) = ifhr-ifincr
3949  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3950  ENDIF
3951  IF (id(18)<0) id(18) = 0
3952 ! print*,'maxval BUCKET GRAUPEL: ', maxval(GRID1)
3953  if(grib=='grib2') then
3954  cfld=cfld+1
3955  fld_info(cfld)%ifld=iavblfld(iget(775))
3956  if(itprec>0) then
3957  fld_info(cfld)%ntrange=(ifhr-id(18))/itprec
3958  else
3959  fld_info(cfld)%ntrange=0
3960  endif
3961  fld_info(cfld)%tinvstat=itprec
3962  if(fld_info(cfld)%ntrange==0) then
3963  if (ifhr==0) then
3964  fld_info(cfld)%tinvstat=0
3965  else
3966  fld_info(cfld)%tinvstat=1
3967  endif
3968  fld_info(cfld)%ntrange=1
3969  end if
3970 !$omp parallel do private(i,j,ii,jj)
3971  do j=1,jend-jsta+1
3972  jj = jsta+j-1
3973  do i=1,iend-ista+1
3974  ii = ista+i-1
3975  datapd(i,j,cfld) = grid1(ii,jj)
3976  enddo
3977  enddo
3978  endif
3979  ENDIF
3980 
3981 ! ERIC JAMES: 10 JUN 2021 -- adding precip comparison to FFG
3982 ! thresholds. 913 is for 1h QPF, 914 for run total QPF.
3983  IF (iget(913).GT.0) THEN
3984  ffgfile='ffg_01h.grib2'
3985  call qpf_comp(913,ffgfile,1)
3986  ENDIF
3987  IF (iget(914).GT.0) THEN
3988  IF (ifhr .EQ. 1) THEN
3989  ffgfile='ffg_01h.grib2'
3990  call qpf_comp(914,ffgfile,1)
3991  ELSEIF (ifhr .EQ. 3) THEN
3992  ffgfile='ffg_03h.grib2'
3993  call qpf_comp(914,ffgfile,3)
3994  ELSEIF (ifhr .EQ. 6) THEN
3995  ffgfile='ffg_06h.grib2'
3996  call qpf_comp(914,ffgfile,6)
3997  ELSEIF (ifhr .EQ. 12) THEN
3998  ffgfile='ffg_12h.grib2'
3999  call qpf_comp(914,ffgfile,12)
4000  ELSE
4001  ffgfile='ffg_01h.grib2'
4002  call qpf_comp(914,ffgfile,0)
4003  ENDIF
4004  ENDIF
4005 
4006 ! ERIC JAMES: 8 OCT 2021 -- adding precip comparison to ARI
4007 ! thresholds. 915 is for 1h QPF, 916 for run total QPF.
4008 
4009  IF (iget(915).GT.0) THEN
4010  arifile='ari2y_01h.grib2'
4011  call qpf_comp(915,arifile,1)
4012  ENDIF
4013  IF (iget(916).GT.0) THEN
4014  IF (ifhr .EQ. 1) THEN
4015  arifile='ari2y_01h.grib2'
4016  call qpf_comp(916,arifile,1)
4017  ELSEIF (ifhr .EQ. 3) THEN
4018  arifile='ari2y_03h.grib2'
4019  call qpf_comp(916,arifile,3)
4020  ELSEIF (ifhr .EQ. 6) THEN
4021  arifile='ari2y_06h.grib2'
4022  call qpf_comp(916,arifile,6)
4023  ELSEIF (ifhr .EQ. 12) THEN
4024  arifile='ari2y_12h.grib2'
4025  call qpf_comp(916,arifile,12)
4026  ELSEIF (ifhr .EQ. 24) THEN
4027  arifile='ari2y_24h.grib2'
4028  call qpf_comp(916,arifile,24)
4029  ELSE
4030  arifile='ari2y_01h.grib2'
4031  call qpf_comp(916,arifile,0)
4032  ENDIF
4033  ENDIF
4034 
4035  IF (iget(917).GT.0) THEN
4036  arifile='ari5y_01h.grib2'
4037  call qpf_comp(917,arifile,1)
4038  ENDIF
4039  IF (iget(918).GT.0) THEN
4040  IF (ifhr .EQ. 1) THEN
4041  arifile='ari5y_01h.grib2'
4042  call qpf_comp(918,arifile,1)
4043  ELSEIF (ifhr .EQ. 3) THEN
4044  arifile='ari5y_03h.grib2'
4045  call qpf_comp(918,arifile,3)
4046  ELSEIF (ifhr .EQ. 6) THEN
4047  arifile='ari5y_06h.grib2'
4048  call qpf_comp(918,arifile,6)
4049  ELSEIF (ifhr .EQ. 12) THEN
4050  arifile='ari5y_12h.grib2'
4051  call qpf_comp(918,arifile,12)
4052  ELSEIF (ifhr .EQ. 24) THEN
4053  arifile='ari5y_24h.grib2'
4054  call qpf_comp(918,arifile,24)
4055  ELSE
4056  arifile='ari5y_01h.grib2'
4057  call qpf_comp(918,arifile,0)
4058  ENDIF
4059  ENDIF
4060 
4061  IF (iget(919).GT.0) THEN
4062  arifile='ari10y_01h.grib2'
4063  call qpf_comp(919,arifile,1)
4064  ENDIF
4065  IF (iget(920).GT.0) THEN
4066  IF (ifhr .EQ. 1) THEN
4067  arifile='ari10y_01h.grib2'
4068  call qpf_comp(920,arifile,1)
4069  ELSEIF (ifhr .EQ. 3) THEN
4070  arifile='ari10y_03h.grib2'
4071  call qpf_comp(920,arifile,3)
4072  ELSEIF (ifhr .EQ. 6) THEN
4073  arifile='ari10y_06h.grib2'
4074  call qpf_comp(920,arifile,6)
4075  ELSEIF (ifhr .EQ. 12) THEN
4076  arifile='ari10y_12h.grib2'
4077  call qpf_comp(920,arifile,12)
4078  ELSEIF (ifhr .EQ. 24) THEN
4079  arifile='ari10y_24h.grib2'
4080  call qpf_comp(920,arifile,24)
4081  ELSE
4082  arifile='ari10y_01h.grib2'
4083  call qpf_comp(920,arifile,0)
4084  ENDIF
4085  ENDIF
4086 
4087  IF (iget(921).GT.0) THEN
4088  arifile='ari100y_01h.grib2'
4089  call qpf_comp(921,arifile,1)
4090  ENDIF
4091  IF (iget(922).GT.0) THEN
4092  IF (ifhr .EQ. 1) THEN
4093  arifile='ari100y_01h.grib2'
4094  call qpf_comp(922,arifile,1)
4095  ELSEIF (ifhr .EQ. 3) THEN
4096  arifile='ari100y_03h.grib2'
4097  call qpf_comp(922,arifile,3)
4098  ELSEIF (ifhr .EQ. 6) THEN
4099  arifile='ari100y_06h.grib2'
4100  call qpf_comp(922,arifile,6)
4101  ELSEIF (ifhr .EQ. 12) THEN
4102  arifile='ari100y_12h.grib2'
4103  call qpf_comp(922,arifile,12)
4104  ELSEIF (ifhr .EQ. 24) THEN
4105  arifile='ari100y_24h.grib2'
4106  call qpf_comp(922,arifile,24)
4107  ELSE
4108  arifile='ari100y_01h.grib2'
4109  call qpf_comp(922,arifile,0)
4110  ENDIF
4111  ENDIF
4112 
4113 ! ERIC JAMES: 10 APR 2019 -- adding 15min precip output for RAP/HRRR
4114 ! PRECIPITATION BUCKETS - accumulated between output times
4115 ! 'BUCKET1 TOTAL PRECIP '
4116  IF (iget(526)>0.) THEN
4117 !$omp parallel do private(i,j)
4118  DO j=jsta,jend
4119  DO i=ista,iend
4120  IF (ifhr == 0 .AND. ifmin == 0) THEN
4121  grid1(i,j) = 0.0
4122  ELSE
4123  grid1(i,j) = pcp_bucket1(i,j)
4124  ENDIF
4125  ENDDO
4126  ENDDO
4127  ifincr = nint(prec_acc_dt1)
4128  if(grib=='grib2') then
4129  cfld=cfld+1
4130  fld_info(cfld)%ifld=iavblfld(iget(518))
4131  if(fld_info(cfld)%ntrange==0) then
4132  if (ifhr==0 .and. ifmin==0) then
4133  fld_info(cfld)%tinvstat=0
4134  else
4135  fld_info(cfld)%tinvstat=ifincr
4136  endif
4137  fld_info(cfld)%ntrange=1
4138  end if
4139 !$omp parallel do private(i,j,ii,jj)
4140  do j=1,jend-jsta+1
4141  jj = jsta+j-1
4142  do i=1,iend-ista+1
4143  ii = ista+i-1
4144  datapd(i,j,cfld) = grid1(ii,jj)
4145  enddo
4146  enddo
4147  endif
4148  ENDIF
4149 ! 'BUCKET1 CONV PRECIP '
4150  IF (iget(527)>0.) THEN
4151 !$omp parallel do private(i,j)
4152  DO j=jsta,jend
4153  DO i=ista,iend
4154  IF (ifhr == 0 .AND. ifmin == 0) THEN
4155  grid1(i,j) = 0.0
4156  ELSE
4157  grid1(i,j) = rainc_bucket1(i,j)
4158  ENDIF
4159  ENDDO
4160  ENDDO
4161  ifincr = nint(prec_acc_dt1)
4162  if(grib=='grib2') then
4163  cfld=cfld+1
4164  fld_info(cfld)%ifld=iavblfld(iget(519))
4165  if(fld_info(cfld)%ntrange==0) then
4166  if (ifhr==0 .and. ifmin==0) then
4167  fld_info(cfld)%tinvstat=0
4168  else
4169  fld_info(cfld)%tinvstat=ifincr
4170  endif
4171  fld_info(cfld)%ntrange=1
4172  end if
4173 !$omp parallel do private(i,j,ii,jj)
4174  do j=1,jend-jsta+1
4175  jj = jsta+j-1
4176  do i=1,iend-ista+1
4177  ii = ista+i-1
4178  datapd(i,j,cfld) = grid1(ii,jj)
4179  enddo
4180  enddo
4181  endif
4182  ENDIF
4183 ! 'BUCKET1 GRDSCALE PRCP'
4184  IF (iget(528)>0.) THEN
4185 !$omp parallel do private(i,j)
4186  DO j=jsta,jend
4187  DO i=ista,iend
4188  IF (ifhr == 0 .AND. ifmin == 0) THEN
4189  grid1(i,j) = 0.0
4190  ELSE
4191  grid1(i,j) = rainnc_bucket1(i,j)
4192  ENDIF
4193  ENDDO
4194  ENDDO
4195  ifincr = nint(prec_acc_dt1)
4196  if(grib=='grib2') then
4197  cfld=cfld+1
4198  fld_info(cfld)%ifld=iavblfld(iget(520))
4199  if(fld_info(cfld)%ntrange==0) then
4200  if (ifhr==0 .and. ifmin==0) then
4201  fld_info(cfld)%tinvstat=0
4202  else
4203  fld_info(cfld)%tinvstat=ifincr
4204  endif
4205  fld_info(cfld)%ntrange=1
4206  end if
4207 !$omp parallel do private(i,j,ii,jj)
4208  do j=1,jend-jsta+1
4209  jj = jsta+j-1
4210  do i=1,iend-ista+1
4211  ii = ista+i-1
4212  datapd(i,j,cfld) = grid1(ii,jj)
4213  enddo
4214  enddo
4215  endif
4216  ENDIF
4217 ! 'BUCKET1 SNOW PRECIP '
4218  IF (iget(529)>0.) THEN
4219 !$omp parallel do private(i,j)
4220  DO j=jsta,jend
4221  DO i=ista,iend
4222  IF (ifhr == 0 .AND. ifmin == 0) THEN
4223  grid1(i,j) = 0.0
4224  ELSE
4225  grid1(i,j) = snow_bucket1(i,j)
4226  ENDIF
4227  ENDDO
4228  ENDDO
4229  ifincr = nint(prec_acc_dt1)
4230 ! if(me==0)print*,'maxval BUCKET1 SNOWFALL: ', maxval(GRID1)
4231  if(grib=='grib2') then
4232  cfld=cfld+1
4233  fld_info(cfld)%ifld=iavblfld(iget(521))
4234  if(fld_info(cfld)%ntrange==0) then
4235  if (ifhr==0 .and. ifmin==0) then
4236  fld_info(cfld)%tinvstat=0
4237  else
4238  fld_info(cfld)%tinvstat=ifincr
4239  endif
4240  fld_info(cfld)%ntrange=1
4241  end if
4242 !$omp parallel do private(i,j,ii,jj)
4243  do j=1,jend-jsta+1
4244  jj = jsta+j-1
4245  do i=1,iend-ista+1
4246  ii = ista+i-1
4247  datapd(i,j,cfld) = grid1(ii,jj)
4248  enddo
4249  enddo
4250  endif
4251  ENDIF
4252 ! 'BUCKET1 GRAUPEL PRECIP '
4253  IF (iget(530)>0.) THEN
4254 !$omp parallel do private(i,j)
4255  DO j=jsta,jend
4256  DO i=ista,iend
4257  IF (ifhr == 0 .AND. ifmin == 0) THEN
4258  grid1(i,j) = 0.0
4259  ELSE
4260  grid1(i,j) = graup_bucket1(i,j)
4261  ENDIF
4262  ENDDO
4263  ENDDO
4264  ifincr = nint(prec_acc_dt1)
4265 ! print*,'maxval BUCKET1 GRAUPEL: ', maxval(GRID1)
4266  if(grib=='grib2') then
4267  cfld=cfld+1
4268  fld_info(cfld)%ifld=iavblfld(iget(522))
4269  if(fld_info(cfld)%ntrange==0) then
4270  if (ifhr==0 .and. ifmin==0) then
4271  fld_info(cfld)%tinvstat=0
4272  else
4273  fld_info(cfld)%tinvstat=ifincr
4274  endif
4275  fld_info(cfld)%ntrange=1
4276  end if
4277 !$omp parallel do private(i,j,ii,jj)
4278  do j=1,jend-jsta+1
4279  jj = jsta+j-1
4280  do i=1,iend-ista+1
4281  ii = ista+i-1
4282  datapd(i,j,cfld) = grid1(ii,jj)
4283  enddo
4284  enddo
4285  endif
4286  ENDIF
4287 !
4288 ! INSTANTANEOUS PRECIPITATION TYPE.
4289 ! print *,'in surfce,iget(160)=',iget(160),'iget(247)=',iget(247)
4290  IF (iget(160)>0 .OR.(iget(247)>0)) THEN
4291 
4292  allocate(sleet(ista:iend,jsta:jend,nalg), rain(ista:iend,jsta:jend,nalg), &
4293  freezr(ista:iend,jsta:jend,nalg), snow(ista:iend,jsta:jend,nalg))
4294  allocate(zwet(ista:iend,jsta:jend))
4295  CALL calwxt_post(t,q,pmid,pint,htm,lmh,prec,zint,iwx1,zwet)
4296 ! write(0,*)' after first CALWXT_POST'
4297 
4298 
4299  IF (iget(160)>0) THEN
4300 !$omp parallel do private(i,j,iwx)
4301  DO j=jsta,jend
4302  DO i=ista,iend
4303  IF(zwet(i,j)<spval)THEN
4304  iwx = iwx1(i,j)
4305  snow(i,j,1) = mod(iwx,2)
4306  sleet(i,j,1) = mod(iwx,4)/2
4307  freezr(i,j,1) = mod(iwx,8)/4
4308  rain(i,j,1) = iwx/8
4309  ELSE
4310  snow(i,j,1) = spval
4311  sleet(i,j,1) = spval
4312  freezr(i,j,1) = spval
4313  rain(i,j,1) = spval
4314  ENDIF
4315  ENDDO
4316  ENDDO
4317  ENDIF
4318 !
4319 ! LOWEST WET BULB ZERO HEIGHT
4320  IF (iget(247)>0) THEN
4321  DO j=jsta,jend
4322  DO i=ista,iend
4323  grid1(i,j) = zwet(i,j)
4324  ENDDO
4325  ENDDO
4326  if(grib=='grib2') then
4327  cfld=cfld+1
4328  fld_info(cfld)%ifld=iavblfld(iget(247))
4329 !$omp parallel do private(i,j,ii,jj)
4330  do j=1,jend-jsta+1
4331  jj = jsta+j-1
4332  do i=1,iend-ista+1
4333  ii = ista+i-1
4334  datapd(i,j,cfld) = grid1(ii,jj)
4335  enddo
4336  enddo
4337  endif
4338  ENDIF
4339 
4340 ! DOMINANT PRECIPITATION TYPE
4341 !GSM IF DOMINANT PRECIP TYPE IS REQUESTED, 4 MORE ALGORITHMS
4342 !GSM WILL BE CALLED. THE TALLIES ARE THEN SUMMED IN
4343 !GSM CALWXT_DOMINANT
4344 
4345  IF (iget(160)>0) THEN
4346 ! RAMER ALGORITHM
4347  CALL calwxt_ramer_post(t,q,pmid,pint,lmh,prec,iwx1)
4348 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4349 
4350 ! DECOMPOSE IWX1 ARRAY
4351 !
4352 !$omp parallel do private(i,j,iwx)
4353  DO j=jsta,jend
4354  DO i=ista,iend
4355  iwx = iwx1(i,j)
4356  snow(i,j,2) = mod(iwx,2)
4357  sleet(i,j,2) = mod(iwx,4)/2
4358  freezr(i,j,2) = mod(iwx,8)/4
4359  rain(i,j,2) = iwx/8
4360  ENDDO
4361  ENDDO
4362 
4363 ! BOURGOUIN ALGORITHM
4364  iseed=44641*(int(sdat(1)-1)*24*31+int(sdat(2))*24+ihrst)+ &
4365  & mod(ifhr*60+ifmin,44641)+4357
4366 ! write(0,*)'in SURFCE,me=',me,'bef 1st CALWXT_BOURG_POST iseed=',iseed
4367  CALL calwxt_bourg_post(im,ista_2l,iend_2u,ista,iend,jm,jsta_2l,jend_2u,jsta,jend,lm,lp1,&
4368  & iseed,g,pthresh, &
4369  & t,q,pmid,pint,lmh,prec,zint,iwx1,me)
4370 ! write(0,*)'in SURFCE,me=',me,'aft 1st CALWXT_BOURG_POST'
4371 ! write(0,*)'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA),'PTHRESH=',PTHRESH
4372 
4373 ! DECOMPOSE IWX1 ARRAY
4374 !
4375 !$omp parallel do private(i,j,iwx)
4376  DO j=jsta,jend
4377  DO i=ista,iend
4378  iwx = iwx1(i,j)
4379  snow(i,j,3) = mod(iwx,2)
4380  sleet(i,j,3) = mod(iwx,4)/2
4381  freezr(i,j,3) = mod(iwx,8)/4
4382  rain(i,j,3) = iwx/8
4383  ENDDO
4384  ENDDO
4385 
4386 ! REVISED NCEP ALGORITHM
4387  CALL calwxt_revised_post(t,q,pmid,pint,htm,lmh,prec,zint,iwx1)
4388 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4389 ! DECOMPOSE IWX1 ARRAY
4390 !
4391 !$omp parallel do private(i,j,iwx)
4392  DO j=jsta,jend
4393  DO i=ista,iend
4394  iwx = iwx1(i,j)
4395  snow(i,j,4) = mod(iwx,2)
4396  sleet(i,j,4) = mod(iwx,4)/2
4397  freezr(i,j,4) = mod(iwx,8)/4
4398  rain(i,j,4) = iwx/8
4399  ENDDO
4400  ENDDO
4401 
4402 ! EXPLICIT ALGORITHM (UNDER 18 NOT ADMITTED WITHOUT PARENT OR GUARDIAN)
4403 
4404  IF(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
4405  CALL calwxt_explicit_post(lmh,ths,pmid,prec,sr,f_rimef,iwx1)
4406  else
4407 !$omp parallel do private(i,j)
4408  DO j=jsta,jend
4409  DO i=ista,iend
4410  iwx1(i,j) = 0
4411  ENDDO
4412  ENDDO
4413  end if
4414 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4415 ! DECOMPOSE IWX1 ARRAY
4416 !
4417 !$omp parallel do private(i,j,iwx)
4418  DO j=jsta,jend
4419  DO i=ista,iend
4420  iwx = iwx1(i,j)
4421  snow(i,j,5) = mod(iwx,2)
4422  sleet(i,j,5) = mod(iwx,4)/2
4423  freezr(i,j,5) = mod(iwx,8)/4
4424  rain(i,j,5) = iwx/8
4425  ENDDO
4426  ENDDO
4427 
4428  allocate(domr(ista:iend,jsta:jend), doms(ista:iend,jsta:jend), &
4429  domzr(ista:iend,jsta:jend), domip(ista:iend,jsta:jend))
4430  CALL calwxt_dominant_post(prec(ista_2l,jsta_2l),rain,freezr,sleet,snow, &
4431  domr,domzr,domip,doms)
4432 ! if ( me==0) print *,'after CALWXT_DOMINANT, no avrg'
4433 ! SNOW.
4434  grid1 = spval
4435 !$omp parallel do private(i,j)
4436  DO j=jsta,jend
4437  DO i=ista,iend
4438  if(prec(i,j) /= spval) grid1(i,j) = doms(i,j)
4439  ENDDO
4440  ENDDO
4441  if(grib=='grib2') then
4442  cfld=cfld+1
4443  fld_info(cfld)%ifld=iavblfld(iget(551))
4444 !$omp parallel do private(i,j,ii,jj)
4445  do j=1,jend-jsta+1
4446  jj = jsta+j-1
4447  do i=1,iend-ista+1
4448  ii = ista+i-1
4449  datapd(i,j,cfld) = grid1(ii,jj)
4450  enddo
4451  enddo
4452  endif
4453 ! ICE PELLETS.
4454  grid1=spval
4455 !$omp parallel do private(i,j)
4456  DO j=jsta,jend
4457  DO i=ista,iend
4458  if(prec(i,j)/=spval) grid1(i,j) = domip(i,j)
4459  ENDDO
4460  ENDDO
4461  if(grib=='grib2') then
4462  cfld=cfld+1
4463  fld_info(cfld)%ifld=iavblfld(iget(552))
4464 !$omp parallel do private(i,j,ii,jj)
4465  do j=1,jend-jsta+1
4466  jj = jsta+j-1
4467  do i=1,iend-ista+1
4468  ii = ista+i-1
4469  datapd(i,j,cfld) = grid1(ii,jj)
4470  enddo
4471  enddo
4472  endif
4473 ! FREEZING RAIN.
4474  grid1=spval
4475 !$omp parallel do private(i,j)
4476  DO j=jsta,jend
4477  DO i=ista,iend
4478 ! if (DOMZR(I,J) == 1) THEN
4479 ! PSFC(I,J)=PINT(I,J,NINT(LMH(I,J))+1)
4480 ! print *, 'aha ', I, J, PSFC(I,J)
4481 ! print *, FREEZR(I,J,1), FREEZR(I,J,2),
4482 ! * FREEZR(I,J,3), FREEZR(I,J,4), FREEZR(I,J,5)
4483 ! endif
4484  if(prec(i,j)/=spval)grid1(i,j) = domzr(i,j)
4485  ENDDO
4486  ENDDO
4487  if(grib=='grib2') then
4488  cfld=cfld+1
4489  fld_info(cfld)%ifld=iavblfld(iget(553))
4490 !$omp parallel do private(i,j,ii,jj)
4491  do j=1,jend-jsta+1
4492  jj = jsta+j-1
4493  do i=1,iend-ista+1
4494  ii = ista+i-1
4495  datapd(i,j,cfld) = grid1(ii,jj)
4496  enddo
4497  enddo
4498  endif
4499 ! RAIN.
4500  grid1=spval
4501 !$omp parallel do private(i,j)
4502  DO j=jsta,jend
4503  DO i=ista,iend
4504  if(prec(i,j)/=spval)grid1(i,j) = domr(i,j)
4505  ENDDO
4506  ENDDO
4507  if(grib=='grib2') then
4508  cfld=cfld+1
4509  fld_info(cfld)%ifld=iavblfld(iget(160))
4510 !$omp parallel do private(i,j,ii,jj)
4511  do j=1,jend-jsta+1
4512  jj = jsta+j-1
4513  do i=1,iend-ista+1
4514  ii = ista+i-1
4515  datapd(i,j,cfld) = grid1(ii,jj)
4516  enddo
4517  enddo
4518  endif
4519  ENDIF
4520  ENDIF
4521 !
4522 ! TIME AVERAGED PRECIPITATION TYPE.
4523  IF (iget(317)>0) THEN
4524 
4525  if (.not. allocated(sleet)) allocate(sleet(ista:iend,jsta:jend,nalg))
4526  if (.not. allocated(rain)) allocate(rain(ista:iend,jsta:jend,nalg))
4527  if (.not. allocated(freezr)) allocate(freezr(ista:iend,jsta:jend,nalg))
4528  if (.not. allocated(snow)) allocate(snow(ista:iend,jsta:jend,nalg))
4529  if (.not. allocated(zwet)) allocate(zwet(ista:iend,jsta:jend))
4530  CALL calwxt_post(t,q,pmid,pint,htm,lmh,avgprec,zint,iwx1,zwet)
4531 
4532 !$omp parallel do private(i,j,iwx)
4533  DO j=jsta,jend
4534  DO i=ista,iend
4535  IF(zwet(i,j)<spval)THEN
4536  iwx = iwx1(i,j)
4537  snow(i,j,1) = mod(iwx,2)
4538  sleet(i,j,1) = mod(iwx,4)/2
4539  freezr(i,j,1) = mod(iwx,8)/4
4540  rain(i,j,1) = iwx/8
4541  ELSE
4542  snow(i,j,1) = spval
4543  sleet(i,j,1) = spval
4544  freezr(i,j,1) = spval
4545  rain(i,j,1) = spval
4546  ENDIF
4547  ENDDO
4548  ENDDO
4549  if (allocated(zwet)) deallocate(zwet)
4550 ! write(0,*)' after second CALWXT_POST me=',me
4551 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4552 
4553 ! DOMINANT PRECIPITATION TYPE
4554 !GSM IF DOMINANT PRECIP TYPE IS REQUESTED, 4 MORE ALGORITHMS
4555 !GSM WILL BE CALLED. THE TALLIES ARE THEN SUMMED IN
4556 !GSM CALWXT_DOMINANT
4557 
4558 ! RAMER ALGORITHM
4559  CALL calwxt_ramer_post(t,q,pmid,pint,lmh,avgprec,iwx1)
4560 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4561 
4562 ! DECOMPOSE IWX1 ARRAY
4563 !
4564 !$omp parallel do private(i,j,iwx)
4565  DO j=jsta,jend
4566  DO i=ista,iend
4567  iwx = iwx1(i,j)
4568  snow(i,j,2) = mod(iwx,2)
4569  sleet(i,j,2) = mod(iwx,4)/2
4570  freezr(i,j,2) = mod(iwx,8)/4
4571  rain(i,j,2) = iwx/8
4572  ENDDO
4573  ENDDO
4574 
4575 ! BOURGOUIN ALGORITHM
4576  iseed=44641*(int(sdat(1)-1)*24*31+int(sdat(2))*24+ihrst)+ &
4577  & mod(ifhr*60+ifmin,44641)+4357
4578 ! write(0,*)'in SURFCE,me=',me,'bef sec CALWXT_BOURG_POST'
4579  CALL calwxt_bourg_post(im,ista_2l,iend_2u,ista,iend,jm,jsta_2l,jend_2u,jsta,jend,lm,lp1,&
4580  & iseed,g,pthresh, &
4581  & t,q,pmid,pint,lmh,avgprec,zint,iwx1,me)
4582 ! write(0,*)'in SURFCE,me=',me,'aft sec CALWXT_BOURG_POST'
4583 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4584 
4585 ! DECOMPOSE IWX1 ARRAY
4586 !
4587 !$omp parallel do private(i,j,iwx)
4588  DO j=jsta,jend
4589  DO i=ista,iend
4590  iwx = iwx1(i,j)
4591  snow(i,j,3) = mod(iwx,2)
4592  sleet(i,j,3) = mod(iwx,4)/2
4593  freezr(i,j,3) = mod(iwx,8)/4
4594  rain(i,j,3) = iwx/8
4595  ENDDO
4596  ENDDO
4597 
4598 ! REVISED NCEP ALGORITHM
4599  CALL calwxt_revised_post(t,q,pmid,pint,htm,lmh,avgprec,zint,iwx1)
4600 ! write(0,*)'in SURFCE,me=',me,'aft sec CALWXT_REVISED_BOURG_POST'
4601 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4602 ! DECOMPOSE IWX1 ARRAY
4603 !
4604 !$omp parallel do private(i,j,iwx)
4605  DO j=jsta,jend
4606  DO i=ista,iend
4607  iwx = iwx1(i,j)
4608  snow(i,j,4) = mod(iwx,2)
4609  sleet(i,j,4) = mod(iwx,4)/2
4610  freezr(i,j,4) = mod(iwx,8)/4
4611  rain(i,j,4) = iwx/8
4612  ENDDO
4613  ENDDO
4614 
4615 ! EXPLICIT ALGORITHM (UNDER 18 NOT ADMITTED WITHOUT PARENT OR GUARDIAN)
4616 
4617 ! write(0,*)'in SURFCE,me=',me,'imp_physics=',imp_physics
4618  IF(imp_physics == 5)then
4619  CALL calwxt_explicit_post(lmh,ths,pmid,avgprec,sr,f_rimef,iwx1)
4620  else
4621 !$omp parallel do private(i,j)
4622  DO j=jsta,jend
4623  DO i=ista,iend
4624  iwx1(i,j) = 0
4625  ENDDO
4626  ENDDO
4627  end if
4628 ! print *,'in SURFCE,me=',me,'IWX1=',IWX1(1:30,JSTA)
4629 ! DECOMPOSE IWX1 ARRAY
4630 !
4631 !$omp parallel do private(i,j,iwx)
4632  DO j=jsta,jend
4633  DO i=ista,iend
4634  iwx = iwx1(i,j)
4635  snow(i,j,5) = mod(iwx,2)
4636  sleet(i,j,5) = mod(iwx,4)/2
4637  freezr(i,j,5) = mod(iwx,8)/4
4638  rain(i,j,5) = iwx/8
4639  ENDDO
4640  ENDDO
4641 
4642 ! print *,'me=',me,'before SNOW=',snow(1:10,JSTA,1:5)
4643 ! print *,'me=',me,'before RAIN=',RAIN(1:10,JSTA,1:5)
4644 ! print *,'me=',me,'before FREEZR=',FREEZR(1:10,JSTA,1:5)
4645 ! print *,'me=',me,'before SLEET=',SLEET(1:10,JSTA,1:5)
4646 
4647  if (.not. allocated(domr)) allocate(domr(ista:iend,jsta:jend))
4648  if (.not. allocated(doms)) allocate(doms(ista:iend,jsta:jend))
4649  if (.not. allocated(domzr)) allocate(domzr(ista:iend,jsta:jend))
4650  if (.not. allocated(domip)) allocate(domip(ista:iend,jsta:jend))
4651 
4652  CALL calwxt_dominant_post(avgprec,rain,freezr,sleet,snow, &
4653  domr,domzr,domip,doms)
4654 
4655  id(1:25) = 0
4656  itprec = nint(tprec)
4657 !mp
4658  if (itprec /= 0) then
4659  ifincr = mod(ifhr,itprec)
4660  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
4661  else
4662  ifincr = 0
4663  endif
4664 !mp
4665  id(18) = 0
4666  id(19) = ifhr
4667  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4668  id(20) = 3
4669  IF (ifincr==0) THEN
4670  id(18) = ifhr-itprec
4671  ELSE
4672  id(18) = ifhr-ifincr
4673  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4674  ENDIF
4675 
4676 ! TPREC,'IFHR=',IFHR,'IFMIN=',IFMIN,'IFINCR=',IFINCR,'ID=',ID
4677 ! SNOW.
4678 
4679  id(8) = 143
4680  grid1=spval
4681 !$omp parallel do private(i,j)
4682  DO j=jsta,jend
4683  DO i=ista,iend
4684  if(avgprec(i,j) /= spval) grid1(i,j) = doms(i,j)
4685  ENDDO
4686  ENDDO
4687 ! print *,'me=',me,'SNOW=',GRID1(1:10,JSTA)
4688  if(grib=='grib2') then
4689  cfld=cfld+1
4690  fld_info(cfld)%ifld=iavblfld(iget(555))
4691  if(itprec==0) then
4692  fld_info(cfld)%ntrange=0
4693  else
4694  fld_info(cfld)%ntrange=1
4695  endif
4696  fld_info(cfld)%tinvstat=ifhr-id(18)
4697 
4698 !$omp parallel do private(i,j,ii,jj)
4699  do j=1,jend-jsta+1
4700  jj = jsta+j-1
4701  do i=1,iend-ista+1
4702  ii = ista+i-1
4703  datapd(i,j,cfld) = grid1(ii,jj)
4704  enddo
4705  enddo
4706  endif
4707 ! ICE PELLETS.
4708  id(8) = 142
4709  itprec = nint(tprec)
4710 !mp
4711  if (itprec /= 0) then
4712  ifincr = mod(ifhr,itprec)
4713  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
4714  else
4715  ifincr = 0
4716  endif
4717 !mp
4718  id(18) = 0
4719  id(19) = ifhr
4720  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4721  id(20) = 3
4722  IF (ifincr==0) THEN
4723  id(18) = ifhr-itprec
4724  ELSE
4725  id(18) = ifhr-ifincr
4726  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4727  ENDIF
4728  grid1=spval
4729 !$omp parallel do private(i,j)
4730  DO j=jsta,jend
4731  DO i=ista,iend
4732  if(avgprec(i,j)/=spval) grid1(i,j) = domip(i,j)
4733  ENDDO
4734  ENDDO
4735  if(grib=='grib2') then
4736  cfld=cfld+1
4737  fld_info(cfld)%ifld=iavblfld(iget(556))
4738  if(itprec==0) then
4739  fld_info(cfld)%ntrange=0
4740  else
4741  fld_info(cfld)%ntrange=1
4742  endif
4743  fld_info(cfld)%tinvstat=ifhr-id(18)
4744 
4745 !$omp parallel do private(i,j,ii,jj)
4746  do j=1,jend-jsta+1
4747  jj = jsta+j-1
4748  do i=1,iend-ista+1
4749  ii = ista+i-1
4750  datapd(i,j,cfld) = grid1(ii,jj)
4751  enddo
4752  enddo
4753  endif
4754 ! FREEZING RAIN.
4755  id(8) = 141
4756 
4757  itprec = nint(tprec)
4758 !mp
4759  if (itprec /= 0) then
4760  ifincr = mod(ifhr,itprec)
4761  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
4762  else
4763  ifincr = 0
4764  endif
4765 !mp
4766  id(18) = 0
4767  id(19) = ifhr
4768  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4769  id(20) = 3
4770  IF (ifincr==0) THEN
4771  id(18) = ifhr-itprec
4772  ELSE
4773  id(18) = ifhr-ifincr
4774  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4775  ENDIF
4776  grid1=spval
4777 !$omp parallel do private(i,j)
4778  DO j=jsta,jend
4779  DO i=ista,iend
4780 ! if (DOMZR(I,J) == 1) THEN
4781 ! PSFC(I,J)=PINT(I,J,NINT(LMH(I,J))+1)
4782 ! print *, 'aha ', I, J, PSFC(I,J)
4783 ! print *, FREEZR(I,J,1), FREEZR(I,J,2),
4784 ! * FREEZR(I,J,3), FREEZR(I,J,4), FREEZR(I,J,5)
4785 ! endif
4786  if(avgprec(i,j)/=spval) grid1(i,j) = domzr(i,j)
4787  ENDDO
4788  ENDDO
4789  if(grib=='grib2') then
4790  cfld=cfld+1
4791  fld_info(cfld)%ifld=iavblfld(iget(557))
4792  if(itprec==0) then
4793  fld_info(cfld)%ntrange=0
4794  else
4795  fld_info(cfld)%ntrange=1
4796  endif
4797  fld_info(cfld)%tinvstat=ifhr-id(18)
4798 
4799 !$omp parallel do private(i,j,ii,jj)
4800  do j=1,jend-jsta+1
4801  jj = jsta+j-1
4802  do i=1,iend-ista+1
4803  ii = ista+i-1
4804  datapd(i,j,cfld) = grid1(ii,jj)
4805  enddo
4806  enddo
4807  endif
4808 ! RAIN.
4809  id(8) = 140
4810 
4811  itprec = nint(tprec)
4812 !mp
4813  if (itprec /= 0) then
4814  ifincr = mod(ifhr,itprec)
4815  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
4816  else
4817  ifincr = 0
4818  endif
4819 !mp:w
4820 
4821  id(18) = 0
4822  id(19) = ifhr
4823  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4824  id(20) = 3
4825  IF (ifincr==0) THEN
4826  id(18) = ifhr-itprec
4827  ELSE
4828  id(18) = ifhr-ifincr
4829  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4830  ENDIF
4831  grid1=spval
4832 !$omp parallel do private(i,j)
4833  DO j=jsta,jend
4834  DO i=ista,iend
4835  if(avgprec(i,j)/=spval) grid1(i,j) = domr(i,j)
4836  ENDDO
4837  ENDDO
4838  if(grib=='grib2') then
4839  cfld=cfld+1
4840  fld_info(cfld)%ifld=iavblfld(iget(317))
4841  if(itprec==0) then
4842  fld_info(cfld)%ntrange=0
4843  else
4844  fld_info(cfld)%ntrange=1
4845  endif
4846  fld_info(cfld)%tinvstat=ifhr-id(18)
4847 
4848 !$omp parallel do private(i,j,ii,jj)
4849  do j=1,jend-jsta+1
4850  jj = jsta+j-1
4851  do i=1,iend-ista+1
4852  ii = ista+i-1
4853  datapd(i,j,cfld) = grid1(ii,jj)
4854  enddo
4855  enddo
4856  endif
4857 
4858  ENDIF
4859 
4860  if (allocated(rain)) deallocate(rain)
4861  if (allocated(snow)) deallocate(snow)
4862  if (allocated(sleet)) deallocate(sleet)
4863  if (allocated(freezr)) deallocate(freezr)
4864 
4865 ! GSD PRECIPITATION TYPE
4866  IF (iget(407)>0 .or. iget(559)>0 .or. &
4867  iget(560)>0 .or. iget(561)>0) THEN
4868 
4869  if (.not. allocated(domr)) allocate(domr(ista:iend,jsta:jend))
4870  if (.not. allocated(doms)) allocate(doms(ista:iend,jsta:jend))
4871  if (.not. allocated(domzr)) allocate(domzr(ista:iend,jsta:jend))
4872  if (.not. allocated(domip)) allocate(domip(ista:iend,jsta:jend))
4873 
4874 !$omp parallel do private(i,j)
4875  DO j=jsta,jend
4876  DO i=ista,iend
4877  doms(i,j) = 0. !-- snow
4878  domr(i,j) = 0. !-- rain
4879  domzr(i,j) = 0. !-- freezing rain
4880  domip(i,j) = 0. !-- ice pellets
4881  ENDDO
4882  ENDDO
4883 
4884  DO j=jsta,jend
4885  DO i=ista,iend
4886 !-- TOTPRCP is total 1-hour accumulated precipitation in [m]
4887  totprcp = (rainc_bucket(i,j) + rainnc_bucket(i,j))*1.e-3
4888  snowratio = 0.0
4889  if(graup_bucket(i,j)*1.e-3 > totprcp)then
4890  print *,'WARNING - Graupel is higher that total precip at point',i,j
4891  print *,'totprcp,graup_bucket(i,j),snow_bucket(i,j),rainnc_bucket',&
4892  totprcp,graup_bucket(i,j),snow_bucket(i,j),rainnc_bucket(i,j)
4893  endif
4894 
4895 ! ---------------------------------------------------------------
4896 ! Minimum 1h precipitation to even consider p-type specification
4897 ! (0.0001 mm in 1h, very light precipitation)
4898 ! ---------------------------------------------------------------
4899  if (totprcp-graup_bucket(i,j)*1.e-3 > 0.0000001) &
4900 ! snowratio = snow_bucket(i,j)*1.e-3/totprcp ! orig
4901 !14aug15 - change from Stan and Trevor
4902 ! ---------------------------------------------------------------
4903 ! Snow-to-total ratio to be used below
4904 ! ---------------------------------------------------------------
4905  snowratio = snow_bucket(i,j)*1.e-3 / (totprcp-graup_bucket(i,j)*1.e-3)
4906 
4907 ! snowratio = SR(i,j)
4908 !-- 2-m temperature
4909  t2 = tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
4910 ! ---------------------------------------------------------------
4911 !--snow (or rain if T2m > 3 C)
4912 ! ---------------------------------------------------------------
4913 !-- SNOW is time step non-convective snow [m]
4914 ! -- based on either instantaneous snowfall or 1h snowfall and
4915 ! snowratio
4916  if( (snownc(i,j)/dt > 0.2e-9 .and. snowratio>=0.25) &
4917  .or. &
4918  (totprcp>0.00001.and.snowratio>=0.25)) then
4919  doms(i,j) = 1.
4920  if (t2>=276.15) then
4921 ! switch snow to rain if 2m temp > 3 deg
4922  domr(i,j) = 1.
4923  doms(i,j) = 0.
4924  end if
4925  end if
4926 
4927 ! ---------------------------------------------------------------
4928 !-- rain/freezing rain
4929 ! ---------------------------------------------------------------
4930 !-- compute RAIN [m/s] from total convective and non-convective precipitation
4931  rainl = (1. - sr(i,j))*prec(i,j)/dt
4932 !-- in RUC RAIN is in cm/h and the limit is 1.e-3,
4933 !-- converted to m/s will be 2.8e-9
4934  if((rainl > 2.8e-9 .and. snowratio<0.60) .or. &
4935  (totprcp>0.00001 .and. snowratio<0.60)) then
4936 
4937  if (t2>=273.15) then
4938 !--rain
4939  domr(i,j) = 1.
4940 ! else if (tmax(i,j)>273.15) then
4941 !14aug15 - stan
4942  else
4943 !-- freezing rain
4944  domzr(i,j) = 1.
4945  endif
4946  endif
4947 
4948 ! ---------------------------------------------------------------
4949 !-- graupel/ice pellets vs. snow or rain
4950 ! ---------------------------------------------------------------
4951 !-- GRAUPEL is time step non-convective graupel in [m]
4952  if(graupelnc(i,j)/dt > 1.e-9) then
4953  if (t2<=276.15) then
4954 ! This T2m test excludes convectively based hail
4955 ! from cold-season ice pellets.
4956 
4957 ! check for max rain mixing ratio
4958 ! if it's > 0.05 g/kg, => ice pellets
4959  if (qrmax(i,j)>0.000005) then
4960  if(graupelnc(i,j) > 0.5*snownc(i,j)) then
4961 ! if (instantaneous graupel fall rate > 0.5*
4962 ! instantaneous snow fall rate, ....
4963 !-- diagnose ice pellets
4964  domip(i,j) = 1.
4965 
4966 ! -- If graupel is greater than rain,
4967 ! report graupel only
4968 ! in RUC --> if (3.6E5*gex2(i,j,8)> gex2(i,j,6)) then
4969  if ((graupelnc(i,j)/dt) > rainl) then
4970  domip(i,j) = 1.
4971  domzr(i,j) = 0.
4972  domr(i,j) = 0.
4973 ! -- If rain is greater than 4x graupel,
4974 ! report rain only
4975 ! in RUC --> else if (gex2(i,j,6)>4.*3.6E5*gex2(i,j,8)) then
4976  else if (rainl > (4.*graupelnc(i,j)/dt)) then
4977  domip(i,j) = 0.
4978  end if
4979 
4980  else ! instantaneous graupel fall rate <
4981  ! 0.5 * instantaneous snow fall rate
4982 ! snow -- ensure snow is diagnosed (no ice pellets)
4983  doms(i,j)=1.
4984  end if
4985  else ! if qrmax is not > 0.00005
4986 ! snow
4987  doms(i,j)=1.
4988  end if
4989 
4990  else ! if t2 >= 3 deg C
4991 ! rain
4992  domr(i,j) = 1.
4993  end if ! End of t2 if/then loop
4994 
4995  end if ! End of GRAUPELNC if/then loop
4996 
4997  ENDDO
4998  ENDDO
4999 
5000 
5001  write (6,*)' Snow/rain ratio'
5002  write (6,*)' max/min 1h-SNOWFALL in [cm]', &
5003  maxval(snow_bucket)*0.1,minval(snow_bucket)*0.1
5004 
5005  DO j=jsta,jend
5006  DO i=ista,iend
5007  do icat=1,10
5008  if (snow_bucket(i,j)*0.1<0.1*float(icat).and. &
5009  snow_bucket(i,j)*0.1>0.1*float(icat-1)) then
5010  cnt_snowratio(icat)=cnt_snowratio(icat)+1
5011  end if
5012  end do
5013  end do
5014  end do
5015 
5016  write (6,*) 'Snow ratio point counts'
5017  do icat=1,10
5018  write (6,*) icat, cnt_snowratio(icat)
5019  end do
5020 
5021  icnt_snow_rain_mixed = 0
5022  DO j=jsta,jend
5023  DO i=ista,iend
5024  if (domr(i,j)==1 .and. doms(i,j)==1) then
5025  icnt_snow_rain_mixed = icnt_snow_rain_mixed + 1
5026  endif
5027  end do
5028  end do
5029 
5030  write (6,*) 'No. of mixed snow/rain p-type diagnosed=', &
5031  icnt_snow_rain_mixed
5032 
5033 
5034 ! SNOW.
5035 !$omp parallel do private(i,j)
5036  DO j=jsta,jend
5037  DO i=ista,iend
5038  grid1(i,j)=doms(i,j)
5039  ENDDO
5040  ENDDO
5041  if(grib=='grib2') then
5042  cfld=cfld+1
5043  fld_info(cfld)%ifld=iavblfld(iget(559))
5044 !$omp parallel do private(i,j,ii,jj)
5045  do j=1,jend-jsta+1
5046  jj = jsta+j-1
5047  do i=1,iend-ista+1
5048  ii = ista+i-1
5049  datapd(i,j,cfld) = grid1(ii,jj)
5050  enddo
5051  enddo
5052  endif
5053 ! ICE PELLETS.
5054 !$omp parallel do private(i,j)
5055  DO j=jsta,jend
5056  DO i=ista,iend
5057  grid1(i,j) = domip(i,j)
5058 ! if (DOMIP(I,J) == 1) THEN
5059 ! print *, 'ICE PELLETS at I,J ', I, J
5060 ! endif
5061  ENDDO
5062  ENDDO
5063  if(grib=='grib2') then
5064  cfld=cfld+1
5065  fld_info(cfld)%ifld=iavblfld(iget(560))
5066 !$omp parallel do private(i,j,ii,jj)
5067  do j=1,jend-jsta+1
5068  jj = jsta+j-1
5069  do i=1,iend-ista+1
5070  ii = ista+i-1
5071  datapd(i,j,cfld) = grid1(ii,jj)
5072  enddo
5073  enddo
5074  endif
5075 ! FREEZING RAIN.
5076 !$omp parallel do private(i,j)
5077  DO j=jsta,jend
5078  DO i=ista,iend
5079 ! if (DOMZR(I,J) == 1) THEN
5080 ! PSFC(I,J)=PINT(I,J,NINT(LMH(I,J))+1)
5081 ! print *, 'FREEZING RAIN AT I,J ', I, J, PSFC(I,J)
5082 ! endif
5083  grid1(i,j) = domzr(i,j)
5084  ENDDO
5085  ENDDO
5086  if(grib=='grib2') then
5087  cfld=cfld+1
5088  fld_info(cfld)%ifld=iavblfld(iget(561))
5089 !$omp parallel do private(i,j,ii,jj)
5090  do j=1,jend-jsta+1
5091  jj = jsta+j-1
5092  do i=1,iend-ista+1
5093  ii = ista+i-1
5094  datapd(i,j,cfld) = grid1(ii,jj)
5095  enddo
5096  enddo
5097  endif
5098 ! RAIN.
5099 !$omp parallel do private(i,j)
5100  DO j=jsta,jend
5101  DO i=ista,iend
5102  grid1(i,j) = domr(i,j)
5103  ENDDO
5104  ENDDO
5105  if(grib=='grib2') then
5106  cfld=cfld+1
5107  fld_info(cfld)%ifld=iavblfld(iget(407))
5108 !$omp parallel do private(i,j,ii,jj)
5109  do j=1,jend-jsta+1
5110  jj = jsta+j-1
5111  do i=1,iend-ista+1
5112  ii = ista+i-1
5113  datapd(i,j,cfld) = grid1(ii,jj)
5114  enddo
5115  enddo
5116  endif
5117 
5118  ENDIF
5119 !
5120  if (allocated(psfc)) deallocate(psfc)
5121  if (allocated(domr)) deallocate(domr)
5122  if (allocated(doms)) deallocate(doms)
5123  if (allocated(domzr)) deallocate(domzr)
5124  if (allocated(domip)) deallocate(domip)
5125 !
5126 !
5127 !*** BLOCK 5. SURFACE EXCHANGE FIELDS.
5128 !
5129 ! TIME AVERAGED SURFACE LATENT HEAT FLUX.
5130  IF (iget(042)>0) THEN
5131  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5132  modelname=='RAPR')THEN
5133  grid1=spval
5134  id(1:25)=0
5135  ELSE
5136  IF(asrfc>0.)THEN
5137  rrnum=1./asrfc
5138  ELSE
5139  rrnum=0.
5140  ENDIF
5141  DO j=jsta,jend
5142  DO i=ista,iend
5143  IF(sfclhx(i,j)/=spval)THEN
5144  grid1(i,j)=-1.*sfclhx(i,j)*rrnum !change the sign to conform with Grib
5145  ELSE
5146  grid1(i,j)=sfclhx(i,j)
5147  END IF
5148  ENDDO
5149  ENDDO
5150  id(1:25) = 0
5151  itsrfc = nint(tsrfc)
5152  IF(itsrfc /= 0) then
5153  ifincr = mod(ifhr,itsrfc)
5154  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5155  ELSE
5156  ifincr = 0
5157  endif
5158  id(19) = ifhr
5159  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5160  id(20) = 3
5161  IF (ifincr==0) THEN
5162  id(18) = ifhr-itsrfc
5163  ELSE
5164  id(18) = ifhr-ifincr
5165  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5166  ENDIF
5167  IF (id(18)<0) id(18) = 0
5168  if(grib=='grib2') then
5169  cfld=cfld+1
5170  fld_info(cfld)%ifld=iavblfld(iget(042))
5171  if(itsrfc>0) then
5172  fld_info(cfld)%ntrange=1
5173  else
5174  fld_info(cfld)%ntrange=0
5175  endif
5176  fld_info(cfld)%tinvstat=ifhr-id(18)
5177  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5178  endif
5179  END IF
5180  ENDIF
5181 !
5182 ! TIME AVERAGED SURFACE SENSIBLE HEAT FLUX.
5183  IF (iget(043)>0) THEN
5184  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5185  modelname=='RAPR')THEN
5186  grid1=spval
5187  id(1:25)=0
5188  ELSE
5189  IF(asrfc>0.)THEN
5190  rrnum=1./asrfc
5191  ELSE
5192  rrnum=0.
5193  ENDIF
5194  DO j=jsta,jend
5195  DO i=ista,iend
5196  IF(sfcshx(i,j)/=spval)THEN
5197  grid1(i,j) = -1.* sfcshx(i,j)*rrnum !change the sign to conform with Grib
5198  ELSE
5199  grid1(i,j)=sfcshx(i,j)
5200  END IF
5201  ENDDO
5202  ENDDO
5203  id(1:25) = 0
5204  itsrfc = nint(tsrfc)
5205  IF(itsrfc /= 0) then
5206  ifincr = mod(ifhr,itsrfc)
5207  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5208  ELSE
5209  ifincr = 0
5210  endif
5211  id(19) = ifhr
5212  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5213  id(20) = 3
5214  IF (ifincr==0) THEN
5215  id(18) = ifhr-itsrfc
5216  ELSE
5217  id(18) = ifhr-ifincr
5218  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5219  ENDIF
5220  IF (id(18)<0) id(18) = 0
5221  END IF
5222  if(grib=='grib2') then
5223  cfld=cfld+1
5224  fld_info(cfld)%ifld=iavblfld(iget(043))
5225  if(itsrfc>0) then
5226  fld_info(cfld)%ntrange=1
5227  else
5228  fld_info(cfld)%ntrange=0
5229  endif
5230  fld_info(cfld)%tinvstat=ifhr-id(18)
5231  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5232  endif
5233  ENDIF
5234 !
5235 ! TIME AVERAGED SUB-SURFACE SENSIBLE HEAT FLUX.
5236  IF (iget(135)>0) THEN
5237  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5238  modelname=='RAPR')THEN
5239  grid1=spval
5240  id(1:25)=0
5241  ELSE
5242  IF(asrfc>0.)THEN
5243  rrnum=1./asrfc
5244  ELSE
5245  rrnum=0.
5246  ENDIF
5247  grid1=spval
5248  DO j=jsta,jend
5249  DO i=ista,iend
5250  if(subshx(i,j)/=spval) grid1(i,j) = subshx(i,j)*rrnum
5251  ENDDO
5252  ENDDO
5253  id(1:25) = 0
5254  itsrfc = nint(tsrfc)
5255  IF(itsrfc /= 0) then
5256  ifincr = mod(ifhr,itsrfc)
5257  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5258  ELSE
5259  ifincr = 0
5260  endif
5261  id(19) = ifhr
5262  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5263  id(20) = 3
5264  IF (ifincr==0) THEN
5265  id(18) = ifhr-itsrfc
5266  ELSE
5267  id(18) = ifhr-ifincr
5268  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5269  ENDIF
5270  IF (id(18)<0) id(18) = 0
5271  END IF
5272  if(grib=='grib2') then
5273  cfld=cfld+1
5274  fld_info(cfld)%ifld=iavblfld(iget(135))
5275  if(itsrfc>0) then
5276  fld_info(cfld)%ntrange=1
5277  else
5278  fld_info(cfld)%ntrange=0
5279  endif
5280  fld_info(cfld)%tinvstat=ifhr-id(18)
5281  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5282  endif
5283  ENDIF
5284 !
5285 ! TIME AVERAGED SNOW PHASE CHANGE HEAT FLUX.
5286  IF (iget(136)>0) THEN
5287  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5288  modelname=='RAPR')THEN
5289  grid1=spval
5290  id(1:25)=0
5291  ELSE
5292  IF(asrfc>0.)THEN
5293  rrnum=1./asrfc
5294  ELSE
5295  rrnum=0.
5296  ENDIF
5297  grid1=spval
5298  DO j=jsta,jend
5299  DO i=ista,iend
5300  if(snopcx(i,j)/=spval) grid1(i,j) = snopcx(i,j)*rrnum
5301  ENDDO
5302  ENDDO
5303  id(1:25) = 0
5304  itsrfc = nint(tsrfc)
5305  IF(itsrfc /= 0) then
5306  ifincr = mod(ifhr,itsrfc)
5307  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5308  ELSE
5309  ifincr = 0
5310  endif
5311  id(19) = ifhr
5312  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5313  id(20) = 3
5314  IF (ifincr==0) THEN
5315  id(18) = ifhr-itsrfc
5316  ELSE
5317  id(18) = ifhr-ifincr
5318  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5319  ENDIF
5320  IF (id(18)<0) id(18) = 0
5321  END IF
5322  if(grib=='grib2') then
5323  cfld=cfld+1
5324  fld_info(cfld)%ifld=iavblfld(iget(136))
5325  if(itsrfc>0) then
5326  fld_info(cfld)%ntrange=1
5327  else
5328  fld_info(cfld)%ntrange=0
5329  endif
5330  fld_info(cfld)%tinvstat=ifhr-id(18)
5331  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5332  endif
5333  ENDIF
5334 !
5335 ! TIME AVERAGED SURFACE MOMENTUM FLUX.
5336  IF (iget(046)>0) THEN
5337  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5338  modelname=='RAPR')THEN
5339  grid1=spval
5340  id(1:25)=0
5341  ELSE
5342  IF(asrfc>0.)THEN
5343  rrnum=1./asrfc
5344  ELSE
5345  rrnum=0.
5346  ENDIF
5347  DO j=jsta,jend
5348  DO i=ista,iend
5349  IF(sfcuvx(i,j)/=spval)THEN
5350  grid1(i,j) = sfcuvx(i,j)*rrnum
5351  ELSE
5352  grid1(i,j) = sfcuvx(i,j)
5353  END IF
5354  ENDDO
5355  ENDDO
5356  id(1:25) = 0
5357  itsrfc = nint(tsrfc)
5358  IF(itsrfc /= 0) then
5359  ifincr = mod(ifhr,itsrfc)
5360  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5361  ELSE
5362  ifincr = 0
5363  endif
5364  id(19) = ifhr
5365  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5366  id(20) = 3
5367  IF (ifincr==0) THEN
5368  id(18) = ifhr-itsrfc
5369  ELSE
5370  id(18) = ifhr-ifincr
5371  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5372  ENDIF
5373  IF (id(18)<0) id(18) = 0
5374  END IF
5375  if(grib=='grib2') then
5376  cfld=cfld+1
5377  fld_info(cfld)%ifld=iavblfld(iget(046))
5378  if(itsrfc>0) then
5379  fld_info(cfld)%ntrange=1
5380  else
5381  fld_info(cfld)%ntrange=0
5382  endif
5383  fld_info(cfld)%tinvstat=ifhr-id(18)
5384  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5385  endif
5386  ENDIF
5387 !
5388 ! TIME AVERAGED SURFACE ZONAL MOMENTUM FLUX.
5389  IF (iget(269)>0) THEN
5390  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5391  modelname=='RAPR')THEN
5392  grid1=spval
5393  id(1:25)=0
5394  ELSE
5395  IF(asrfc>0.)THEN
5396  rrnum=1./asrfc
5397  ELSE
5398  rrnum=0.
5399  ENDIF
5400  grid1=spval
5401  DO j=jsta,jend
5402  DO i=ista,iend
5403  if(sfcux(i,j)/=spval) grid1(i,j) = sfcux(i,j)*rrnum
5404  ENDDO
5405  ENDDO
5406  id(1:25) = 0
5407  itsrfc = nint(tsrfc)
5408  IF(itsrfc /= 0) then
5409  ifincr = mod(ifhr,itsrfc)
5410  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5411  ELSE
5412  ifincr = 0
5413  endif
5414  id(19) = ifhr
5415  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5416  id(20) = 3
5417  IF (ifincr==0) THEN
5418  id(18) = ifhr-itsrfc
5419  ELSE
5420  id(18) = ifhr-ifincr
5421  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5422  ENDIF
5423  IF (id(18)<0) id(18) = 0
5424  END IF
5425  if(grib=='grib2') then
5426  cfld=cfld+1
5427  fld_info(cfld)%ifld=iavblfld(iget(269))
5428  if(itsrfc>0) then
5429  fld_info(cfld)%ntrange=1
5430  else
5431  fld_info(cfld)%ntrange=0
5432  endif
5433  fld_info(cfld)%tinvstat=ifhr-id(18)
5434  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5435  endif
5436  ENDIF
5437 !
5438 ! TIME AVERAGED SURFACE MOMENTUM FLUX.
5439  IF (iget(270)>0) THEN
5440  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. &
5441  modelname=='RAPR')THEN
5442  grid1=spval
5443  id(1:25)=0
5444  ELSE
5445  IF(asrfc>0.)THEN
5446  rrnum=1./asrfc
5447  ELSE
5448  rrnum=0.
5449  ENDIF
5450  grid1=spval
5451  DO j=jsta,jend
5452  DO i=ista,iend
5453  if(sfcvx(i,j)/=spval) grid1(i,j) = sfcvx(i,j)*rrnum
5454  ENDDO
5455  ENDDO
5456  id(1:25) = 0
5457  itsrfc = nint(tsrfc)
5458  IF(itsrfc /= 0) then
5459  ifincr = mod(ifhr,itsrfc)
5460  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5461  ELSE
5462  ifincr = 0
5463  endif
5464  id(19) = ifhr
5465  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5466  id(20) = 3
5467  IF (ifincr==0) THEN
5468  id(18) = ifhr-itsrfc
5469  ELSE
5470  id(18) = ifhr-ifincr
5471  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5472  ENDIF
5473  IF (id(18)<0) id(18) = 0
5474  END IF
5475  if(grib=='grib2') then
5476  cfld=cfld+1
5477  fld_info(cfld)%ifld=iavblfld(iget(270))
5478  if(itsrfc>0) then
5479  fld_info(cfld)%ntrange=1
5480  else
5481  fld_info(cfld)%ntrange=0
5482  endif
5483  fld_info(cfld)%tinvstat=ifhr-id(18)
5484  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5485  endif
5486  ENDIF
5487 !
5488 ! ACCUMULATED SURFACE EVAPORATION
5489  IF (iget(047)>0) THEN
5490  grid1=spval
5491  DO j=jsta,jend
5492  DO i=ista,iend
5493  if(sfcevp(i,j)/=spval) grid1(i,j) = sfcevp(i,j)*1000.
5494  ENDDO
5495  ENDDO
5496  id(1:25) = 0
5497  itprec = nint(tprec)
5498 !mp
5499  if (itprec /= 0) then
5500  ifincr = mod(ifhr,itprec)
5501  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
5502  else
5503  ifincr = 0
5504  endif
5505 !mp
5506  id(18) = 0
5507  id(19) = ifhr
5508  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5509  id(20) = 4
5510  IF (ifincr==0) THEN
5511  id(18) = ifhr-itprec
5512  ELSE
5513  id(18) = ifhr-ifincr
5514  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5515  ENDIF
5516  IF (id(18)<0) id(18) = 0
5517  if(grib=='grib2') then
5518  cfld=cfld+1
5519  fld_info(cfld)%ifld=iavblfld(iget(047))
5520  if(itprec>0) then
5521  fld_info(cfld)%ntrange=1
5522  else
5523  fld_info(cfld)%ntrange=0
5524  endif
5525  fld_info(cfld)%tinvstat=ifhr-id(18)
5526  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5527 
5528  endif
5529  ENDIF
5530 !
5531 ! ACCUMULATED POTENTIAL EVAPORATION
5532  IF (iget(137)>0) THEN
5533  grid1=spval
5534  DO j=jsta,jend
5535  DO i=ista,iend
5536  if(potevp(i,j)/=spval) grid1(i,j) = potevp(i,j)*1000.
5537  ENDDO
5538  ENDDO
5539  id(1:25) = 0
5540  itprec = nint(tprec)
5541 !mp
5542  if (itprec /= 0) then
5543  ifincr = mod(ifhr,itprec)
5544  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
5545  else
5546  ifincr = 0
5547  endif
5548 !mp
5549  id(18) = 0
5550  id(19) = ifhr
5551  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5552  id(20) = 4
5553  IF (ifincr==0) THEN
5554  id(18) = ifhr-itprec
5555  ELSE
5556  id(18) = ifhr-ifincr
5557  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5558  ENDIF
5559  IF (id(18)<0) id(18) = 0
5560  if(grib=='grib2') then
5561  cfld=cfld+1
5562  fld_info(cfld)%ifld=iavblfld(iget(137))
5563  if(itprec>0) then
5564  fld_info(cfld)%ntrange=1
5565  else
5566  fld_info(cfld)%ntrange=0
5567  endif
5568  fld_info(cfld)%tinvstat=ifhr-id(18)
5569  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5570  endif
5571  ENDIF
5572 !
5573 ! ROUGHNESS LENGTH.
5574  IF (iget(044)>0) THEN
5575  DO j=jsta,jend
5576  DO i=ista,iend
5577  grid1(i,j) = z0(i,j)
5578  ENDDO
5579  ENDDO
5580  if(grib=='grib2') then
5581  cfld=cfld+1
5582  fld_info(cfld)%ifld=iavblfld(iget(044))
5583  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5584  endif
5585  ENDIF
5586 !
5587 ! FRICTION VELOCITY.
5588  IF (iget(045)>0) THEN
5589  DO j=jsta,jend
5590  DO i=ista,iend
5591  grid1(i,j) = ustar(i,j)
5592  ENDDO
5593  ENDDO
5594  if(grib=='grib2') then
5595  cfld=cfld+1
5596  fld_info(cfld)%ifld=iavblfld(iget(045))
5597  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5598  endif
5599  ENDIF
5600 !
5601 ! SURFACE DRAG COEFFICIENT.
5602 ! dong add missing value for cd
5603  IF (iget(132)>0) THEN
5604  grid1=spval
5605  CALL caldrg(egrid1(ista_2l:iend_2u,jsta_2l:jend_2u))
5606  DO j=jsta,jend
5607  DO i=ista,iend
5608  IF(ustar(i,j) < spval) grid1(i,j)=egrid1(i,j)
5609  ENDDO
5610  ENDDO
5611  if(grib=='grib2') then
5612  cfld=cfld+1
5613  fld_info(cfld)%ifld=iavblfld(iget(132))
5614  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5615  endif
5616  ENDIF
5617 
5618  write_cd: IF(iget(922)>0) THEN
5619  DO j=jsta,jend
5620  DO i=ista,iend
5621  grid1(i,j)=cd10(i,j)
5622  ENDDO
5623  ENDDO
5624  if(grib=='grib2') then
5625  cfld=cfld+1
5626  fld_info(cfld)%ifld=iavblfld(iget(922))
5627  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5628  endif
5629  ENDIF write_cd
5630  write_ch: IF(iget(923)>0) THEN
5631  DO j=jsta,jend
5632  DO i=ista,iend
5633  grid1(i,j)=ch10(i,j)
5634  ENDDO
5635  ENDDO
5636  if(grib=='grib2') then
5637  cfld=cfld+1
5638  fld_info(cfld)%ifld=iavblfld(iget(923))
5639  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5640  endif
5641  ENDIF write_ch
5642 !
5643 ! MODEL OUTPUT SURFACE U AND/OR V COMPONENT WIND STRESS
5644  IF ( (iget(900)>0) .OR. (iget(901)>0) ) THEN
5645 !
5646 ! MODEL OUTPUT SURFACE U COMPONENT WIND STRESS.
5647  IF (iget(900)>0) THEN
5648  DO j=jsta,jend
5649  DO i=ista,iend
5650  grid1(i,j)=mdltaux(i,j)
5651  ENDDO
5652  ENDDO
5653  if(grib=='grib2') then
5654  cfld=cfld+1
5655  fld_info(cfld)%ifld=iavblfld(iget(900))
5656  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5657  endif
5658 
5659  ENDIF
5660 !
5661 ! MODEL OUTPUT SURFACE V COMPONENT WIND STRESS
5662  IF (iget(901)>0) THEN
5663  DO j=jsta,jend
5664  DO i=ista,iend
5665  grid1(i,j)=mdltauy(i,j)
5666  ENDDO
5667  ENDDO
5668  if(grib=='grib2') then
5669  cfld=cfld+1
5670  fld_info(cfld)%ifld=iavblfld(iget(901))
5671  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5672  endif
5673  ENDIF
5674  ENDIF
5675 !
5676 ! SURFACE U AND/OR V COMPONENT WIND STRESS
5677  IF ( (iget(133)>0) .OR. (iget(134)>0) ) THEN
5678 ! dong add missing value
5679  grid1 = spval
5680  IF(modelname /= 'FV3R') &
5681  CALL caltau(egrid1(ista:iend,jsta:jend),egrid2(ista:iend,jsta:jend))
5682 !
5683 ! SURFACE U COMPONENT WIND STRESS.
5684 ! dong for FV3, directly use model output
5685  IF (iget(133)>0) THEN
5686  DO j=jsta,jend
5687  DO i=ista,iend
5688  IF(modelname == 'FV3R') THEN
5689  grid1(i,j)=sfcuxi(i,j)
5690  ELSE
5691  grid1(i,j)=egrid1(i,j)
5692  ENDIF
5693  ENDDO
5694  ENDDO
5695 !
5696  if(grib=='grib2') then
5697  cfld=cfld+1
5698  fld_info(cfld)%ifld=iavblfld(iget(133))
5699  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5700  endif
5701  ENDIF
5702 !
5703 ! SURFACE V COMPONENT WIND STRESS
5704  IF (iget(134)>0) THEN
5705  DO j=jsta,jend
5706  DO i=ista,iend
5707  IF(modelname == 'FV3R') THEN
5708  grid1(i,j)=sfcvxi(i,j)
5709  ELSE
5710  grid1(i,j)=egrid2(i,j)
5711  END IF
5712  ENDDO
5713  ENDDO
5714  if(grib=='grib2') then
5715  cfld=cfld+1
5716  fld_info(cfld)%ifld=iavblfld(iget(134))
5717  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5718  endif
5719  ENDIF
5720  ENDIF
5721 !
5722 ! GRAVITY U AND/OR V COMPONENT STRESS
5723  IF ( (iget(315)>0) .OR. (iget(316)>0) ) THEN
5724 !
5725 ! GRAVITY U COMPONENT WIND STRESS.
5726  IF (iget(315)>0) THEN
5727  DO j=jsta,jend
5728  DO i=ista,iend
5729  grid1(i,j) = gtaux(i,j)
5730  ENDDO
5731  ENDDO
5732  id(1:25) = 0
5733  itsrfc = nint(tsrfc)
5734  IF(itsrfc /= 0) then
5735  ifincr = mod(ifhr,itsrfc)
5736  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5737  ELSE
5738  ifincr = 0
5739  endif
5740  id(19) = ifhr
5741  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5742  id(20) = 3
5743  IF (ifincr==0) THEN
5744  id(18) = ifhr-itsrfc
5745  ELSE
5746  id(18) = ifhr-ifincr
5747  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5748  ENDIF
5749  IF (id(18)<0) id(18) = 0
5750  if(grib=='grib2') then
5751  cfld=cfld+1
5752  fld_info(cfld)%ifld=iavblfld(iget(315))
5753  if(itsrfc==0) then
5754  fld_info(cfld)%ntrange=0
5755  else
5756  fld_info(cfld)%ntrange=1
5757  endif
5758  fld_info(cfld)%tinvstat=ifhr-id(18)
5759  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5760  endif
5761  ENDIF
5762 !
5763 ! SURFACE V COMPONENT WIND STRESS
5764  IF (iget(316)>0) THEN
5765  DO j=jsta,jend
5766  DO i=ista,iend
5767  grid1(i,j)=gtauy(i,j)
5768  ENDDO
5769  ENDDO
5770  id(1:25) = 0
5771  itsrfc = nint(tsrfc)
5772  IF(itsrfc /= 0) then
5773  ifincr = mod(ifhr,itsrfc)
5774  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
5775  ELSE
5776  ifincr = 0
5777  endif
5778  id(19) = ifhr
5779  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
5780  id(20) = 3
5781  IF (ifincr==0) THEN
5782  id(18) = ifhr-itsrfc
5783  ELSE
5784  id(18) = ifhr-ifincr
5785  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
5786  ENDIF
5787  IF (id(18)<0) id(18) = 0
5788  if(grib=='grib2') then
5789  cfld=cfld+1
5790  fld_info(cfld)%ifld=iavblfld(iget(316))
5791  if(itsrfc==0) then
5792  fld_info(cfld)%ntrange=0
5793  else
5794  fld_info(cfld)%ntrange=1
5795  endif
5796  fld_info(cfld)%tinvstat=ifhr-id(18)
5797  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5798  endif
5799  ENDIF
5800  ENDIF
5801 !
5802 ! INSTANTANEOUS SENSIBLE HEAT FLUX
5803  IF (iget(154)>0) THEN
5804 ! dong add missing value to shtfl
5805  grid1 = spval
5806  IF(modelname=='NCAR'.OR.modelname=='RSM' .OR. &
5807  modelname=='RAPR')THEN
5808 !4omp parallel do private(i,j)
5809  DO j=jsta,jend
5810  DO i=ista,iend
5811  grid1(i,j) = twbs(i,j)
5812  ENDDO
5813  ENDDO
5814  ELSE
5815 !4omp parallel do private(i,j)
5816  DO j=jsta,jend
5817  DO i=ista,iend
5818  IF(twbs(i,j) < spval) grid1(i,j) = -twbs(i,j)
5819  ENDDO
5820  ENDDO
5821  END IF
5822  if(grib=='grib2') then
5823  cfld=cfld+1
5824  fld_info(cfld)%ifld=iavblfld(iget(154))
5825  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5826  endif
5827  ENDIF
5828 !
5829 ! INSTANTANEOUS LATENT HEAT FLUX
5830  IF (iget(155)>0) THEN
5831 ! dong add missing value to lhtfl
5832  grid1 = spval
5833  IF(modelname=='NCAR'.OR.modelname=='RSM' .OR. &
5834  modelname=='RAPR')THEN
5835 !4omp parallel do private(i,j)
5836  DO j=jsta,jend
5837  DO i=ista,iend
5838  grid1(i,j) = qwbs(i,j)
5839  ENDDO
5840  ENDDO
5841  ELSE
5842 !4omp parallel do private(i,j)
5843  DO j=jsta,jend
5844  DO i=ista,iend
5845  IF(qwbs(i,j) < spval) grid1(i,j) = -qwbs(i,j)
5846  ENDDO
5847  ENDDO
5848  END IF
5849  if(grib=='grib2') then
5850  cfld=cfld+1
5851  fld_info(cfld)%ifld=iavblfld(iget(155))
5852  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5853  endif
5854  ENDIF
5855 !
5856 ! SURFACE EXCHANGE COEFF
5857  IF (iget(169)>0) THEN
5858  DO j=jsta,jend
5859  DO i=ista,iend
5860  grid1(i,j)=sfcexc(i,j)
5861  ENDDO
5862  ENDDO
5863  if(grib=='grib2') then
5864  cfld=cfld+1
5865  fld_info(cfld)%ifld=iavblfld(iget(169))
5866  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5867  endif
5868  ENDIF
5869 !
5870 ! GREEN VEG FRACTION
5871  IF (iget(170)>0) THEN
5872  grid1=spval
5873  DO j=jsta,jend
5874  DO i=ista,iend
5875  if(vegfrc(i,j)/=spval) grid1(i,j)=vegfrc(i,j)*100.
5876  ENDDO
5877  ENDDO
5878  if(grib=='grib2') then
5879  cfld=cfld+1
5880  fld_info(cfld)%ifld=iavblfld(iget(170))
5881  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5882  endif
5883  ENDIF
5884 
5885 !
5886 ! MIN GREEN VEG FRACTION
5887  IF (iget(726)>0) THEN
5888  grid1=spval
5889  DO j=jsta,jend
5890  DO i=ista,iend
5891  if(shdmin(i,j)/=spval) grid1(i,j)=shdmin(i,j)*100.
5892  ENDDO
5893  ENDDO
5894  if(grib=='grib2') then
5895  cfld=cfld+1
5896  fld_info(cfld)%ifld=iavblfld(iget(726))
5897  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5898  endif
5899  ENDIF
5900 !
5901 ! MAX GREEN VEG FRACTION
5902  IF (iget(729)>0) THEN
5903  grid1=spval
5904  DO j=jsta,jend
5905  DO i=ista,iend
5906  if(shdmax(i,j)/=spval) grid1(i,j)=shdmax(i,j)*100.
5907  ENDDO
5908  ENDDO
5909  if(grib=='grib2') then
5910  cfld=cfld+1
5911  fld_info(cfld)%ifld=iavblfld(iget(729))
5912  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5913  endif
5914  ENDIF
5915 !
5916 ! LEAF AREA INDEX
5917  IF (modelname == 'NCAR'.OR.modelname=='NMM' .OR. &
5918  modelname == 'FV3R' .OR. modelname=='RAPR')THEN
5919  IF (isf_surface_physics == 2 .OR. modelname=='RAPR') THEN
5920  IF (iget(254)>0) THEN
5921  DO j=jsta,jend
5922  DO i=ista,iend
5923  IF (modelname=='RAPR')THEN
5924  grid1(i,j)=lai(i,j)
5925  ELSE
5926  grid1(i,j) = xlai
5927  ENDIF
5928  ENDDO
5929  ENDDO
5930  if(grib=='grib2') then
5931  cfld=cfld+1
5932  fld_info(cfld)%ifld=iavblfld(iget(254))
5933  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5934  endif
5935  ENDIF
5936  ENDIF
5937  ENDIF
5938 !
5939 ! INSTANTANEOUS GROUND HEAT FLUX
5940  IF (iget(152)>0) THEN
5941  DO j=jsta,jend
5942  DO i=ista,iend
5943  grid1(i,j)=grnflx(i,j)
5944  ENDDO
5945  ENDDO
5946  if(grib=='grib2') then
5947  cfld=cfld+1
5948  fld_info(cfld)%ifld=iavblfld(iget(152))
5949  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5950  endif
5951  ENDIF
5952 ! VEGETATION TYPE
5953  IF (iget(218)>0) THEN
5954  DO j=jsta,jend
5955  DO i=ista,iend
5956  grid1(i,j) = float(ivgtyp(i,j))
5957  ENDDO
5958  ENDDO
5959  if(grib=='grib2') then
5960  cfld=cfld+1
5961  fld_info(cfld)%ifld=iavblfld(iget(218))
5962  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5963  endif
5964  ENDIF
5965 !
5966 ! SOIL TYPE
5967  IF (iget(219)>0) THEN
5968  DO j=jsta,jend
5969  DO i=ista,iend
5970  grid1(i,j) = float(isltyp(i,j))
5971  ENDDO
5972  ENDDO
5973  if(grib=='grib2') then
5974  cfld=cfld+1
5975  fld_info(cfld)%ifld=iavblfld(iget(219))
5976  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5977  endif
5978  ENDIF
5979 ! SLOPE TYPE
5980  IF (iget(223)>0) THEN
5981  DO j=jsta,jend
5982  DO i=ista,iend
5983  grid1(i,j) = float(islope(i,j))
5984  ENDDO
5985  ENDDO
5986  if(grib=='grib2') then
5987  cfld=cfld+1
5988  fld_info(cfld)%ifld=iavblfld(iget(223))
5989  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5990  endif
5991  ENDIF
5992 ! if (me==0)print*,'starting computing canopy conductance'
5993 !
5994 ! CANOPY CONDUCTANCE
5995 ! ONLY OUTPUT NEW LSM FIELDS FOR NMM AND ARW BECAUSE RSM USES OLD SOIL TYPES
5996  IF (modelname == 'NCAR'.OR.modelname=='NMM' .OR. &
5997  modelname == 'FV3R' .OR. modelname=='RAPR')THEN
5998  IF (iget(220)>0 .OR. iget(234)>0 &
5999  & .OR. iget(235)>0 .OR. iget(236)>0 &
6000  & .OR. iget(237)>0 .OR. iget(238)>0 &
6001  & .OR. iget(239)>0 .OR. iget(240)>0 &
6002  & .OR. iget(241)>0 ) THEN
6003  IF (isf_surface_physics == 2) THEN !NSOIL == 4
6004 ! if(me==0)print*,'starting computing canopy conductance'
6005  allocate(rsmin(ista:iend,jsta:jend), smcref(ista:iend,jsta:jend), gc(ista:iend,jsta:jend), &
6006  rcq(ista:iend,jsta:jend), rct(ista:iend,jsta:jend), rcsoil(ista:iend,jsta:jend), rcs(ista:iend,jsta:jend))
6007  DO j=jsta,jend
6008  DO i=ista,iend
6009  IF( (abs(sm(i,j)-0.) < 1.0e-5) .AND. &
6010  & (abs(sice(i,j)-0.) < 1.0e-5) ) THEN
6011  IF(czmean(i,j)>1.e-6) THEN
6012  factrs = czen(i,j)/czmean(i,j)
6013  ELSE
6014  factrs = 0.0
6015  ENDIF
6016 ! SOLAR=HBM2(I,J)*RSWIN(I,J)*FACTRS
6017  llmh = nint(lmh(i,j))
6018  solar = rswin(i,j)*factrs
6019  sfctmp = t(i,j,llmh)
6020  sfcq = q(i,j,llmh)
6021  sfcprs = pint(i,j,llmh+1)
6022 ! IF(IVGTYP(I,J)==0)PRINT*,'IVGTYP ZERO AT ',I,J
6023 ! & ,SM(I,J)
6024  ivg = ivgtyp(i,j)
6025 ! IF(IVGTYP(I,J)==0)IVG=7
6026 ! CALL CANRES(SOLAR,SFCTMP,SFCQ,SFCPRS
6027 ! & ,SMC(I,J,1:NSOIL),GC(I,J),RC,IVG,ISLTYP(I,J))
6028 !
6029  CALL canres(solar,sfctmp,sfcq,sfcprs &
6030  & ,sh2o(i,j,1:nsoil),gc(i,j),rc,ivg,isltyp(i,j) &
6031  & ,rsmin(i,j),nroots(i,j),smcwlt(i,j),smcref(i,j) &
6032  & ,rcs(i,j),rcq(i,j),rct(i,j),rcsoil(i,j),sldpth)
6033  IF(abs(smcwlt(i,j)-0.5)<1.e-5)print*, &
6034  & 'LARGE SMCWLT',i,j,sm(i,j),isltyp(i,j),smcwlt(i,j)
6035  ELSE
6036  gc(i,j) = 0.
6037  rsmin(i,j) = 0.
6038  nroots(i,j) = 0
6039  smcwlt(i,j) = 0.
6040  smcref(i,j) = 0.
6041  rcs(i,j) = 0.
6042  rcq(i,j) = 0.
6043  rct(i,j) = 0.
6044  rcsoil(i,j) = 0.
6045  END IF
6046  ENDDO
6047  ENDDO
6048 
6049  IF (iget(220)>0 )THEN
6050  DO j=jsta,jend
6051  DO i=ista,iend
6052  grid1(i,j) = gc(i,j)
6053  ENDDO
6054  ENDDO
6055  if(grib=='grib2') then
6056  cfld=cfld+1
6057  fld_info(cfld)%ifld=iavblfld(iget(220))
6058  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6059  endif
6060  ENDIF
6061 
6062  IF (iget(234)>0 )THEN
6063  DO j=jsta,jend
6064  DO i=ista,iend
6065  grid1(i,j) = rsmin(i,j)
6066  ENDDO
6067  ENDDO
6068  if(grib=='grib2') then
6069  cfld=cfld+1
6070  fld_info(cfld)%ifld=iavblfld(iget(234))
6071  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6072  endif
6073  ENDIF
6074 
6075  IF (iget(235)>0 )THEN
6076  DO j=jsta,jend
6077  DO i=ista,iend
6078  grid1(i,j) = float(nroots(i,j))
6079  ENDDO
6080  ENDDO
6081  if(grib=='grib2') then
6082  cfld=cfld+1
6083  fld_info(cfld)%ifld=iavblfld(iget(235))
6084  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6085  endif
6086  ENDIF
6087 
6088  IF (iget(236)>0 )THEN
6089  DO j=jsta,jend
6090  DO i=ista,iend
6091  grid1(i,j) = smcwlt(i,j)
6092  ENDDO
6093  ENDDO
6094  if(grib=='grib2') then
6095  cfld=cfld+1
6096  fld_info(cfld)%ifld=iavblfld(iget(236))
6097  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6098  endif
6099  ENDIF
6100 
6101  IF (iget(237)>0 )THEN
6102  DO j=jsta,jend
6103  DO i=ista,iend
6104  grid1(i,j) = smcref(i,j)
6105  ENDDO
6106  ENDDO
6107  if(grib=='grib2') then
6108  cfld=cfld+1
6109  fld_info(cfld)%ifld=iavblfld(iget(237))
6110  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6111  endif
6112  ENDIF
6113 
6114  IF (iget(238)>0 )THEN
6115  DO j=jsta,jend
6116  DO i=ista,iend
6117  grid1(i,j) = rcs(i,j)
6118  ENDDO
6119  ENDDO
6120  if(grib=='grib2') then
6121  cfld=cfld+1
6122  fld_info(cfld)%ifld=iavblfld(iget(238))
6123  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6124  endif
6125  ENDIF
6126 
6127  IF (iget(239)>0 )THEN
6128  DO j=jsta,jend
6129  DO i=ista,iend
6130  grid1(i,j) = rct(i,j)
6131  ENDDO
6132  ENDDO
6133  if(grib=='grib2') then
6134  cfld=cfld+1
6135  fld_info(cfld)%ifld=iavblfld(iget(239))
6136  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6137  endif
6138  ENDIF
6139 
6140  IF (iget(240)>0 )THEN
6141  DO j=jsta,jend
6142  DO i=ista,iend
6143  grid1(i,j) = rcq(i,j)
6144  ENDDO
6145  ENDDO
6146  if(grib=='grib2') then
6147  cfld=cfld+1
6148  fld_info(cfld)%ifld=iavblfld(iget(240))
6149  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6150  endif
6151  ENDIF
6152 
6153  IF (iget(241)>0 )THEN
6154  DO j=jsta,jend
6155  DO i=ista,iend
6156  grid1(i,j) = rcsoil(i,j)
6157  ENDDO
6158  ENDDO
6159  if(grib=='grib2') then
6160  cfld=cfld+1
6161  fld_info(cfld)%ifld=iavblfld(iget(241))
6162  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6163  endif
6164  ENDIF
6165 
6166  if (allocated(rsmin)) deallocate(rsmin)
6167  if (allocated(smcref)) deallocate(smcref)
6168  if (allocated(rcq)) deallocate(rcq)
6169  if (allocated(rct)) deallocate(rct)
6170  if (allocated(rcsoil)) deallocate(rcsoil)
6171  if (allocated(rcs)) deallocate(rcs)
6172  if (allocated(gc)) deallocate(gc)
6173 
6174 
6175  ENDIF
6176  END IF
6177 !GPL added endif here
6178  ENDIF
6179  IF(modelname == 'GFS')THEN
6180 ! Outputting wilting point and field capacity for TIGGE
6181  IF(iget(236)>0)THEN
6182 !$omp parallel do private(i,j)
6183  DO j=jsta,jend
6184  DO i=ista,iend
6185  grid1(i,j) = smcwlt(i,j)
6186 ! IF(isltyp(i,j)/=0)THEN
6187 ! GRID1(I,J) = WLTSMC(isltyp(i,j))
6188 ! ELSE
6189 ! GRID1(I,J) = spval
6190 ! END IF
6191  ENDDO
6192  ENDDO
6193  if(grib=='grib2') then
6194  cfld=cfld+1
6195  fld_info(cfld)%ifld=iavblfld(iget(236))
6196 !$omp parallel do private(i,j,ii,jj)
6197  do j=1,jend-jsta+1
6198  jj = jsta+j-1
6199  do i=1,iend-ista+1
6200  ii = ista+i-1
6201  datapd(i,j,cfld) = grid1(ii,jj)
6202  enddo
6203  enddo
6204  endif
6205  ENDIF
6206 
6207  IF(iget(397)>0)THEN
6208 !$omp parallel do private(i,j)
6209  DO j=jsta,jend
6210  DO i=ista,iend
6211  grid1(i,j) = fieldcapa(i,j)
6212 ! IF(isltyp(i,j)/=0)THEN
6213 ! GRID1(I,J) = REFSMC(isltyp(i,j))
6214 ! ELSE
6215 ! GRID1(I,J) = spval
6216 ! END IF
6217  ENDDO
6218  ENDDO
6219  if(grib=='grib2') then
6220  cfld=cfld+1
6221  fld_info(cfld)%ifld=iavblfld(iget(397))
6222 !$omp parallel do private(i,j,ii,jj)
6223  do j=1,jend-jsta+1
6224  jj = jsta+j-1
6225  do i=1,iend-ista+1
6226  ii = ista+i-1
6227  datapd(i,j,cfld) = grid1(ii,jj)
6228  enddo
6229  enddo
6230  endif
6231  ENDIF
6232  END IF
6233  IF(iget(396)>0)THEN
6234 !$omp parallel do private(i,j)
6235  DO j=jsta,jend
6236  DO i=ista,iend
6237  grid1(i,j) = suntime(i,j)
6238  ENDDO
6239  ENDDO
6240  id(1:25) = 0
6241  itsrfc = nint(tsrfc)
6242  IF(itsrfc /= 0) then
6243  ifincr = mod(ifhr,itsrfc)
6244  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
6245  ELSE
6246  ifincr = 0
6247  endif
6248  id(19) = ifhr
6249  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
6250  id(20) = 3
6251  IF (ifincr==0) THEN
6252  id(18) = ifhr-itsrfc
6253  ELSE
6254  id(18) = ifhr-ifincr
6255  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
6256  ENDIF
6257  IF (id(18)<0) id(18) = 0
6258  if(grib=='grib2') then
6259  cfld=cfld+1
6260  fld_info(cfld)%ifld=iavblfld(iget(396))
6261  if(itsrfc>0) then
6262  fld_info(cfld)%ntrange=1
6263  else
6264  fld_info(cfld)%ntrange=0
6265  endif
6266  fld_info(cfld)%tinvstat=ifhr-id(18)
6267 !$omp parallel do private(i,j,ii,jj)
6268  do j=1,jend-jsta+1
6269  jj = jsta+j-1
6270  do i=1,iend-ista+1
6271  ii = ista+i-1
6272  datapd(i,j,cfld) = grid1(ii,jj)
6273  enddo
6274  enddo
6275  endif
6276  ENDIF
6277 
6278  IF(iget(517)>0)THEN
6279 !$omp parallel do private(i,j)
6280  DO j=jsta,jend
6281  DO i=ista,iend
6282  grid1(i,j) = avgpotevp(i,j)
6283  ENDDO
6284  ENDDO
6285  id(1:25) = 0
6286  itsrfc = nint(tsrfc)
6287  IF(itsrfc /= 0) then
6288  ifincr = mod(ifhr,itsrfc)
6289  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itsrfc*60)
6290  ELSE
6291  ifincr = 0
6292  endif
6293  id(19) = ifhr
6294  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
6295  id(20) = 3
6296  IF (ifincr==0) THEN
6297  id(18) = ifhr-itsrfc
6298  ELSE
6299  id(18) = ifhr-ifincr
6300  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
6301  ENDIF
6302  IF (id(18)<0) id(18) = 0
6303  if(grib=='grib2') then
6304  cfld=cfld+1
6305  fld_info(cfld)%ifld=iavblfld(iget(517))
6306  if(itsrfc>0) then
6307  fld_info(cfld)%ntrange=1
6308  else
6309  fld_info(cfld)%ntrange=0
6310  endif
6311  fld_info(cfld)%tinvstat=ifhr-id(18)
6312 !$omp parallel do private(i,j,ii,jj)
6313  do j=1,jend-jsta+1
6314  jj = jsta+j-1
6315  do i=1,iend-ista+1
6316  ii = ista+i-1
6317  datapd(i,j,cfld) = grid1(ii,jj)
6318  enddo
6319  enddo
6320  endif
6321  ENDIF
6322 
6323 !
6324 !
6325 ! MODEL TOP REQUESTED BY CMAQ
6326  IF (iget(282)>0) THEN
6327 !$omp parallel do private(i,j)
6328  DO j=jsta,jend
6329  DO i=ista,iend
6330  grid1(i,j) = pt
6331  ENDDO
6332  ENDDO
6333  if(grib=='grib2') then
6334  cfld=cfld+1
6335  fld_info(cfld)%ifld=iavblfld(iget(282))
6336  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6337  endif
6338  ENDIF
6339 !
6340 ! PRESSURE THICKNESS REQUESTED BY CMAQ
6341  IF (iget(283)>0) THEN
6342  DO j=jsta,jend
6343  DO i=ista,iend
6344  grid1(i,j)=pdtop
6345  ENDDO
6346  ENDDO
6347  id(1:25) = 0
6348  IF(me == 0)THEN
6349  DO l=1,lm
6350  IF(pmid(1,1,l)>=(pdtop+pt))EXIT
6351  END DO
6352 ! PRINT*,'hybrid boundary ',L
6353  END IF
6354  CALL mpi_bcast(l,1,mpi_integer,0,mpi_comm_comp,irtn)
6355  if(grib=='grib2') then
6356  cfld=cfld+1
6357  fld_info(cfld)%ifld=iavblfld(iget(283))
6358  fld_info(cfld)%lvl1=1
6359  fld_info(cfld)%lvl2=l
6360  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6361  endif
6362  ENDIF
6363 !
6364 ! SIGMA PRESSURE THICKNESS REQUESTED BY CMAQ
6365  IF (iget(273)>0) THEN
6366  DO j=jsta,jend
6367  DO i=ista,iend
6368  grid1(i,j)=pd(i,j)
6369  ENDDO
6370  ENDDO
6371  IF(me == 0)THEN
6372  DO l=1,lm
6373 ! print*,'Debug CMAQ: ',L,PINT(1,1,LM+1),PD(1,1),PINT(1,1,L)
6374  IF((pint(1,1,lm+1)-pd(1,1))<=(pint(1,1,l)+1.00))EXIT
6375  END DO
6376 ! PRINT*,'hybrid boundary ',L
6377  END IF
6378  CALL mpi_bcast(l,1,mpi_integer,0,mpi_comm_comp,irtn)
6379  if(grib=='grib2') then
6380  cfld=cfld+1
6381  fld_info(cfld)%ifld=iavblfld(iget(273))
6382  fld_info(cfld)%lvl1=l
6383  fld_info(cfld)%lvl2=lm+1
6384  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6385  endif
6386  ENDIF
6387 
6388 
6389 ! TIME-AVERAGED EXCHANGE COEFFICIENTS FOR MASS REQUESTED FOR CMAQ
6390  IF (iget(503)>0) THEN
6391  DO j=jsta,jend
6392  DO i=ista,iend
6393  grid1(i,j)=akhsavg(i,j)
6394  ENDDO
6395  ENDDO
6396  id(1:25) = 0
6397  id(02)= 133
6398  id(19) = ifhr
6399  IF (ifhr==0) THEN
6400  id(18) = 0
6401  ELSE
6402  id(18) = ifhr - 1
6403  ENDIF
6404  id(20) = 3
6405  if(grib=='grib2') then
6406  cfld=cfld+1
6407  fld_info(cfld)%ifld=iavblfld(iget(503))
6408  fld_info(cfld)%ntrange=ifhr-id(18)
6409  fld_info(cfld)%tinvstat=1
6410  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6411  endif
6412  ENDIF
6413 
6414 ! TIME-AVERAGED EXCHANGE COEFFICIENTS FOR WIND REQUESTED FOR CMAQ
6415  IF (iget(504)>0) THEN
6416  DO j=jsta,jend
6417  DO i=ista,iend
6418  grid1(i,j)=akmsavg(i,j)
6419  ENDDO
6420  ENDDO
6421  id(1:25) = 0
6422  id(02)= 133
6423  id(19) = ifhr
6424  IF (ifhr==0) THEN
6425  id(18) = 0
6426  ELSE
6427  id(18) = ifhr - 1
6428  ENDIF
6429  id(20) = 3
6430  if(grib=='grib2') then
6431  cfld=cfld+1
6432  fld_info(cfld)%ifld=iavblfld(iget(504))
6433  fld_info(cfld)%ntrange=ifhr-id(18)
6434  fld_info(cfld)%tinvstat=1
6435  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
6436  endif
6437 
6438  ENDIF
6439 
6440  RETURN
6441  END
6442 
6443  subroutine qpf_comp(igetfld,compfile,fcst)
6444 ! Read in QPF threshold for exceedance grid.
6445 ! Calculate exceedance grid.
6446 ! compfile: file name for reference grid.
6447 ! fcst: forecast length in hours.
6448  use ctlblk_mod, only: spval,jsta,jend,im,dtq2,ifhr,ifmin,tprec,grib, &
6449  modelname,jm,cfld,datapd,fld_info,jsta_2l,jend_2u,&
6450  ista,iend,ista_2l,iend_2u
6451  use rqstfld_mod, only: iget, id, lvls, iavblfld
6452  use grib2_module, only: read_grib2_head, read_grib2_sngle
6453  use vrbls2d, only: avgprec, avgprec_cont
6454  implicit none
6455  character(len=256), intent(in) :: compfile
6456  integer, intent(in) :: igetfld,fcst
6457  integer :: trange,invstat
6458  real, dimension(IM,JM) :: outgrid
6459 
6460  real, allocatable, dimension(:,:) :: mscvalue
6461 
6462  integer :: nx, ny, nz, ntot, mscnlon, mscnlat, height
6463  integer :: itprec, ifincr
6464  real :: rlonmin, rlatmax
6465  real*8 rdx, rdy
6466 
6467  logical :: file_exists
6468 
6469  integer :: i, j, k, ii, jj
6470 
6471 ! Read in reference grid.
6472  INQUIRE(file=compfile, exist=file_exists)
6473  if (file_exists) then
6474  call read_grib2_head(compfile,nx,ny,nz,rlonmin,rlatmax,&
6475  rdx,rdy)
6476  mscnlon=nx
6477  mscnlat=ny
6478  if (.not. allocated(mscvalue)) then
6479  allocate(mscvalue(mscnlon,mscnlat))
6480  endif
6481  ntot = nx*ny
6482  call read_grib2_sngle(compfile,ntot,height,mscvalue)
6483  else
6484  write(*,*) 'WARNING: FFG file not available for hour: ', fcst
6485  endif
6486 
6487 ! Set GRIB variables.
6488  id(1:25) = 0
6489  itprec = nint(tprec)
6490  if (itprec /= 0) then
6491  ifincr = mod(ifhr,itprec)
6492  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itprec*60)
6493  else
6494  ifincr = 0
6495  endif
6496  id(18) = 0
6497  id(19) = ifhr
6498  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
6499  id(20) = 4
6500  IF (ifincr==0) THEN
6501  id(18) = ifhr-itprec
6502  ELSE
6503  id(18) = ifhr-ifincr
6504  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
6505  ENDIF
6506 
6507 ! Calculate exceedance grid.
6508  IF(modelname == 'GFS' .OR. modelname == 'FV3R') THEN
6509 ! !$omp parallel do private(i,j)
6510  IF (file_exists) THEN
6511  DO j=jsta,jend
6512  DO i=ista,iend
6513  IF (ifhr .EQ. 0 .OR. fcst .EQ. 0) THEN
6514  outgrid(i,j) = 0.0
6515  ELSE IF (mscvalue(i,j) .LE. 0.0) THEN
6516  outgrid(i,j) = 0.0
6517  ELSE IF (fcst .EQ. 1 .AND. avgprec(i,j)*float(id(19)-id(18))*3600.*1000./dtq2 .GT. mscvalue(i,j)) THEN
6518  outgrid(i,j) = 1.0
6519  ELSE IF (fcst .GT. 1 .AND. avgprec_cont(i,j)*float(ifhr)*3600.*1000./dtq2 .GT. mscvalue(i,j)) THEN
6520  outgrid(i,j) = 1.0
6521  ELSE
6522  outgrid(i,j) = 0.0
6523  ENDIF
6524  ENDDO
6525  ENDDO
6526  ELSE
6527  outgrid = 0.0*avgprec
6528  ENDIF
6529  ENDIF
6530 ! write(*,*) 'FFG MAX, MIN:', &
6531 ! maxval(mscValue),minval(mscValue)
6532  IF (id(18).LT.0) id(18) = 0
6533 
6534 ! Set GRIB2 variables.
6535  IF(fcst .EQ. 1) THEN
6536  IF(itprec>0) THEN
6537  trange = (ifhr-id(18))/itprec
6538  ELSE
6539  trange = 0
6540  ENDIF
6541  invstat = itprec
6542  IF(trange .EQ. 0) THEN
6543  IF (ifhr .EQ. 0) THEN
6544  invstat = 0
6545  ELSE
6546  invstat = 1
6547  ENDIF
6548  trange = 1
6549  ENDIF
6550  ELSE
6551  trange = 1
6552  IF (ifhr .EQ. fcst) THEN
6553  invstat = fcst
6554  ELSE
6555  invstat = 0
6556  ENDIF
6557  ENDIF
6558 
6559  IF(grib=='grib2') then
6560  cfld=cfld+1
6561  fld_info(cfld)%ifld=iavblfld(iget(igetfld))
6562  fld_info(cfld)%ntrange=trange
6563  fld_info(cfld)%tinvstat=invstat
6564 !$omp parallel do private(i,j,ii,jj)
6565  do j=1,jend-jsta+1
6566  jj = jsta+j-1
6567  do i=1,iend-ista+1
6568  ii = ista+i-1
6569  datapd(i,j,cfld) = outgrid(ii,jj)
6570  enddo
6571  enddo
6572  endif
6573 
6574  RETURN
6575 
6576  end subroutine qpf_comp
masks
Definition: MASKS_mod.f:1
physcons_post
Definition: physcons.f:1
params_mod
Definition: params.F:1
grib2_module
Definition: grib2_module.f:1
soil
Definition: SOIL_mod.f:1
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:345
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27