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