UPP (develop)
Loading...
Searching...
No Matches
MDL2AGL.f
Go to the documentation of this file.
1
2!
50 SUBROUTINE mdl2agl
51
52!
53!
54 use vrbls3d, only: zmid, zint, dbz, dbzr, dbzi, dbzc, uh, vh, pmid, t, q, ref_10cm
55 use vrbls2d, only: refd_max, up_heli_max, up_heli_max16, grpl_max, &
56 up_heli_min, up_heli_min16, up_heli_max02, &
57 up_heli_min02, up_heli_max03, up_heli_min03, &
58 rel_vort_max, rel_vort_max01, hail_max2d, hail_maxk1,&
59 hail_maxhailcast,refdm10c_max, rel_vort_maxhy1, &
60 ltg1_max, ltg2_max, ltg3_max, up_heli, up_heli16, &
61 nci_ltg, nca_ltg, nci_wq, nca_wq, nci_refd, nca_refd,&
62 u10, v10, u10h, v10h
63 use masks, only: lmh, lmv
64 use params_mod, only: dbzmin, small, eps, rd
65 use ctlblk_mod, only: spval, lm, modelname, grib, cfld, fld_info, datapd,&
66 ifhr, global, jsta_m, jend_m, mpi_comm_comp, &
67 jsta_2l, jend_2u, im, jm, jsta, jend, imp_physics, &
68 ista, iend, ista_2l, iend_2u, ista_m, iend_m
69 use rqstfld_mod, only: iget, lvls, iavblfld, lvlsxml, id
70 use gridspec_mod, only: gridtype
71!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 implicit none
73 include "mpif.h"
74!
75! INCLUDE MODEL DIMENSIONS. SET/DERIVE OTHER PARAMETERS.
76! GAMMA AND RGAMOG ARE USED IN THE EXTRAPOLATION OF VIRTUAL
77! TEMPERATURES BEYOND THE UPPER OF LOWER LIMITS OF DATA.
78!
79 integer,PARAMETER :: LAGL=2,lagl2=1
80!
81! DECLARE VARIABLES.
82!
83 LOGICAL IOOMG,IOALL
84 REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: grid1
85 REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: UAGL, VAGL, tagl, pagl, qagl
86!
87 INTEGER,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: NL1X
88 integer,dimension(jm) :: IHE, IHW
89 INTEGER LXXX,IERR, maxll, minll
90 INTEGER ISTART,ISTOP,JSTART,JSTOP
91!
92!
93!--- Definition of the following 2D (horizontal) dummy variables
94!
95! C1D - total condensate
96! QW1 - cloud water mixing ratio
97! QI1 - cloud ice mixing ratio
98! QR1 - rain mixing ratio
99! QS1 - snow mixing ratio
100! DBZ1 - radar reflectivity
101! DBZR1 - radar reflectivity from rain
102! DBZI1 - radar reflectivity from ice (snow + graupel + sleet)
103! DBZC1 - radar reflectivity from parameterized convection (bogused)
104!
105! REAL C1D(IM,JM),QW1(IM,JM),QI1(IM,JM),QR1(IM,JM)
106! &, QS1(IM,JM) ,DBZ1(IM,JM)
107 REAL,dimension(ista:iend,jsta:jend) :: DBZ1, DBZR1, DBZI1, DBZC1, dbz1log
108 real,dimension(lagl) :: ZAGL
109 real,dimension(lagl2) :: ZAGL2, ZAGL3
110 real PAGLU,PAGLL,TAGLU,TAGLL,QAGLU,QAGLL, pv, rho
111
112 integer I,J,L,II,JJ,LP,LL,LLMH,ie,iw,jn,js,iget1,iget2,iget3,iget4
113 real UAGLL,UAGLU,VAGLL,VAGLU,FACT,ZDUM
114!
115!
116!******************************************************************************
117!
118! START MDL2P.
119!
120! SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
121!
122!---------------------------------------------------------------
123
124
125 zagl(1) = 4000.
126 zagl(2) = 1000.
127 zagl2(1) = 609.6 ! 2000 ft
128! CRA
129 zagl3(1) = 80.
130! CRA
131
132!
133! *** PART I ***
134!
135! VERTICAL INTERPOLATION OF EVERYTHING ELSE. EXECUTE ONLY
136! IF THERE'S SOMETHING WE WANT.
137!
138 IF (iget(253)>0 .OR. iget(279)>0 .OR. iget(280)>0 .OR. &
139 & iget(281)>0 ) THEN
140!
141!---------------------------------------------------------------------
142!***
143!*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL
144!*** INTERPOLATION ABOVE GROUND NOW.
145!***
146!
147 DO 310 lp=1,lagl
148 iget1 = -1 ; iget2 = -1 ; iget3 = -1 ; iget4 = -1
149 if (iget(253) > 0) iget1 = lvls(lp,iget(253))
150 if (iget(279) > 0) iget2 = lvls(lp,iget(279))
151 if (iget(280) > 0) iget3 = lvls(lp,iget(280))
152 if (iget(281) > 0) iget4 = lvls(lp,iget(281))
153 IF (iget1 > 0 .or. iget2 > 0 .or. iget3 > 0 .or. iget4 > 0) then
154!
155 jj=float(jsta+jend)/2.0
156 ii=float(ista+iend)/3.0
157
158 DO j=jsta,jend
159 DO i=ista,iend
160 dbz1(i,j) = spval
161 dbzr1(i,j) = spval
162 dbzi1(i,j) = spval
163 dbzc1(i,j) = spval
164!
165!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
166!*** THE AGL LEVEL TO WHICH WE ARE INTERPOLATING.
167!
168 llmh = nint(lmh(i,j))
169 nl1x(i,j) = llmh+1
170 DO l=llmh,2,-1
171 zdum = zmid(i,j,l)-zint(i,j,llmh+1)
172 IF(zdum >= zagl(lp)) THEN
173 nl1x(i,j) = l+1
174 exit
175 ENDIF
176 ENDDO
177!
178! IF THE AGL LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
179! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
180! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
181! WILL EXTRAPOLATE TO THAT POINT
182!
183 IF(nl1x(i,j) == (llmh+1) .AND. zagl(lp) > 0.) THEN
184 nl1x(i,j) = lm
185 ENDIF
186!
187! if(NL1X(I,J)==LMP1)print*,'Debug: NL1X=LMP1 AT '
188! 1 ,i,j,lp
189 ENDDO
190 ENDDO
191!
192!mptest IF(NHOLD==0)GO TO 310
193!
194!!$omp parallel do
195!!$omp& private(nn,i,j,ll,fact,qsat,rhl)
196!hc DO 220 NN=1,NHOLD
197!hc I=IHOLD(NN)
198!hc J=JHOLD(NN)
199! DO 220 J=JSTA,JEND
200
201 DO j=jsta,jend
202 DO i=ista,iend
203 ll = nl1x(i,j)
204!---------------------------------------------------------------------
205!*** VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
206!*** HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
207!---------------------------------------------------------------------
208!
209!HC IF(NL1X(I,J)<=LM)THEN
210 llmh = nint(lmh(i,j))
211 IF(nl1x(i,j)<=llmh)THEN
212 IF(zmid(i,j,ll)<spval.and.zmid(i,j,ll-1)<spval)THEN
213!
214!---------------------------------------------------------------------
215! INTERPOLATE LINEARLY IN LOG(P)
216!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
217!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
218!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
219!---------------------------------------------------------------------
220!
221
222! FACT=(ALSL(LP)-ALOG(PMID(I,J,LL)))/
223! & (ALOG(PMID(I,J,LL))-ALOG(PMID(I,J,LL-1)))
224 zdum=zagl(lp)+zint(i,j,nint(lmh(i,j))+1)
225 fact=(zdum-zmid(i,j,ll))/(zmid(i,j,ll)-zmid(i,j,ll-1))
226!
227! KRF: Use arw/nmm output if thompson
228 if (imp_physics==8) then
229 dbz1(i,j)=ref_10cm(i,j,ll)+(ref_10cm(i,j,ll)-ref_10cm(i,j,ll-1))*fact
230 else
231 dbz1(i,j)=dbz(i,j,ll)+(dbz(i,j,ll)-dbz(i,j,ll-1))*fact
232 end if
233 ! DBZ1(I,J) = DBZ(I,J,LL) + (DBZ(I,J,LL)-DBZ(I,J,LL-1))*FACT
234 dbzr1(i,j) = dbzr(i,j,ll) + (dbzr(i,j,ll)-dbzr(i,j,ll-1))*fact
235 dbzi1(i,j) = dbzi(i,j,ll) + (dbzi(i,j,ll)-dbzi(i,j,ll-1))*fact
236 dbzc1(i,j) = dbzc(i,j,ll) + (dbzc(i,j,ll)-dbzc(i,j,ll-1))*fact
237 if(modelname=='RAPR') then
238 if(dbz1(i,j)>0.) then
239 dbz1log(i,j)= 10.*log10(dbz1(i,j))
240 else
241 dbz1log(i,j)= -100.
242 endif
243 endif
244! IF(I==ii.and.j==jj)print*,'Debug AGL RADAR REF',
245! & i,j,ll,zagl(lp),ZINT(I,J,NINT(LMH(I,J))+1)
246! & ,ZMID(I,J,LL-1),ZMID(I,J,LL)
247! & ,DBZ(I,J,LL-1),DBZ(I,J,LL),DBZ1(I,J)
248! & ,DBZR(I,J,LL-1),DBZR(I,J,LL),DBZR1(I,J)
249! & ,DBZI(I,J,LL-1),DBZI(I,J,LL),DBZI1(I,J)
250! & ,DBZC(I,J,LL-1),DBZC(I,J,LL),DBZC1(I,J)
251 if(modelname=='RAPR') then
252 dbz1log(i,j)=max(dbz1log(i,j),dbzmin)
253 else
254 dbz1(i,j)=max(dbz1(i,j),dbzmin)
255 endif
256 dbzr1(i,j) = max(dbzr1(i,j),dbzmin)
257 dbzi1(i,j) = max(dbzi1(i,j),dbzmin)
258 dbzc1(i,j) = max(dbzc1(i,j),dbzmin)
259 ENDIF !end ZMID(I,J,LL)<spval
260!
261! FOR UNDERGROUND AGL LEVELS, ASSUME TEMPERATURE TO CHANGE
262! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
263! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
264! GOUND
265 ELSE
266 dbz1log(i,j) = dbzmin
267 dbzr1(i,j) = dbzmin
268 dbzi1(i,j) = dbzmin
269 dbzc1(i,j) = dbzmin
270 END IF
271 enddo
272 enddo
273!
274!
275!---------------------------------------------------------------------
276! *** PART II ***
277!---------------------------------------------------------------------
278!
279! OUTPUT SELECTED FIELDS.
280!
281!---------------------------------------------------------------------
282!
283!
284!--- Radar Reflectivity
285 IF((iget(253)>0) )THEN
286 if(modelname=='RAPR') then
287 DO j=jsta,jend
288 DO i=ista,iend
289 grid1(i,j)=dbz1log(i,j)
290 ENDDO
291 ENDDO
292 else
293 DO j=jsta,jend
294 DO i=ista,iend
295 grid1(i,j)=dbz1(i,j)
296 ENDDO
297 ENDDO
298 endif
299 if(grib=='grib2') then
300 cfld=cfld+1
301 fld_info(cfld)%ifld=iavblfld(iget(253))
302 fld_info(cfld)%lvl=lvlsxml(lp,iget(253))
303 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
304 endif
305 END IF
306!--- Radar reflectivity from rain
307 IF((iget(279)>0) )THEN
308 DO j=jsta,jend
309 DO i=ista,iend
310 grid1(i,j)=dbzr1(i,j)
311 ENDDO
312 ENDDO
313 if(grib=='grib2') then
314 cfld=cfld+1
315 fld_info(cfld)%ifld=iavblfld(iget(279))
316 fld_info(cfld)%lvl=lvlsxml(lp,iget(279))
317 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
318 endif
319 END IF
320!--- Radar reflectivity from all ice habits (snow + graupel + sleet, etc.)
321 IF((iget(280)>0) )THEN
322 DO j=jsta,jend
323 DO i=ista,iend
324 grid1(i,j)=dbzi1(i,j)
325 ENDDO
326 ENDDO
327 if(grib=='grib2') then
328 cfld=cfld+1
329 fld_info(cfld)%ifld=iavblfld(iget(280))
330 fld_info(cfld)%lvl=lvlsxml(lp,iget(280))
331 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
332 endif
333 END IF
334!--- Radar reflectivity from parameterized convection
335 IF((iget(281)>0) )THEN
336 DO j=jsta,jend
337 DO i=ista,iend
338 grid1(i,j)=dbzc1(i,j)
339 ENDDO
340 ENDDO
341 if(grib=='grib2') then
342 cfld=cfld+1
343 fld_info(cfld)%ifld=iavblfld(iget(281))
344 fld_info(cfld)%lvl=lvlsxml(lp,iget(281))
345 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
346 endif
347 END IF
348!
349 ENDIF ! FOR LEVEL
350!
351!*** END OF MAIN VERTICAL LOOP
352!
353 310 CONTINUE
354!*** ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
355!
356 ENDIF
357! SRD
358 lp=1
359!--- Max Derived Radar Reflectivity
360 IF((iget(421)>0) )THEN
361 DO j=jsta,jend
362 DO i=ista,iend
363 grid1(i,j)=refd_max(i,j)
364 ENDDO
365 ENDDO
366 if(grib=='grib2') then
367 cfld=cfld+1
368 fld_info(cfld)%ifld=iavblfld(iget(421))
369 fld_info(cfld)%lvl=lvlsxml(lp,iget(421))
370 fld_info(cfld)%tinvstat=1
371 if (ifhr > 0) then
372 fld_info(cfld)%tinvstat=1
373 else
374 fld_info(cfld)%tinvstat=0
375 endif
376 fld_info(cfld)%ntrange=1
377 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
378 endif
379 END IF
380
381!--- Max Derived Radar Reflectivity at -10C
382 IF((iget(785)>0) )THEN
383 DO j=jsta,jend
384 DO i=ista,iend
385 grid1(i,j)=refdm10c_max(i,j)
386 ENDDO
387 ENDDO
388 if(grib=='grib2') then
389 cfld=cfld+1
390 fld_info(cfld)%ifld=iavblfld(iget(785))
391 fld_info(cfld)%lvl=lvlsxml(lp,iget(785))
392 if (ifhr > 0) then
393 fld_info(cfld)%tinvstat=1
394 else
395 fld_info(cfld)%tinvstat=0
396 endif
397 fld_info(cfld)%ntrange=1
398 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
399 endif
400 END IF
401
402!--- Max Updraft Helicity
403 IF((iget(420)>0) )THEN
404 DO j=jsta,jend
405 DO i=ista,iend
406 grid1(i,j)=up_heli_max(i,j)
407 ENDDO
408 ENDDO
409 if(grib=='grib2') then
410 cfld=cfld+1
411 fld_info(cfld)%ifld=iavblfld(iget(420))
412 fld_info(cfld)%lvl=lvlsxml(lp,iget(420))
413 if (ifhr > 0) then
414 fld_info(cfld)%tinvstat = 1
415 else
416 fld_info(cfld)%tinvstat = 0
417 endif
418 fld_info(cfld)%ntrange = 1
419 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
420 endif
421 END IF
422
423!--- Max Updraft Helicity 1-6 km
424 IF((iget(700)>0) )THEN
425 DO j=jsta,jend
426 DO i=ista,iend
427 grid1(i,j)=up_heli_max16(i,j)
428 ENDDO
429 ENDDO
430 if(grib=='grib2') then
431 cfld=cfld+1
432 fld_info(cfld)%ifld=iavblfld(iget(700))
433 fld_info(cfld)%lvl=lvlsxml(lp,iget(700))
434 if (ifhr == 0) then
435 fld_info(cfld)%tinvstat = 0
436 else
437 fld_info(cfld)%tinvstat = 1
438 endif
439 fld_info(cfld)%ntrange = 1
440 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
441 endif
442 END IF
443
444!--- Min Updraft Helicity
445 IF((iget(786)>0) )THEN
446 DO j=jsta,jend
447 DO i=ista,iend
448 grid1(i,j)=up_heli_min(i,j)
449 ENDDO
450 ENDDO
451 if(grib=='grib2') then
452 cfld=cfld+1
453 fld_info(cfld)%ifld=iavblfld(iget(786))
454 fld_info(cfld)%lvl=lvlsxml(lp,iget(786))
455 if (ifhr > 0) then
456 fld_info(cfld)%tinvstat = 1
457 else
458 fld_info(cfld)%tinvstat = 0
459 endif
460 fld_info(cfld)%ntrange = 1
461 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
462 endif
463 END IF
464
465!--- Min Updraft Helicity 1-6 km
466 IF((iget(787)>0) )THEN
467 DO j=jsta,jend
468 DO i=ista,iend
469 grid1(i,j)=up_heli_min16(i,j)
470 ENDDO
471 ENDDO
472 if(grib=='grib2') then
473 cfld=cfld+1
474 fld_info(cfld)%ifld=iavblfld(iget(787))
475 fld_info(cfld)%lvl=lvlsxml(lp,iget(787))
476 if (ifhr == 0) then
477 fld_info(cfld)%tinvstat = 0
478 else
479 fld_info(cfld)%tinvstat = 1
480 endif
481 fld_info(cfld)%ntrange = 1
482 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
483 endif
484 END IF
485
486!--- Max Updraft Helicity 0-2 km
487 IF((iget(788)>0) )THEN
488 DO j=jsta,jend
489 DO i=ista,iend
490 grid1(i,j)=up_heli_max02(i,j)
491 ENDDO
492 ENDDO
493 if(grib=='grib2') then
494 cfld=cfld+1
495 fld_info(cfld)%ifld=iavblfld(iget(788))
496 fld_info(cfld)%lvl=lvlsxml(lp,iget(788))
497 if (ifhr > 0) then
498 fld_info(cfld)%tinvstat = 1
499 else
500 fld_info(cfld)%tinvstat = 0
501 endif
502 fld_info(cfld)%ntrange = 1
503 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
504 endif
505 END IF
506!--- Min Updraft Helicity 0-2 km
507 IF((iget(789)>0) )THEN
508 DO j=jsta,jend
509 DO i=ista,iend
510 grid1(i,j)=up_heli_min02(i,j)
511 ENDDO
512 ENDDO
513 if(grib=='grib2') then
514 cfld=cfld+1
515 fld_info(cfld)%ifld=iavblfld(iget(789))
516 fld_info(cfld)%lvl=lvlsxml(lp,iget(789))
517 if (ifhr == 0) then
518 fld_info(cfld)%tinvstat = 0
519 else
520 fld_info(cfld)%tinvstat = 1
521 endif
522 fld_info(cfld)%ntrange = 1
523 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
524 endif
525 END IF
526
527!--- Max Updraft Helicity 0-3 km
528 IF((iget(790)>0) )THEN
529 DO j=jsta,jend
530 DO i=ista,iend
531 grid1(i,j)=up_heli_max03(i,j)
532 ENDDO
533 ENDDO
534 if(grib=='grib2') then
535 cfld=cfld+1
536 fld_info(cfld)%ifld=iavblfld(iget(790))
537 fld_info(cfld)%lvl=lvlsxml(lp,iget(790))
538 if (ifhr > 0) then
539 fld_info(cfld)%tinvstat = 1
540 else
541 fld_info(cfld)%tinvstat = 0
542 endif
543 fld_info(cfld)%ntrange = 1
544 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
545 endif
546 END IF
547
548!--- Min Updraft Helicity 0-3 km
549 IF((iget(791)>0) )THEN
550 DO j=jsta,jend
551 DO i=ista,iend
552 grid1(i,j)=up_heli_min03(i,j)
553 ENDDO
554 ENDDO
555 if(grib=='grib2') then
556 cfld=cfld+1
557 fld_info(cfld)%ifld=iavblfld(iget(791))
558 fld_info(cfld)%lvl=lvlsxml(lp,iget(791))
559 if (ifhr == 0) then
560 fld_info(cfld)%tinvstat = 0
561 else
562 fld_info(cfld)%tinvstat = 1
563 endif
564 fld_info(cfld)%ntrange = 1
565 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
566 endif
567 END IF
568
569!--- Max Relative Vertical Vorticity 0-2 km
570 IF((iget(792)>0) )THEN
571 DO j=jsta,jend
572 DO i=ista,iend
573 grid1(i,j)=rel_vort_max(i,j)
574 ENDDO
575 ENDDO
576 if(grib=='grib2') then
577 cfld=cfld+1
578 fld_info(cfld)%ifld=iavblfld(iget(792))
579 fld_info(cfld)%lvl=lvlsxml(lp,iget(792))
580 if (ifhr > 0) then
581 fld_info(cfld)%tinvstat = 1
582 else
583 fld_info(cfld)%tinvstat = 0
584 endif
585 fld_info(cfld)%ntrange = 1
586 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
587 endif
588 END IF
589
590!--- Max Relative Vertical Vorticity 0-1 km
591 IF((iget(793)>0) )THEN
592 DO j=jsta,jend
593 DO i=ista,iend
594 grid1(i,j)=rel_vort_max01(i,j)
595 ENDDO
596 ENDDO
597 if(grib=='grib2') then
598 cfld=cfld+1
599 fld_info(cfld)%ifld=iavblfld(iget(793))
600 fld_info(cfld)%lvl=lvlsxml(lp,iget(793))
601 if (ifhr > 0) then
602 fld_info(cfld)%tinvstat = 1
603 else
604 fld_info(cfld)%tinvstat = 0
605 endif
606 fld_info(cfld)%ntrange = 1
607 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
608 endif
609 END IF
610!--- Max Relative Vertical Vorticity @ hybrid level 1
611 IF((iget(890)>0) )THEN
612 DO j=jsta,jend
613 DO i=ista,iend
614 grid1(i,j)=rel_vort_maxhy1(i,j)
615 ENDDO
616 ENDDO
617 if(grib=='grib2') then
618 cfld=cfld+1
619 fld_info(cfld)%ifld=iavblfld(iget(890))
620 fld_info(cfld)%lvl=lvlsxml(lp,iget(890))
621 if (ifhr > 0) then
622 fld_info(cfld)%tinvstat = 1
623 else
624 fld_info(cfld)%tinvstat = 0
625 endif
626 fld_info(cfld)%ntrange = 1
627 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
628 endif
629 END IF
630
631!--- Max Hail Diameter in Column
632 IF((iget(794)>0) )THEN
633 DO j=jsta,jend
634 DO i=ista,iend
635 grid1(i,j)=hail_max2d(i,j)
636 ENDDO
637 ENDDO
638 if(grib=='grib2') then
639 cfld=cfld+1
640 fld_info(cfld)%ifld=iavblfld(iget(794))
641 fld_info(cfld)%lvl=lvlsxml(lp,iget(794))
642 if (ifhr == 0) then
643 fld_info(cfld)%tinvstat = 0
644 else
645 fld_info(cfld)%tinvstat = 1
646 endif
647 fld_info(cfld)%ntrange = 1
648 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
649 endif
650 END IF
651
652!--- Max Hail Diameter at k=1
653 IF((iget(795)>0) )THEN
654 DO j=jsta,jend
655 DO i=ista,iend
656 grid1(i,j)=hail_maxk1(i,j)
657 ENDDO
658 ENDDO
659 if(grib=='grib2') then
660 cfld=cfld+1
661 fld_info(cfld)%ifld=iavblfld(iget(795))
662 fld_info(cfld)%lvl=lvlsxml(lp,iget(795))
663 if (ifhr == 0) then
664 fld_info(cfld)%tinvstat = 0
665 else
666 fld_info(cfld)%tinvstat = 1
667 endif
668 fld_info(cfld)%ntrange = 1
669 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
670 endif
671 END IF
672
673!--- Max hail diameter at surface from HAILCAST algorithm (HRRR and RRFS applications)
674 IF((iget(728)>0) )THEN
675 DO j=jsta,jend
676 DO i=ista,iend
677 grid1(i,j)=hail_maxhailcast(i,j)/1000.0 ! convert mm to m
678 ENDDO
679 ENDDO
680 if(grib=='grib2') then
681 cfld=cfld+1
682 fld_info(cfld)%ifld=iavblfld(iget(728))
683 fld_info(cfld)%lvl=lvlsxml(lp,iget(728))
684 if (ifhr == 0) then
685 fld_info(cfld)%tinvstat = 0
686 else
687 fld_info(cfld)%tinvstat = 1
688 endif
689 fld_info(cfld)%ntrange = 1
690 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
691 endif
692 END IF
693
694!--- Max Column Integrated Graupel
695 IF((iget(429)>0) )THEN
696 DO j=jsta,jend
697 DO i=ista,iend
698 grid1(i,j)=grpl_max(i,j)
699 ENDDO
700 ENDDO
701 if(grib=='grib2') then
702 cfld=cfld+1
703 fld_info(cfld)%ifld=iavblfld(iget(429))
704 fld_info(cfld)%lvl=lvlsxml(lp,iget(429))
705 if (ifhr == 0) then
706 fld_info(cfld)%tinvstat = 0
707 else
708 fld_info(cfld)%tinvstat = 1
709 endif
710 fld_info(cfld)%ntrange = 1
711 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
712 endif
713 END IF
714
715!--- Max Lightning Threat 1
716 IF((iget(702)>0) )THEN
717 DO j=jsta,jend
718 DO i=ista,iend
719 grid1(i,j)=ltg1_max(i,j)
720 ENDDO
721 ENDDO
722 if(grib=='grib2') then
723 cfld=cfld+1
724 fld_info(cfld)%ifld=iavblfld(iget(702))
725 fld_info(cfld)%lvl=lvlsxml(lp,iget(702))
726 if (ifhr == 0) then
727 fld_info(cfld)%tinvstat = 0
728 else
729 fld_info(cfld)%tinvstat = 1
730 endif
731 fld_info(cfld)%ntrange = 1
732 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
733 endif
734 END IF
735
736!--- Max Lightning Threat 2
737 IF((iget(703)>0) )THEN
738 DO j=jsta,jend
739 DO i=ista,iend
740 grid1(i,j)=ltg2_max(i,j)
741 ENDDO
742 ENDDO
743 if(grib=='grib2') then
744 cfld=cfld+1
745 fld_info(cfld)%ifld=iavblfld(iget(703))
746 fld_info(cfld)%lvl=lvlsxml(lp,iget(703))
747 if (ifhr == 0) then
748 fld_info(cfld)%tinvstat = 0
749 else
750 fld_info(cfld)%tinvstat = 1
751 endif
752 fld_info(cfld)%ntrange = 1
753 datapd(1:iend-ista+1,1:jend-jsta+1,cfld) = grid1(ista:iend,jsta:jend)
754 endif
755 END IF
756
757!--- Max Lightning Threat 3
758 IF((iget(704)>0) )THEN
759 DO j=jsta,jend
760 DO i=ista,iend
761 grid1(i,j)=ltg3_max(i,j)
762 ENDDO
763 ENDDO
764 if(grib=='grib2') then
765 cfld=cfld+1
766 fld_info(cfld)%ifld=iavblfld(iget(704))
767 fld_info(cfld)%lvl=lvlsxml(lp,iget(704))
768 if (ifhr == 0) then
769 fld_info(cfld)%tinvstat = 0
770 else
771 fld_info(cfld)%tinvstat = 1
772 endif
773 fld_info(cfld)%ntrange = 1
774 do j=1,jend-jsta+1
775 jj = jsta+j-1
776 do i=1,iend-ista+1
777 ii = ista+i-1
778 datapd(i,j,cfld) = grid1(ii,jj)
779 enddo
780 enddo
781 endif
782 !datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
783 END IF
784
785!--- GSD Updraft Helicity
786 IF((iget(727)>0) )THEN
787 DO j=jsta,jend
788 DO i=ista,iend
789 grid1(i,j)=up_heli(i,j)
790 ENDDO
791 ENDDO
792 if(grib=='grib2') then
793 cfld=cfld+1
794 fld_info(cfld)%ifld=iavblfld(iget(727))
795 fld_info(cfld)%lvl=lvlsxml(lp,iget(727))
796 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
797 endif
798 END IF
799
800!--- Updraft Helicity 1-6 km layer
801 IF((iget(701)>0) )THEN
802 DO j=jsta,jend
803 DO i=ista,iend
804 grid1(i,j)=up_heli16(i,j)
805 ENDDO
806 ENDDO
807 if(grib=='grib2') then
808 cfld=cfld+1
809 fld_info(cfld)%ifld=iavblfld(iget(701))
810 fld_info(cfld)%lvl=lvlsxml(lp,iget(701))
811 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
812 endif
813 END IF
814
815!--- Convective Initiation Lightning
816 IF((iget(705)>0) )THEN
817 DO j=jsta,jend
818 DO i=ista,iend
819 grid1(i,j)=nci_ltg(i,j)/60.0
820 ENDDO
821 ENDDO
822 if(grib=='grib2') then
823 cfld=cfld+1
824 fld_info(cfld)%ifld=iavblfld(iget(705))
825 fld_info(cfld)%lvl=lvlsxml(lp,iget(705))
826 if (ifhr == 0) then
827 fld_info(cfld)%tinvstat = 0
828 else
829 fld_info(cfld)%tinvstat = 1
830 endif
831 fld_info(cfld)%ntrange = 1
832 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
833 endif
834 END IF
835
836!--- Convective Activity Lightning
837 IF((iget(706)>0) )THEN
838 DO j=jsta,jend
839 DO i=ista,iend
840 grid1(i,j)=nca_ltg(i,j)/60.0
841 ENDDO
842 ENDDO
843 if(grib=='grib2') then
844 cfld=cfld+1
845 fld_info(cfld)%ifld=iavblfld(iget(706))
846 fld_info(cfld)%lvl=lvlsxml(lp,iget(706))
847 if (ifhr == 0) then
848 fld_info(cfld)%tinvstat = 0
849 else
850 fld_info(cfld)%tinvstat = 1
851 endif
852 fld_info(cfld)%ntrange = 1
853 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
854 endif
855 END IF
856
857!--- Convective Initiation Vertical Hydrometeor Flux
858 IF((iget(707)>0) )THEN
859 DO j=jsta,jend
860 DO i=ista,iend
861 grid1(i,j)=nci_wq(i,j)/60.0
862 ENDDO
863 ENDDO
864 if(grib=='grib2') then
865 cfld=cfld+1
866 fld_info(cfld)%ifld=iavblfld(iget(707))
867 fld_info(cfld)%lvl=lvlsxml(lp,iget(707))
868 if (ifhr == 0) then
869 fld_info(cfld)%tinvstat = 0
870 else
871 fld_info(cfld)%tinvstat = 1
872 endif
873 fld_info(cfld)%ntrange = 1
874 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
875 endif
876 END IF
877
878!--- Convective Activity Vertical Hydrometeor Flux
879 IF((iget(708)>0) )THEN
880 DO j=jsta,jend
881 DO i=ista,iend
882 grid1(i,j)=nca_wq(i,j)/60.0
883 ENDDO
884 ENDDO
885 if(grib=='grib2') then
886 cfld=cfld+1
887 fld_info(cfld)%ifld=iavblfld(iget(708))
888 fld_info(cfld)%lvl=lvlsxml(lp,iget(708))
889 if (ifhr == 0) then
890 fld_info(cfld)%tinvstat = 0
891 else
892 fld_info(cfld)%tinvstat = 1
893 endif
894 fld_info(cfld)%ntrange = 1
895 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
896 endif
897 END IF
898
899!--- Convective Initiation Reflectivity
900 IF((iget(709)>0) )THEN
901 DO j=jsta,jend
902 DO i=ista,iend
903 grid1(i,j)=nci_refd(i,j)/60.0
904 ENDDO
905 ENDDO
906 if(grib=='grib2') then
907 cfld=cfld+1
908 fld_info(cfld)%ifld=iavblfld(iget(709))
909 fld_info(cfld)%lvl=lvlsxml(lp,iget(709))
910 if (ifhr == 0) then
911 fld_info(cfld)%tinvstat = 0
912 else
913 fld_info(cfld)%tinvstat = 1
914 endif
915 fld_info(cfld)%ntrange = 1
916 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
917 endif
918 END IF
919
920!--- Convective Activity Reflectivity
921 IF((iget(710)>0) )THEN
922 DO j=jsta,jend
923 DO i=ista,iend
924 grid1(i,j)=nca_refd(i,j)/60.0
925 ENDDO
926 ENDDO
927 if(grib=='grib2') then
928 cfld=cfld+1
929 fld_info(cfld)%ifld=iavblfld(iget(710))
930 fld_info(cfld)%lvl=lvlsxml(lp,iget(710))
931 if (ifhr == 0) then
932 fld_info(cfld)%tinvstat = 0
933 else
934 fld_info(cfld)%tinvstat = 1
935 endif
936 fld_info(cfld)%ntrange = 1
937 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
938 endif
939 END IF
940!
941! SRD
942
943!
944 IF((iget(259)>0) )THEN
945!
946!---------------------------------------------------------------------
947!***
948!*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL
949!*** INTERPOLATION ABOVE GROUND NOW.
950!***
951!
952 iget2 = -1
953 if (iget(253) > 0 ) iget2 = iavblfld(iget(253))
954 iget2 = iget(253)
955 DO 320 lp=1,lagl2
956 iget1 = -1
957 if (iget(259) > 0 ) iget1 = lvls(lp,iget(259))
958 IF(iget1 > 0 .or. iget2 > 0) THEN
959!
960 jj=(jsta+jend)/2
961 ii=(ista+iend)/2
962 DO j=jsta,jend
963 DO i=ista,iend
964 uagl(i,j) = spval
965 vagl(i,j) = spval
966!
967!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
968!*** THE AGL LEVEL TO WHICH WE ARE INTERPOLATING.
969!
970 llmh=nint(lmh(i,j))
971 nl1x(i,j) = llmh+1
972 DO l=llmh,2,-1
973 zdum=zmid(i,j,l)-zint(i,j,llmh+1)
974 IF(zdum >= zagl2(lp))THEN
975 nl1x(i,j)=l+1
976 exit
977 ENDIF
978 ENDDO
979!
980! IF THE AGL LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
981! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
982! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
983! WILL EXTRAPOLATE TO THAT POINT
984!
985 IF(nl1x(i,j) == (llmh+1) .AND. zagl2(lp) > 0.) THEN
986 nl1x(i,j)=lm
987 ENDIF
988!
989! if(NL1X(I,J)==LMP1)print*,'Debug: NL1X=LMP1 AT '
990! 1 ,i,j,lp
991 ENDDO
992 ENDDO
993!
994!mptest IF(NHOLD==0)GO TO 310
995!
996!!$omp parallel do
997!!$omp& private(nn,i,j,ll,fact,qsat,rhl)
998!hc DO 220 NN=1,NHOLD
999!hc I=IHOLD(NN)
1000!hc J=JHOLD(NN)
1001! DO 220 J=JSTA,JEND
1002 DO j=jsta,jend
1003 IF(gridtype=='A')THEN
1004 ihw(j)=-1
1005 ihe(j)=1
1006 ELSE IF(gridtype=='E')THEN
1007 ihw(j)=-mod(j,2)
1008 ihe(j)=ihw(j)+1
1009 END IF
1010 ENDDO
1011 IF(global)then
1012 istart=ista
1013 istop=iend
1014 jstart=jsta
1015 jstop=jend
1016 ELSE
1017 istart=ista_m
1018 istop=iend_m
1019 jstart=jsta_m
1020 jstop=jend_m
1021 END IF
1022 IF(gridtype/='A')THEN
1023! MAXLL=maxval(NL1X)
1024 minll=minval(nl1x)
1025! print*,'MINLL before all reduce= ',MINLL
1026 CALL mpi_allreduce(minll,lxxx,1,mpi_integer,mpi_min,mpi_comm_comp,ierr)
1027 minll=lxxx
1028! print*,'exchange wind in MDL2AGL from ',MINLL
1029 DO ll=minll,lm
1030 call exch(uh(ista_2l:iend_2u,jsta_2l:jend_2u,ll))
1031 call exch(vh(ista_2l:iend_2u,jsta_2l:jend_2u,ll))
1032 END DO
1033 END IF
1034 DO 230 j=jstart,jstop
1035 DO 230 i=istart,istop
1036 ll=nl1x(i,j)
1037!---------------------------------------------------------------------
1038!*** VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
1039!*** HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
1040!---------------------------------------------------------------------
1041!
1042!HC IF(NL1X(I,J)<=LM)THEN
1043 llmh = nint(lmh(i,j))
1044 IF(nl1x(i,j)<=llmh)THEN
1045!
1046!---------------------------------------------------------------------
1047! INTERPOLATE LINEARLY IN LOG(P)
1048!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
1049!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
1050!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
1051!---------------------------------------------------------------------
1052!
1053! FACT=(ALSL(LP)-ALOG(PMID(I,J,LL)))/
1054! & (ALOG(PMID(I,J,LL))-ALOG(PMID(I,J,LL-1)))
1055 zdum=zagl2(lp)+zint(i,j,nint(lmh(i,j))+1)
1056 fact=(zdum-zmid(i,j,ll))/(zmid(i,j,ll)-zmid(i,j,ll-1))
1057!
1058 IF(gridtype=='A')THEN
1059 uaglu=uh(i,j,ll-1)
1060 uagll=uh(i,j,ll)
1061
1062 vaglu=vh(i,j,ll-1)
1063 vagll=vh(i,j,ll)
1064 ELSE IF(gridtype=='E')THEN
1065 uaglu=(uh(i+ihe(j),j,ll-1)+uh(i+ihw(j),j,ll-1)+ &
1066 & uh(i,j-1,ll-1)+uh(i,j+1,ll-1))/4.0
1067 uagll=(uh(i+ihe(j),j,ll)+uh(i+ihw(j),j,ll)+ &
1068 & uh(i,j-1,ll)+uh(i,j+1,ll))/4.0
1069
1070 vaglu=(vh(i+ihe(j),j,ll-1)+vh(i+ihw(j),j,ll-1)+ &
1071 & vh(i,j-1,ll-1)+vh(i,j+1,ll-1))/4.0
1072 vagll=(vh(i+ihe(j),j,ll)+vh(i+ihw(j),j,ll)+ &
1073 & vh(i,j-1,ll)+vh(i,j+1,ll))/4.0
1074 ELSE IF(gridtype=='B')THEN
1075 ie=i
1076 iw=i-1
1077 jn=j
1078 js=j-1
1079 uaglu=(uh(ie,j,ll-1)+uh(iw,j,ll-1)+ &
1080 & uh(ie,js,ll-1)+uh(iw,js,ll-1))/4.0
1081 uagll=(uh(ie,j,ll)+uh(iw,j,ll)+ &
1082 & uh(ie,js,ll)+uh(iw,js,ll))/4.0
1083
1084 vaglu=(vh(ie,j,ll-1)+vh(iw,j,ll-1)+ &
1085 & vh(ie,js,ll-1)+vh(iw,js,ll-1))/4.0
1086 vagll=(vh(ie,j,ll)+vh(iw,j,ll)+ &
1087 & vh(ie,js,ll)+vh(iw,js,ll))/4.0
1088 END IF
1089 uagl(i,j)=uagll+(uagll-uaglu)*fact
1090 vagl(i,j)=vagll+(vagll-vaglu)*fact
1091 IF(i==ii.AND.j==jj)print*, &
1092 & 'DEBUG LLWS: I,J,NL1X,UU,UL,VU,VL,ZSFC,ZMIDU,ZMIDL,U,V= ' &
1093 &, i,j,ll,uaglu,uagll,vaglu,vagll,zint(i,j,nint(lmh(i,j))+1)&
1094 &, zmid(i,j,ll-1),zmid(i,j,ll),uagl(i,j),vagl(i,j) &
1095 &, u10(i,j),v10(i,j)
1096!
1097! FOR UNDERGROUND AGL LEVELS, ASSUME TEMPERATURE TO CHANGE
1098! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
1099! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
1100! GOUND
1101 ELSE
1102 IF(gridtype=='A')THEN
1103 uagl(i,j)=uh(i,j,nint(lmv(i,j)))
1104 vagl(i,j)=vh(i,j,nint(lmv(i,j)))
1105 ELSE IF(gridtype=='E')THEN
1106 uagl(i,j)=(uh(i+ihe(j),j,nint(lmv(i+ihe(j),j))) &
1107 & +uh(i+ihw(j),j,nint(lmv(i+ihw(j),j)))+ &
1108 & uh(i,j-1,nint(lmv(i,j-1)))+uh(i,j+1,nint(lmv(i,j+1))))/4.0
1109 vagl(i,j)=(vh(i+ihe(j),j,nint(lmv(i+ihe(j),j))) &
1110 & +vh(i+ihw(j),j,nint(lmv(i+ihw(j),j)))+ &
1111 & vh(i,j-1,nint(lmv(i,j-1)))+vh(i,j+1,nint(lmv(i,j+1))))/4.0
1112 ELSE IF(gridtype=='B')THEN
1113 ie=i
1114 iw=i-1
1115 jn=j
1116 js=j-1
1117 uagl(i,j)=(uh(ie,j,nint(lmv(ie,j))) &
1118 & +uh(iw,j,nint(lmv(iw,j)))+ &
1119 & uh(ie,js,nint(lmv(ie,js)))+uh(iw,js,nint(lmv(iw,js))))/4.0
1120 vagl(i,j)=(vh(ie,j,nint(lmv(ie,j))) &
1121 & +vh(iw,j,nint(lmv(iw,j)))+ &
1122 & vh(ie,js,nint(lmv(ie,js)))+vh(iw,js,nint(lmv(iw,js))))/4.0
1123 END IF
1124 END IF
1125 230 CONTINUE
1126!
1127!
1128!---------------------------------------------------------------------
1129! *** PART II ***
1130!---------------------------------------------------------------------
1131!
1132! OUTPUT SELECTED FIELDS.
1133!
1134!---------------------------------------------------------------------
1135!
1136!
1137!--- Wind Shear (wind speed difference in knots between sfc and 2000 ft)
1138
1139 DO j=jsta,jend
1140 DO i=ista,iend
1141 IF(abs(uagl(i,j)-spval)>small .AND. &
1142 abs(vagl(i,j)-spval)>small)THEN
1143 IF(gridtype=='B' .OR. gridtype=='E')THEN
1144 grid1(i,j)=sqrt((uagl(i,j)-u10h(i,j))**2+ &
1145 (vagl(i,j)-v10h(i,j))**2)*1.943*zagl2(lp)/ &
1146 (zagl2(lp)-10.)
1147 ELSE
1148 grid1(i,j)=sqrt((uagl(i,j)-u10(i,j))**2+ &
1149 (vagl(i,j)-v10(i,j))**2)*1.943*zagl2(lp)/ &
1150 (zagl2(lp)-10.)
1151 END IF
1152 ELSE
1153 grid1(i,j)=spval
1154 END IF
1155 ENDDO
1156 ENDDO
1157 if(grib=="grib2" )then
1158 cfld=cfld+1
1159 fld_info(cfld)%ifld=iavblfld(iget(259))
1160 fld_info(cfld)%lvl=lvlsxml(lp,iget(259))
1161 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1162 endif
1163!
1164 ENDIF ! FOR LEVEL
1165!
1166!*** END OF MAIN VERTICAL LOOP
1167!
1168 320 CONTINUE
1169!*** ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
1170!
1171 ENDIF
1172! CRA
1173 IF (iget(411)>0 .OR. iget(412)>0 .OR. iget(413)>0) THEN
1174!
1175!---------------------------------------------------------------------
1176!***
1177!*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL
1178!*** INTERPOLATION ABOVE GROUND NOW.
1179!***
1180!
1181 DO 330 lp=1,lagl2
1182 iget1 = -1 ; iget2 = -1 ; iget3 = -1
1183 if (iget(411) > 0) iget1 = lvls(lp,iget(411))
1184 if (iget(412) > 0) iget2 = lvls(lp,iget(412))
1185 if (iget(413) > 0) iget3 = lvls(lp,iget(413))
1186 IF (iget1 > 0 .or. iget2 > 0 .or. iget3 > 0) then
1187
1188!
1189 jj = float(jsta+jend)/2.0
1190 ii = float(ista+iend)/3.0
1191 DO j=jsta_2l,jend_2u
1192 DO i=ista_2l,iend_2u
1193!
1194 pagl(i,j) = spval
1195 tagl(i,j) = spval
1196 qagl(i,j) = spval
1197 uagl(i,j) = spval
1198 vagl(i,j) = spval
1199!
1200!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
1201!*** THE AGL LEVEL TO WHICH WE ARE INTERPOLATING.
1202!
1203 llmh = nint(lmh(i,j))
1204 nl1x(i,j) = llmh+1
1205 DO l=llmh,2,-1
1206 zdum = zmid(i,j,l)-zint(i,j,llmh+1)
1207 IF(zdum >= zagl3(lp))THEN
1208 nl1x(i,j) = l+1
1209 exit
1210 ENDIF
1211 ENDDO
1212!
1213! IF THE AGL LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
1214! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
1215! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
1216! WILL EXTRAPOLATE TO THAT POINT
1217!
1218 IF(nl1x(i,j)==(llmh+1) .AND. zagl3(lp)>0.)THEN
1219 nl1x(i,j) = lm
1220 ENDIF
1221!
1222! if(NL1X(I,J)==LMP1)print*,'Debug: NL1X=LMP1 AT '
1223! 1 ,i,j,lp
1224 ENDDO
1225 ENDDO
1226!
1227!mptest IF(NHOLD==0)GO TO 310
1228!
1229!!$omp parallel do
1230!!$omp& private(nn,i,j,ll,fact,qsat,rhl)
1231!chc DO 220 NN=1,NHOLD
1232!chc I=IHOLD(NN)
1233!chc J=JHOLD(NN)
1234! DO 220 J=JSTA,JEND
1235 DO 240 j=jsta_2l,jend_2u
1236 DO 240 i=ista_2l,iend_2u
1237 ll = nl1x(i,j)
1238!---------------------------------------------------------------------
1239!*** VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
1240!*** HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
1241!---------------------------------------------------------------------
1242!
1243!CHC IF(NL1X(I,J)<=LM)THEN
1244 llmh = nint(lmh(i,j))
1245 IF(nl1x(i,j)<=llmh)THEN
1246!
1247!---------------------------------------------------------------------
1248! INTERPOLATE LINEARLY IN LOG(P)
1249!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
1250!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
1251!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
1252!---------------------------------------------------------------------
1253!
1254! FACT=(ALSL(LP)-ALOG(PMID(I,J,LL)))/
1255! & (ALOG(PMID(I,J,LL))-ALOG(PMID(I,J,LL-1)))
1256 zdum=zagl3(lp)+zint(i,j,nint(lmh(i,j))+1)
1257 fact = (zdum-zmid(i,j,ll)) &
1258 / (zmid(i,j,ll)-zmid(i,j,ll-1))
1259!
1260 paglu = log(pmid(i,j,ll-1))
1261 pagll = log(pmid(i,j,ll))
1262
1263 taglu = t(i,j,ll-1)
1264 tagll = t(i,j,ll)
1265
1266 qaglu = q(i,j,ll-1)
1267 qagll = q(i,j,ll)
1268
1269 uaglu = uh(i,j,ll-1)
1270 uagll = uh(i,j,ll)
1271
1272 vaglu = vh(i,j,ll-1)
1273 vagll = vh(i,j,ll)
1274
1275 pagl(i,j) = exp(pagll+(pagll-paglu)*fact)
1276 tagl(i,j) = tagll+(tagll-taglu)*fact
1277 qagl(i,j) = qagll+(qagll-qaglu)*fact
1278 uagl(i,j) = uagll+(uagll-uaglu)*fact
1279 vagl(i,j) = vagll+(vagll-vaglu)*fact
1280!
1281! FOR UNDERGROUND AGL LEVELS, ASSUME TEMPERATURE TO CHANGE
1282! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
1283! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
1284! GOUND
1285 ELSE
1286 pagl(i,j) = pmid(i,j,nint(lmv(i,j)))
1287 tagl(i,j) = t(i,j,nint(lmv(i,j)))
1288 qagl(i,j) = q(i,j,nint(lmv(i,j)))
1289 uagl(i,j) = uh(i,j,nint(lmv(i,j)))
1290 vagl(i,j) = vh(i,j,nint(lmv(i,j)))
1291 END IF
1292 240 CONTINUE
1293!
1294!
1295!---------------------------------------------------------------------
1296! *** PART II ***
1297!---------------------------------------------------------------------
1298!
1299! OUTPUT SELECTED FIELDS.
1300!
1301!---------------------------------------------------------------------
1302!
1303!
1304!--- Wind Energy Potential -- 0.5 * moist air density * wind speed^3
1305 IF((iget(411)>0) ) THEN
1306 DO j=jsta,jend
1307 DO i=ista,iend
1308 IF(qagl(i,j)<spval.and.pagl(i,j)<spval.and.tagl(i,j)<spval.and.&
1309 uagl(i,j)<spval.and.vagl(i,j)<spval)THEN
1310 qagl(i,j)=qagl(i,j)/1000.0
1311 pv=qagl(i,j)*pagl(i,j)/(eps*(1-qagl(i,j)) + qagl(i,j))
1312 rho=(1/tagl(i,j))*(((pagl(i,j)-pv)/rd) + pv/461.495)
1313 grid1(i,j)=0.5*rho*(sqrt(uagl(i,j)**2+vagl(i,j)**2))**3
1314 ELSE
1315 grid1(i,j)=spval
1316 ENDIF
1317 ENDDO
1318 ENDDO
1319 if(grib=="grib2" )then
1320 cfld=cfld+1
1321 fld_info(cfld)%ifld=iavblfld(iget(411))
1322 fld_info(cfld)%lvl=lvlsxml(lp,iget(411))
1323 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1324 endif
1325 ENDIF
1326!--- U Component of wind
1327 IF((iget(412)>0) ) THEN
1328 DO j=jsta,jend
1329 DO i=ista,iend
1330 grid1(i,j)=uagl(i,j)
1331 ENDDO
1332 ENDDO
1333 if(grib=="grib2" )then
1334 cfld=cfld+1
1335 fld_info(cfld)%ifld=iavblfld(iget(412))
1336 fld_info(cfld)%lvl=lvlsxml(lp,iget(412))
1337 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1338 endif
1339 ENDIF
1340!--- V Component of wind
1341 IF((iget(413)>0) ) THEN
1342 DO j=jsta,jend
1343 DO i=ista,iend
1344 grid1(i,j)=vagl(i,j)
1345 ENDDO
1346 ENDDO
1347 if(grib=="grib2" )then
1348 cfld=cfld+1
1349 fld_info(cfld)%ifld=iavblfld(iget(413))
1350 fld_info(cfld)%lvl=lvlsxml(lp,iget(413))
1351 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1352 endif
1353 ENDIF
1354!
1355 ENDIF ! FOR LEVEL
1356
1357!
1358!*** END OF MAIN VERTICAL LOOP
1359!
1360 330 CONTINUE
1361!*** ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
1362!
1363 ENDIF
1364! CRA
1365!
1366! END OF ROUTINE.
1367!
1368 RETURN
1369 END
1370
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27
mdl2agl
subroutine mdl2agl
SUBPROGRAM: MDL2P VERT INTRP OF MODEL LVLS TO AGL HEIGHT PRGRMMR: CHUANG ORG: W/NP22 DATE: 05-05-23
Definition MDL2AGL.f:51