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