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