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