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