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