UPP  001
 All Data Structures Files Functions Pages
INITPOST.F
Go to the documentation of this file.
1 
5 
24  SUBROUTINE initpost
25 
26  use vrbls4d, only: dust, smoke
27  use vrbls3d, only: t, u, uh, v, vh, wh, q, pmid, t, omga, pint, alpint, &
28  qqr, qqs, qqi, qqg, qqnw, qqni,qqnr, cwm, qqw, qqi, qqr, qqs, extcof55,&
29  f_ice, f_rain, f_rimef, q2, zint, zmid, ttnd, cfr, cfr_raw, qc_bl, ref_10cm, &
30  qqnwfa,qqnifa,taod5503d,aextc55
31  use vrbls2d, only: tmax, qrmax, htop, hbot, cuppt, fis, cfrach, cfracl, &
32  sr, cfrach, cfracm, wspd10max, w_up_max, w_dn_max, w_mean, refd_max, &
33  up_heli_max, up_heli_max16, grpl_max, up_heli, up_heli16, &
34  up_heli_min,up_heli_min16,up_heli_max02,up_heli_min02, &
35  up_heli_max03,up_heli_min03,rel_vort_max,rel_vort_max01, &
36  wspd10umax,wspd10vmax,refdm10c_max, &
37  hail_max2d,hail_maxk1,hail_maxhailcast,ltg1_max, &
38  ltg2_max, ltg3_max, nci_ltg, nca_ltg, nci_wq, nca_wq, nci_refd, &
39  u10, v10, th10, q10, tshltr, mrshltr, &
40  nca_refd, qv2m, qshltr, smstav, smstot, ssroff, bgroff, sfcevp, &
41  sfcexc, vegfrc, acsnow, cmc, sst, thz0, qz0, uz0, vz0, qs, qvg, &
42  z0, ustar, akhs, akms, radot, ths, acsnom, cuprec, ancprc, acprec, &
43  rainc_bucket, pcp_bucket, cprate, prec, snownc, snow_bucket, &
44  graup_bucket, swddni, swddif, mean_frp, acgraup, acfrain, &
45  graupelnc, albedo, rswin, rswout, swdnbc, swddnic, &
46  swddifc, swupbc, swupt, czen, czmean, rlwin, lwdnbc, lwupbc, &
47  rainnc_bucket, taod5502d, aerasy2d, aerssa2d, lwp, iwp, &
48  sigt4, rlwtoa, rswinc, aswin, aswout, alwin, alwout, alwtoa, aswtoa, &
49  tg, soiltb, twbs, qwbs,grnflx, sfcshx, sfclhx, subshx, snopcx, &
50  sfcuvx, potevp, ncfrcv, ncfrst, sno, si, pctsno, snonc, tsnow, &
51  ivgtyp, isltyp, islope, pblh, pblhgust, f, &
52  qvl1,refc_10cm,ref1km_10cm,ref4km_10cm, &
53  swradmean,u10mean,v10mean,spduv10mean,swnormmean,snfden,sndepac, &
54  hbotd,hbots,rainc_bucket1,rainnc_bucket1,pcp_bucket1,snow_bucket1, &
55  graup_bucket1, shdmin, shdmax, lai, htopd,htops
56  use soil, only: smc, sh2o, stc, sldpth, sllevel
57  use masks, only: lmv, lmh, vtm, sice, gdlat, gdlon, sm, dx, dy, htm
58  use ctlblk_mod, only: jsta_2l, jend_2u, filename, datahandle, datestr, &
59  ihrst, imin, idat, sdat, ifhr, ifmin, imp_physics, jsta, jend, &
60  spval,gdsdegr, modelname, pt, icu_physics, jsta_m, jend_m, nsoil, &
61  isf_surface_physics, nsoil, ardlw, ardsw, asrfc, me, mpi_comm_comp, &
62  nphs, smflag, spl, lsm, dt, prec_acc_dt, dtq2, tsrfc, trdlw, &
63  trdsw, theat, tclod, tprec, nprec, alsl, im, jm, lm, grib, &
64  prec_acc_dt1, submodelname
65  use params_mod, only: capa, g, rd, d608, tfrz, ad05, cft0, stbol, &
66  p1000, pi, rtd, lheat, dtr, erad
67  use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, &
68  qs0, sqs, sthe, the0, ttblq, rdpq, rdtheq, stheq, the0q
69  use gridspec_mod, only: gridtype, dxval, latstart, latlast, lonstart, &
70  lonlast, dyval, cenlat, cenlon, maptype, truelat1, truelat2, &
71  standlon, psmapf
72  use wrf_io_flags_mod, only:
73 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74  implicit none
75 !
76 ! INCLUDE/SET PARAMETERS.
77 !
78  include "mpif.h"
79 !
80 ! This version of INITPOST shows how to initialize, open, read from, and
81 ! close a NetCDF dataset. In order to change it to read an internal (binary)
82 ! dataset, do a global replacement of _ncd_ with _int_.
83 
84  character(len=31) :: varname
85  integer :: status
86  character startdate*19,sysdepinfo*80
87 !
88 ! NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
89 ! AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
90 !
91 ! ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
92 ! INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
93 
94  INTEGER idate(8),jdate(8)
95 !
96 ! DECLARE VARIABLES.
97 !
98  REAL rinc(5)
99  REAL dummy ( im, jm )
100  REAL dummy2 ( im, jm )
101  real, allocatable:: msft(:,:)
102  INTEGER idummy ( im, jm )
103  REAL,allocatable :: dum3d ( :, :, : )
104 ! REAL DUM3D2 ( IM+1, JM+1, LM+1 )
105  real, allocatable:: pvapor(:,:)
106  real, allocatable:: pvapor_orig(:,:)
107  REAL, ALLOCATABLE :: thv(:,:,:)
108 
109  integer js,je,jev,iyear,imn,iday,itmp,ioutcount,istatus, &
110  ii,jj,ll,i,j,l,nrdlw,nrdsw,n,igdout,irtn,idyvald, &
111  idxvald,nsrfc , lflip, k, k1
112  real dz,tsph,tmp,qmean,pvapornew,dumcst,tlmh,rho,zsf,zpbltop
113  real t2,th2,x2m,p2m,tsk, fact, temp
114  real lat
115 
116  integer jdn, numr, ic, jc, ierr
117  integer, external :: iw3jdn
118  real sun_zenith,sun_azimuth, ptop_low, ptop_mid, ptop_high
119  real delta_theta4gust
120 !
121 !
122 !***********************************************************************
123 ! START INIT HERE.
124 !
125  ALLOCATE ( thv(im,jsta_2l:jend_2u,lm) )
126  ALLOCATE ( dum3d( im+1, jm+1, lm+1 ) )
127  WRITE(6,*)'INITPOST: ENTER INITPOST'
128 !
129  gridtype='A'
130  hbotd=0
131  hbots=0
132  htopd=0
133  htops=0
134  ttnd=0
135  qs=0
136 ! STEP 1. READ MODEL OUTPUT FILE
137 !
138 !
139 !***
140 !
141 ! LMH always = LM for sigma-type vert coord
142 ! LMV always = LM for sigma-type vert coord
143 
144  do j = jsta_2l, jend_2u
145  do i = 1, im
146  lmv( i, j ) = lm
147  lmh( i, j ) = lm
148  end do
149  end do
150 
151 
152 ! HTM VTM all 1 for sigma-type vert coord
153 
154  do l = 1, lm
155  do j = jsta_2l, jend_2u
156  do i = 1, im
157  htm( i, j, l ) = 1.0
158  vtm( i, j, l ) = 1.0
159  end do
160  end do
161  end do
162 !
163 ! how do I get the filename?
164 ! fileName = '/ptmp/wx20mb/wrfout_01_030500'
165 ! DateStr = '2002-03-05_18:00:00'
166 ! how do I get the filename?
167  call ext_ncd_ioinit(sysdepinfo,status)
168  print*,'called ioinit', status
169  call ext_ncd_open_for_read( trim(filename), 0, 0, " ", &
170  datahandle, status)
171  print*,'called open for read', status
172  if ( status /= 0 ) then
173  print*,'error opening ',filename, ' Status = ', status ; stop
174  endif
175 ! get date/time info
176 ! this routine will get the next time from the file, not using it
177  print *,'DateStr before calling ext_ncd_get_next_time=',datestr
178 ! call ext_ncd_get_next_time(DataHandle, DateStr, Status)
179  print *,'DateStri,Status,DataHandle = ',datestr,status,datahandle
180 
181 ! The end j row is going to be jend_2u for all variables except for V.
182  js=jsta_2l
183  je=jend_2u
184  IF (jend_2u==jm) THEN
185  jev=jend_2u+1
186  ELSE
187  jev=jend_2u
188  ENDIF
189 !
190 ! Getting start time
191 #ifdef COMMCODE
192  call ext_ncd_get_dom_ti_char(datahandle,'SIMULATION_START_DATE', &
193  startdate, status )
194 #else
195  call ext_ncd_get_dom_ti_char(datahandle,'START_DATE',startdate, &
196  status )
197 #endif
198  print*,'startdate= ',startdate
199  jdate=0
200  idate=0
201  read(startdate,15)iyear,imn,iday,ihrst,imin
202  15 format(i4,1x,i2,1x,i2,1x,i2,1x,i2)
203  print*,'start yr mo day hr min=',iyear,imn,iday,ihrst,imin
204  print*,'processing yr mo day hr min=' &
205  ,idat(3),idat(1),idat(2),idat(4),idat(5)
206  idate(1)=iyear
207  idate(2)=imn
208  idate(3)=iday
209  idate(5)=ihrst
210  idate(6)=imin
211  sdat(1)=imn
212  sdat(2)=iday
213  sdat(3)=iyear
214  jdate(1)=idat(3)
215  jdate(2)=idat(1)
216  jdate(3)=idat(2)
217  jdate(5)=idat(4)
218  jdate(6)=idat(5)
219 ! CALL W3DIFDAT(JDATE,IDATE,2,RINC)
220 ! ifhr=nint(rinc(2))
221  CALL w3difdat(jdate,idate,0,rinc)
222  ifhr=nint(rinc(2)+rinc(1)*24.)
223  ifmin=nint(rinc(3))
224  print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,filename
225 ! OK, since all of the variables are dimensioned/allocated to be
226 ! the same size, this means we have to be careful int getVariable
227 ! to not try to get too much data. For example,
228 ! DUM3D is dimensioned IM+1,JM+1,LM+1 but there might actually
229 ! only be im,jm,lm points of data available for a particular variable.
230 
231  call ext_ncd_get_dom_ti_integer(datahandle,'MP_PHYSICS' &
232  ,itmp,1,ioutcount,istatus)
233  imp_physics=itmp
234  print*,'MP_PHYSICS= ',itmp
235 
236 ! Initializes constants for Ferrier microphysics
237  if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
238  CALL microinit(imp_physics)
239  end if
240 
241  call ext_ncd_get_dom_ti_integer(datahandle,'CU_PHYSICS' &
242  ,itmp,1,ioutcount,istatus)
243  icu_physics=itmp
244  print*,'CU_PHYSICS= ',icu_physics
245 
246 ! get 3-D variables
247  print*,'im,jm,lm= ',im,jm,lm
248  ii=im/2
249  jj=(jsta+jend)/2
250  ll=lm
251 
252  varname='T'
253  call getvariable(filename,datestr,datahandle,varname,dum3d, &
254  im+1,1,jm+1,lm+1,im,js,je,lm)
255  do l = 1, lm
256  do j = jsta_2l, jend_2u
257  do i = 1, im
258  t( i, j, l ) = dum3d( i, j, l ) + 300.
259 !MEB this is theta the 300 is my guess at what T0 is
260  end do
261  end do
262  end do
263 
264  do l=1,lm
265  end do
266  varname='U'
267  call getvariable(filename,datestr,datahandle,varname,dum3d, &
268  im+1,1,jm+1,lm+1,im+1,js,je,lm)
269  do l = 1, lm
270  do j = jsta_2l, jend_2u
271  do i = 1, im+1
272  u( i, j, l ) = dum3d( i, j, l )
273  end do
274  end do
275 ! fill up UH which is U at P-points including 2 row halo
276  do j = jsta_2l, jend_2u
277  do i = 1, im
278  uh(i,j,l) = (dum3d(i,j,l)+dum3d(i+1,j,l))*0.5
279  end do
280  end do
281  end do
282 ! if(jj>= jsta .and. jj<=jend)print*,'sample U= ',U(ii,jj,ll)
283  varname='V'
284  call getvariable(filename,datestr,datahandle,varname,dum3d, &
285  im+1,1,jm+1,lm+1,im, js,jev,lm)
286  do l = 1, lm
287  do j = jsta_2l, jev
288  do i = 1, im
289  v( i, j, l ) = dum3d( i, j, l )
290  end do
291  end do
292 ! fill up VH which is V at P-points including 2 row halo
293  do j = jsta_2l, jend_2u
294  do i = 1, im
295  vh(i,j,l) = (dum3d(i,j,l)+dum3d(i,j+1,l))*0.5
296  end do
297  end do
298  end do
299 ! if(jj>= jsta .and. jj<=jend)print*,'sample V= ',V(ii,jj,ll)
300 
301  varname='W'
302  call getvariable(filename,datestr,datahandle,varname,dum3d, &
303  im+1,1,jm+1,lm+1,im, js,je,lm+1)
304 ! do l = 1, lm+1
305 ! do j = jsta_2l, jend_2u
306 ! do i = 1, im
307 ! w ( i, j, l ) = dum3d ( i, j, l )
308 ! end do
309 ! end do
310 ! end do
311 ! fill up WH which is W at P-points including 2 row halo
312  DO l=1,lm
313  DO i=1,im
314  DO j=jsta_2l,jend_2u
315  wh(i,j,l) = (dum3d(i,j,l)+dum3d(i,j,l+1))*0.5
316  ENDDO
317  ENDDO
318  ENDDO
319  print*,'finish reading W'
320 
321  varname='QVAPOR'
322  call getvariable(filename,datestr,datahandle,varname,dum3d, &
323  im+1,1,jm+1,lm+1,im,js,je,lm)
324  do l = 1, lm
325  do j = jsta_2l, jend_2u
326  do i = 1, im
327 !HC q ( i, j, l ) = dum3d ( i, j, l )
328 !HC CONVERT MIXING RATIO TO SPECIFIC HUMIDITY
329 !mhu check !!!! if (dum3d(i,j,l) < 10E-12) dum3d(i,j,l) = 10E-12
330  q( i, j, l ) = dum3d( i, j, l )/(1.0+dum3d( i, j, l ))
331  end do
332  end do
333  end do
334  print*,'finish reading mixing ratio'
335 ! if(jj>= jsta .and. jj<=jend)print*,'sample Q= ',Q(ii,jj,ll)
336 
337 ! DCD 4/3/13
338 ! previously initialized PMID from sum of base-state (PB) and
339 ! perturbation (P) pressure now initializing PMID from
340 ! hydrostatic pressure (P_HYD), motivated particularly by 0h HRRR
341 ! VarName='PB'
342  varname='P_HYD'
343  call getvariable(filename,datestr,datahandle,varname,dum3d, &
344  im+1,1,jm+1,lm+1,im, js,je,lm)
345 ! VarName='P'
346 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D2, &
347 ! IM+1,1,JM+1,LM+1,IM, JS,JE,LM)
348  do l = 1, lm
349  do j = jsta_2l, jend_2u
350  do i = 1, im
351 ! PMID(I,J,L)=DUM3D(I,J,L)+DUM3D2(I,J,L)
352  pmid(i,j,l)=dum3d(i,j,l)
353  thv( i, j, l ) = t(i,j,l)*(q(i,j,l)*0.608+1.)
354 ! now that I have P, convert theta to t
355  t( i, j, l ) = t(i,j,l)*(pmid(i,j,l)*1.e-5)**capa
356 ! now that I have T,q,P compute omega from wh
357  if(abs(t( i, j, l ))>1.0e-3) &
358  omga(i,j,l) = -wh(i,j,l)*pmid(i,j,l)*g/ &
359  (rd*t(i,j,l)*(1.+d608*q(i,j,l)))
360 ! seperate rain from snow and cloud water from cloud ice for WSM3 scheme
361 ! if(imp_physics == 3)then
362 ! if(t(i,j,l) < TFRZ)then
363 ! qqs(i,j,l)=qqr(i,j,l)
364 ! qqi(i,j,l)=qqw(i,j,l)
365 ! end if
366 ! end if
367  end do
368  end do
369  end do
370 !tgs - 6 June 2012 - added check on monotonic PMID
371  do l = 2, lm-1
372  ll=lm-l+1
373  do j = jsta_2l, jend_2u
374  do i = 1, im
375  if((pmid(i,j,ll-1) - pmid(i,j,ll))>=0.) then
376  write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll
377  write(*,*) 'PMID: ll-1,ll,ll+1', pmid(i,j,ll-1) &
378  ,pmid(i,j,ll), pmid(i,j,ll+1)
379 
380  pmid(i,j,ll)=0.5*(pmid(i,j,ll+1)+pmid(i,j,ll-1))
381 
382  write(*,*) 'after adjustment-i,j,ll ', i,j,ll
383  write(*,*) 'PMID: ll-1,ll,ll+1', pmid(i,j,ll-1) &
384  ,pmid(i,j,ll), pmid(i,j,ll+1)
385  endif
386  end do
387  end do
388  end do
389 !
390 ! check the lowest level too: set l=1, then extroplate the lowest level
391 ! based on znu
392 ! znu(lm)=0.9990000, znu(lm-1)=0.9960001, znu(lm-2)=0.9905000
393 ! P(lm)=p(lm-2) + (p(lm-1)-p(lm-2))*(znu(lm)-znu(lm-2))/ &
394 ! (znu(lm-1)-znu(lm-2))
395 ! where: (znu(lm)-znu(lm-2))/(znu(lm-1)-znu(lm-2))=17/11
396 ! Thus: P(lm)=p(lm-2) + (p(lm-1)-p(lm-2))*17/11
397 ! =p(lm-1) + (p(lm-1)-p(lm-2))*6/11
398  fact=6.0/11.0
399  ll=lm
400  do j = jsta_2l, jend_2u
401  do i = 1, im
402  if((pmid(i,j,ll-1) - pmid(i,j,ll))>=0.) then
403  write(*,*) 'non-monotonic PMID, i,j,ll ', i,j,ll
404  write(*,*) 'PMID: ll-2,ll-1,ll', pmid(i,j,ll-2) &
405  ,pmid(i,j,ll-1), pmid(i,j,ll)
406 
407  pmid(i,j,ll)=pmid(i,j,ll-1) + &
408  fact*(pmid(i,j,ll-1)-pmid(i,j,ll-2))
409 
410  write(*,*) 'after adjustment-i,j,ll ', i,j,ll
411  write(*,*) 'PMID: ll-2,ll-1,ll', pmid(i,j,ll-2) &
412  ,pmid(i,j,ll-1), pmid(i,j,ll)
413  endif
414  end do
415  end do
416 
417 ! end check
418 
419  DO l=2,lm
420  DO i=1,im
421  DO j=jsta_2l,jend_2u
422  pint(i,j,l)=(pmid(i,j,l-1)+pmid(i,j,l))*0.5
423  alpint(i,j,l)=alog(pint(i,j,l))
424  ENDDO
425  ENDDO
426  ENDDO
427 !--- Compute max temperature in the column up to level 20
428 !--- to be used later in precip type computation
429  do j = jsta_2l, jend_2u
430  do i = 1, im
431  tmax(i,j)=0.
432  end do
433  end do
434 
435  do l = 2,20
436  lflip = lm - l + 1
437  do j = jsta_2l, jend_2u
438  do i = 1, im
439  tmax(i,j)=max(tmax(i,j),t(i,j,lflip))
440  end do
441  end do
442  end do
443 
444 
445 ! Brad comment out the output of individual species for Ferriers scheme within
446 ! ARW in Registry file
447 
448  qqw=0.
449  qqr=0.
450  qqs=0.
451  qqi=0.
452  qqg=0.
453  qqni=0.
454  qqnr=0.
455  cwm=0.
456 
457 ! extinction coef for aerosol
458  extcof55=0.
459  aextc55=0.
460 
461  if(imp_physics/=5 .and. imp_physics/=0)then
462  varname='QCLOUD'
463  call getvariable(filename,datestr,datahandle,varname,dum3d, &
464  im+1,1,jm+1,lm+1,im, js,je,lm)
465  do l = 1, lm
466  do j = jsta_2l, jend_2u
467  do i = 1, im
468 ! partition cloud water and ice for WSM3
469  if(imp_physics==3)then
470  if(t(i,j,l) >= tfrz)then
471  qqw( i, j, l ) = dum3d( i, j, l )
472  else
473  qqi( i, j, l ) = dum3d( i, j, l )
474  end if
475  else ! bug fix provided by J CASE
476  qqw( i, j, l ) = dum3d( i, j, l )
477  end if
478  end do
479  end do
480  end do
481 ! if(jj>= jsta .and. jj<=jend)
482 ! + print*,'sample QCLOUD= ',QQW(ii,jj,ll)
483 ! print*,'finish reading cloud mixing ratio'
484  end if
485 
486 
487 
488  if(imp_physics/=1 .and. imp_physics/=3 &
489  .and. imp_physics/=5 .and. imp_physics/=0)then
490  varname='QICE'
491  call getvariable(filename,datestr,datahandle,varname,dum3d, &
492  im+1,1,jm+1,lm+1,im, js,je,lm)
493  do l = 1, lm
494  do j = jsta_2l, jend_2u
495  do i = 1, im
496  qqi( i, j, l ) = dum3d( i, j, l )
497  end do
498  end do
499  end do
500 ! if(jj>= jsta .and. jj<=jend)
501 ! + print*,'sample QICE= ',qqi(ii,jj,ll)
502  end if
503 
504 
505  if(imp_physics/=5 .and. imp_physics/=0)then
506  varname='QRAIN'
507  call getvariable(filename,datestr,datahandle,varname,dum3d, &
508  im+1,1,jm+1,lm+1,im, js,je,lm)
509  do l = 1, lm
510  do j = jsta_2l, jend_2u
511  do i = 1, im
512 ! partition rain and snow for WSM3
513  if(imp_physics == 3)then
514  if(t(i,j,l) >= tfrz)then
515  qqr( i, j, l ) = dum3d( i, j, l )
516  else
517  qqs( i, j, l ) = dum3d( i, j, l )
518  end if
519  else ! bug fix provided by J CASE
520  qqr( i, j, l ) = dum3d( i, j, l )
521  end if
522  dummy(i,j)=dum3d(i,j,l)
523  end do
524  end do
525 ! print*,'max rain water= ',l,maxval(dummy)
526  end do
527 ! if(jj>= jsta .and. jj<=jend)
528 ! + print*,'sample QRAIN= ',qqr(ii,jj,ll)
529 !tgs
530 ! Compute max QRAIN in the column to be used later in precip type computation
531  do j = jsta_2l, jend_2u
532  do i = 1, im
533  qrmax(i,j)=0.
534  end do
535  end do
536 
537  do l = 1, lm
538  do j = jsta_2l, jend_2u
539  do i = 1, im
540  qrmax(i,j)=max(qrmax(i,j),qqr(i,j,l))
541  end do
542  end do
543  end do
544 
545  end if
546 
547 
548  if(imp_physics/=1 .and. imp_physics/=3 .and. &
549  imp_physics/=5 .and. imp_physics/=0)then
550  varname='QSNOW'
551  call getvariable(filename,datestr,datahandle,varname,dum3d, &
552  im+1,1,jm+1,lm+1,im, js,je,lm)
553  do l = 1, lm
554  do j = jsta_2l, jend_2u
555  do i = 1, im
556  qqs( i, j, l ) = dum3d( i, j, l )
557  dummy(i,j)=dum3d(i,j,l)
558  end do
559  end do
560 ! print*,'max snow= ',l,maxval(dummy)
561  end do
562  end if
563 
564 
565  if(imp_physics==2 .or. imp_physics==6 .or. &
566  imp_physics==8 .or. imp_physics==9 .or. imp_physics==28)then
567  varname='QGRAUP'
568  call getvariable(filename,datestr,datahandle,varname,dum3d, &
569  im+1,1,jm+1,lm+1,im, js,je,lm)
570  do l = 1, lm
571  do j = jsta_2l, jend_2u
572  do i = 1, im
573  qqg( i, j, l ) = dum3d( i, j, l )
574  end do
575  end do
576  end do
577 ! if(jj>= jsta .and. jj<=jend)
578 ! + print*,'sample QGRAUP= ',qqg(ii,jj,ll)
579  end if
580 
581  if(imp_physics==8 .or. imp_physics==9 .or.imp_physics==28)then
582  varname='QNICE'
583  call getvariable(filename,datestr,datahandle,varname,dum3d, &
584  im+1,1,jm+1,lm+1,im, js,je,lm)
585  do l = 1, lm
586  do j = jsta_2l, jend_2u
587  do i = 1, im
588  qqni( i, j, l ) = dum3d( i, j, l )
589  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample QQNI= ', &
590  i,j,l,qqni( i, j, l )
591  end do
592  end do
593  end do
594  varname='QNRAIN'
595  call getvariable(filename,datestr,datahandle,varname,dum3d, &
596  im+1,1,jm+1,lm+1,im, js,je,lm)
597  do l = 1, lm
598  do j = jsta_2l, jend_2u
599  do i = 1, im
600  qqnr( i, j, l ) = dum3d( i, j, l )
601  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample QQNR= ', &
602  i,j,l,qqnr( i, j, l )
603  end do
604  end do
605  end do
606  end if
607 
608 ! For aerosol aware microphyscis
609  if(imp_physics==28) then
610  varname='QNCLOUD'
611  call getvariable(filename,datestr,datahandle,varname,dum3d, &
612  im+1,1,jm+1,lm+1,im, js,je,lm)
613  do l = 1, lm
614  do j = jsta_2l, jend_2u
615  do i = 1, im
616  qqnw( i, j, l ) = dum3d( i, j, l )
617  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample QQNW= ', &
618  i,j,l,qqnw( i, j, l )
619  end do
620  end do
621  end do
622  varname='QNWFA'
623  call getvariable(filename,datestr,datahandle,varname,dum3d, &
624  im+1,1,jm+1,lm+1,im, js,je,lm)
625  do l = 1, lm
626  do j = jsta_2l, jend_2u
627  do i = 1, im
628  qqnwfa( i, j, l ) = dum3d( i, j, l )
629  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample QQNWFA= ', &
630  i,j,l,qqnwfa( i, j, l )
631  end do
632  end do
633  end do
634  varname='QNIFA'
635  call getvariable(filename,datestr,datahandle,varname,dum3d, &
636  im+1,1,jm+1,lm+1,im, js,je,lm)
637  do l = 1, lm
638  do j = jsta_2l, jend_2u
639  do i = 1, im
640  qqnifa( i, j, l ) = dum3d( i, j, l )
641  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample QQNIFA= ', &
642  i,j,l,qqnifa( i, j, l )
643  end do
644  end do
645  end do
646  end if
647 
648 ! Read in extinction coefficient for aerosol at 550 nm
649 ! VarName='EXTCOF55'
650 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D, &
651 ! IM+1,1,JM+1,LM+1,IM,JS,JE,LM)
652 ! do l = 1, lm
653 ! do j = jsta_2l, jend_2u
654 ! do i = 1, im
655 ! extcof55 ( i, j, l ) = dum3d ( i, j, l )
656 ! end do
657 ! end do
658 ! end do
659 ! print*,'finish reading EXTCOF55'
660 
661  if(imp_physics/=5)then
662 !HC SUM UP ALL CONDENSATE FOR CWM
663  do l = 1, lm
664  do j = jsta_2l, jend_2u
665  do i = 1, im
666  IF(qqr(i,j,l)<spval)THEN
667  cwm(i,j,l)=qqr(i,j,l)
668  END IF
669  IF(qqi(i,j,l)<spval)THEN
670  cwm(i,j,l)=cwm(i,j,l)+qqi(i,j,l)
671  END IF
672  IF(qqw(i,j,l)<spval)THEN
673  cwm(i,j,l)=cwm(i,j,l)+qqw(i,j,l)
674  END IF
675  IF(qqs(i,j,l)<spval)THEN
676  cwm(i,j,l)=cwm(i,j,l)+qqs(i,j,l)
677  END IF
678  IF(qqg(i,j,l)<spval)THEN
679  cwm(i,j,l)=cwm(i,j,l)+qqg(i,j,l)
680  END IF
681  end do
682  end do
683  end do
684  else
685 
686  varname='CWM'
687  call getvariable(filename,datestr,datahandle,varname,dum3d, &
688  im+1,1,jm+1,lm+1,im, js,je,lm)
689  do l = 1, lm
690  do j = jsta_2l, jend_2u
691  do i = 1, im
692  cwm( i, j, l ) = dum3d( i, j, l )
693  end do
694  end do
695  end do
696 
697  varname='F_ICE_PHY'
698  call getvariable(filename,datestr,datahandle,varname,dum3d, &
699  im+1,1,jm+1,lm+1,im, js,je,lm)
700  do l = 1, lm
701  do j = jsta_2l, jend_2u
702  do i = 1, im
703  f_ice( i, j, l ) = dum3d( i, j, l )
704  end do
705  end do
706  end do
707 
708  varname='F_RAIN_PHY'
709  call getvariable(filename,datestr,datahandle,varname,dum3d, &
710  im+1,1,jm+1,lm+1,im, js,je,lm)
711  do l = 1, lm
712  do j = jsta_2l, jend_2u
713  do i = 1, im
714  f_rain( i, j, l ) = dum3d( i, j, l )
715  end do
716  end do
717  end do
718 
719  varname='F_RIMEF_PHY'
720  call getvariable(filename,datestr,datahandle,varname,dum3d, &
721  im+1,1,jm+1,lm+1,im, js,je,lm)
722  do l = 1, lm
723  do j = jsta_2l, jend_2u
724  do i = 1, im
725  f_rimef( i, j, l ) = dum3d( i, j, l )
726  end do
727  end do
728  end do
729 
730  end if
731 
732  varname='HTOP'
733  IF(icu_physics == 3 .or. icu_physics == 5) varname='CUTOP'
734  call getvariable(filename,datestr,datahandle,varname,dummy, &
735  im,1,jm,1,im,js,je,1)
736  do j = jsta_2l, jend_2u
737  do i = 1, im
738  htop( i, j ) = float(lm)-dummy(i,j)+1.0
739  end do
740  end do
741  varname='HBOT'
742  IF(icu_physics == 3 .or. icu_physics == 5) varname='CUBOT'
743  call getvariable(filename,datestr,datahandle,varname,dummy, &
744  im,1,jm,1,im,js,je,1)
745  do j = jsta_2l, jend_2u
746  do i = 1, im
747  hbot( i, j ) = float(lm)-dummy(i,j)+1.0
748  end do
749  end do
750 
751  varname='CUPPT'
752  call getvariable(filename,datestr,datahandle,varname,dummy, &
753  im,1,jm,1,im,js,je,1)
754  do j = jsta_2l, jend_2u
755  do i = 1, im
756  cuppt( i, j ) = dummy( i, j )
757  end do
758  end do
759 
760 
761  IF(modelname == 'RAPR')THEN
762  call getvariable(filename,datestr,datahandle,'QKE',dum3d, &
763  im+1,1,jm+1,lm+1,im,js,je,lm)
764  do l = 1, lm
765  do j = jsta_2l, jend_2u
766  do i = 1, im
767  q2( i, j, l ) = dum3d( i, j, l ) / 2.0
768  end do
769  end do
770  end do
771  ELSE
772  call getvariable(filename,datestr,datahandle,'TKE',dum3d, &
773  im+1,1,jm+1,lm+1,im,js,je,lm)
774  do l = 1, lm
775  do j = jsta_2l, jend_2u
776  do i = 1, im
777  q2( i, j, l ) = dum3d( i, j, l )
778  end do
779  end do
780  end do
781  ENDIF
782 
783 
784 !MEB call getVariable(fileName,DateStr,DataHandle,'QRAIN',new)
785 
786 ! get sfc pressure
787  varname='MU'
788  call getvariable(filename,datestr,datahandle,varname,dummy, &
789  im,1,jm,1,im,js,je,1)
790  varname='MUB'
791  call getvariable(filename,datestr,datahandle,varname,dummy2, &
792  im,1,jm,1,im,js,je,1)
793  varname='P_TOP'
794  call getvariable(filename,datestr,datahandle,varname,pt, &
795  1,1,1,1,1,1,1,1)
796 
797  DO i=1,im
798  DO j=js,je
799  pint(i,j,lm+1) = dummy(i,j)+dummy2(i,j)+pt
800  pint(i,j,1) = pt
801  alpint(i,j,lm+1)=alog(pint(i,j,lm+1))
802  alpint(i,j,1)=alog(pint(i,j,1))
803  ENDDO
804  ENDDO
805 
806 
807  varname='HGT'
808  call getvariable(filename,datestr,datahandle,varname,dummy, &
809  im,1,jm,1,im,js,je,1)
810  do j = jsta_2l, jend_2u
811  do i = 1, im
812  fis( i, j ) = dummy( i, j ) * g
813  end do
814  end do
815 !
816  varname='PHB'
817  call getvariable(filename,datestr,datahandle,varname,dum3d, &
818  im+1,1,jm+1,lm+1,im,js,je,lm+1)
819  DO l=1,lm+1
820  DO i=1,im
821  DO j=js,je
822  zint(i,j,l)=dum3d(i,j,l)
823  ENDDO
824  ENDDO
825  ENDDO
826  varname='PH'
827  call getvariable(filename,datestr,datahandle,varname,dum3d, &
828  im+1,1,jm+1,lm+1,im,js,je,lm+1)
829 
830  print*,'finish reading geopotential'
831 ! ph/phb are geopotential z=(ph+phb)/9.801
832  DO l=1,lm+1
833  DO i=1,im
834  DO j=js,je
835  zint(i,j,l)=(zint(i,j,l)+dum3d(i,j,l))/g
836  ENDDO
837  ENDDO
838  ENDDO
839 
840  IF(modelname == 'RAPR')THEN
841 
842  varname='PSFC'
843  call getvariable(filename,datestr,datahandle,varname,dummy, &
844  im,1,jm,1,im,js,je,1)
845 
846  DO j=jsta,jend
847  DO i=1,im
848  if((pint(i,j,lm) - dummy(i,j))>=0.) then
849  write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm
850  write(*,*) 'PINT: lm,lm+1, PMID: lm', pint(i,j,lm),dummy(i,j), pmid(i,j,lm)
851  dummy(i,j)=pmid(i,j,lm)*1.001
852  write(*,*) 'after adjustment-i,j,lm+1 ', i,j,lm+1,dummy(i,j)
853  endif
854  pint(i,j,lm+1)=dummy(i,j)
855  alpint(i,j,lm+1)=alog(pint(i,j,lm+1))
856  ENDDO
857  ENDDO
858 
859  ELSE
860 !!!!!!!!!!!!!
861 ! Pyle's and Chuang's fixes for ARW SLP
862 
863  allocate(pvapor(im,jsta_2l:jend_2u))
864  allocate(pvapor_orig(im,jsta_2l:jend_2u))
865  DO j=jsta,jend
866  DO i=1,im
867 
868 
869  pvapor(i,j)=0.
870  do l=1,lm
871  dz=zint(i,j,l)-zint(i,j,l+1)
872  rho=pmid(i,j,l)/(rd*t(i,j,l))
873 
874 
875  if (l <= lm-1) then
876  qmean=0.5*(q(i,j,l)+q(i,j,l+1))
877  else
878  qmean=q(i,j,l)
879  endif
880 
881 
882  pvapor(i,j)=pvapor(i,j)+g*rho*dz*qmean
883  enddo
884 
885 ! test elim
886 ! pvapor(I,J)=0.
887 
888 
889  pvapor_orig(i,j)=pvapor(i,j)
890 
891 
892  ENDDO
893  ENDDO
894 
895  do l=1,405
896  call exch(pvapor(1,jsta_2l))
897  do j=jsta_m,jend_m
898  do i=2,im-1
899 
900  pvapornew=ad05*(4.*(pvapor(i-1,j)+pvapor(i+1,j) &
901  +pvapor(i,j-1)+pvapor(i,j+1)) &
902  +pvapor(i-1,j-1)+pvapor(i+1,j-1) &
903  +pvapor(i-1,j+1)+pvapor(i+1,j+1)) &
904  -cft0*pvapor(i,j)
905 
906  pvapor(i,j)=pvapornew
907 
908  enddo
909  enddo
910  enddo ! iteration loop
911 
912 ! southern boundary
913  if (js == 1) then
914  j=1
915  do i=2,im-1
916  pvapor(i,j)=pvapor_orig(i,j)+(pvapor(i,j+1)-pvapor_orig(i,j+1))
917  enddo
918  endif
919 
920 ! northern boundary
921 
922  if (je == jm) then
923  j=jm
924  do i=2,im-1
925  pvapor(i,j)=pvapor_orig(i,j)+(pvapor(i,j-1)-pvapor_orig(i,j-1))
926  enddo
927  endif
928 
929 ! western boundary
930  i=1
931  do j=js,je
932  pvapor(i,j)=pvapor_orig(i,j)+(pvapor(i+1,j)-pvapor_orig(i+1,j))
933  enddo
934 
935 ! eastern boundary
936  i=im
937  do j=js,je
938  pvapor(i,j)=pvapor_orig(i,j)+(pvapor(i-1,j)-pvapor_orig(i-1,j))
939  enddo
940 
941  DO j=jsta,jend
942  DO i=1,im
943  pint(i,j,lm+1)=pint(i,j,lm+1)+pvapor(i,j)
944 ! KRF - check surface pressure for monotonic correctness
945 
946  if((pint(i,j,lm) - pint(i,j,lm+1))>=0. ) then
947  write(*,*) 'non-monotonic PINT, i,j,lm ', i,j,lm
948  write(*,*) 'PINT: lm,lm+1, PMID: lm', pint(i,j,lm), pint(i,j,lm+1), pmid(i,j,lm)
949  pint(i,j,lm+1) = pint(i,j,lm)*1.001
950  write(*,*) 'after adjustment-i,j,lm+1 PINT ', i,j,lm+1, pint(i,j,lm+1)
951  endif
952  alpint(i,j,lm+1)=alog(pint(i,j,lm+1))
953  ENDDO
954  ENDDO
955 
956  deallocate(pvapor)
957  deallocate(pvapor_orig)
958 
959  ENDIF ! IF(MODELNAME == 'RAPR')THEN
960 
961 !!!!!!!!!!!!!
962  IF(modelname == 'RAPR')THEN
963 !integrate heights hydrostatically
964  do j = js, je
965  do i = 1, im
966  zint(i,j,lm+1)=fis(i,j)/g
967  dummy(i,j)=fis(i,j)
968  if(i==im/2.and.j==(jsta+jend)/2) &
969  print*,'i,j,L,ZINT from unipost= ',i,j,lm+1,zint(i,j,lm+1) &
970  , alpint(i,j,lm+1),alpint(i,j,lm)
971  end do
972  end do
973  DO l=lm,1,-1
974  do j = js, je
975  do i = 1, im
976  dummy2(i,j)=htm(i,j,l)*t(i,j,l)*(q(i,j,l)*d608+1.0)*rd* &
977  (alpint(i,j,l+1)-alpint(i,j,l))+dummy(i,j)
978 ! compute difference between model and unipost heights:
979  dum3d(i,j,l)=zint(i,j,l)-dummy2(i,j)/g
980 ! now replace model heights with unipost heights
981  zint(i,j,l)=dummy2(i,j)/g
982  if(i==im/2.and.j==(jsta+jend)/2) &
983  print*,'i,j,L,ZINT from unipost= ',i,j,l,zint(i,j,l)
984  dummy(i,j)=dummy2(i,j)
985  ENDDO
986  ENDDO
987  END DO
988  DO l=lm,1,-1
989  do j = js, je
990  do i = 1, im
991  if(i==im/2.and.j==(jsta+jend)/2) then
992  print*,'DIFF heights model-unipost= ', &
993  i,j,l,dum3d(i,j,l)
994  endif
995  ENDDO
996  ENDDO
997  END DO
998 
999  print*,'finish deriving geopotential in ARW'
1000 
1001  ENDIF ! IF(MODELNAME == 'RAPR')THEN
1002 
1003 
1004  IF(modelname == 'RAPR')THEN
1005 
1006  DO l=1,lm-1
1007  DO i=1,im
1008  DO j=js,je
1009  fact=(alog(pmid(i,j,l))-alpint(i,j,l))/ &
1010  max(1.e-6,(alpint(i,j,l+1)-alpint(i,j,l)))
1011  zmid(i,j,l)=zint(i,j,l)+(zint(i,j,l+1)-zint(i,j,l))*fact
1012  dummy(i,j)=zmid(i,j,l)
1013  if((alpint(i,j,l+1)-alpint(i,j,l)) < 1.e-6) print*, &
1014  'P(K+1) and P(K) are too close, i,j,L,', &
1015  'ALPINT(I,J,L+1),ALPINT(I,J,L),ZMID = ', &
1016  i,j,l,alpint(i,j,l+1),alpint(i,j,l),zmid(i,j,l)
1017  ENDDO
1018  ENDDO
1019  print*,'max/min ZMID= ',l,maxval(dummy),minval(dummy)
1020  ENDDO
1021 
1022  DO i=1,im
1023  DO j=js,je
1024  DO l=1,lm
1025  zint(i,j,l+1) =amin1(zint(i,j,l)-2.,zint(i,j,l+1))
1026  zmid(i,j,l)=(zint(i,j,l+1)+zint(i,j,l))*0.5 ! ave of z
1027  ENDDO
1028  dummy(i,j)=zmid(i,j,lm)
1029  ENDDO
1030  ENDDO
1031  print*,'max/min ZMID= ',lm,maxval(zmid(1:im,js:je,lm)), &
1032  minval(zmid(1:im,js:je,lm))
1033  ELSE
1034 
1035  DO l=1,lm
1036  DO i=1,im
1037  DO j=js,je
1038  zmid(i,j,l)=(zint(i,j,l+1)+zint(i,j,l))*0.5 ! ave of z
1039 ! if(i==297.and.j==273) &
1040 ! print*,'i,j,L,ZMID = ', &
1041 ! i,j,l,ZMID(I,J,L)
1042  ENDDO
1043  ENDDO
1044  ENDDO
1045  ENDIF ! IF(MODELNAME == 'RAPR')THEN
1046 
1047 !
1048 ! E. James - 8 Dec 2017: this is for HRRR-smoke; it needs to be after ZINT
1049 ! is defined.
1050 !
1051  if(imp_physics==28) then
1052  varname='AOD3D_SMOKE'
1053  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1054  im+1,1,jm+1,lm+1,im, js,je,lm)
1055  do l = 1, lm
1056  do j = jsta_2l, jend_2u
1057  do i = 1, im
1058  taod5503d( i, j, l ) = dum3d( i, j, l )
1059  dz = zint( i, j, l ) - zint( i, j, l+1 )
1060  aextc55( i, j, l ) = taod5503d( i, j, l ) / dz
1061  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample TAOD5503D= ', &
1062  i,j,l,taod5503d( i, j, l )
1063  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample dz= ', &
1064  dz
1065  if(i==im/2.and.j==(jsta+jend)/2)print*,'sample AEXTC55= ', &
1066  i,j,l,aextc55( i, j, l )
1067  end do
1068  end do
1069  end do
1070  end if
1071 
1072 ! get 3-d soil variables
1073  varname='SMOIS'
1074  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1075  im+1,1,jm+1,lm+1,im,js,je,nsoil)
1076  do l = 1, nsoil
1077  do j = jsta_2l, jend_2u
1078  do i = 1, im
1079 ! smc ( i, j, l ) = dum3d ( i, j, l )
1080 ! flip soil layer again because wrf soil variable vertical indexing
1081 ! is the same with eta and vertical indexing was flipped for both
1082 ! atmospheric and soil layers within getVariable
1083  smc( i, j, l ) = dum3d( i, j, nsoil-l+1)
1084  end do
1085  end do
1086  end do
1087 
1088  varname='SH2O'
1089  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1090  im+1,1,jm+1,lm+1,im,js,je,nsoil)
1091  do l = 1, nsoil
1092  do j = jsta_2l, jend_2u
1093  do i = 1, im
1094  sh2o( i, j, l ) = dum3d( i, j, nsoil-l+1)
1095  end do
1096  end do
1097  end do
1098 
1099  varname='SEAICE'
1100  call getvariable(filename,datestr,datahandle,varname,dummy, &
1101  im,1,jm,1,im,js,je,1)
1102 
1103  do j = jsta_2l, jend_2u
1104  do i = 1, im
1105  sice( i, j ) = dummy( i, j )
1106  end do
1107  end do
1108 
1109  varname='TSLB'
1110  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1111  im+1,1,jm+1,lm+1,im,js,je,nsoil)
1112  do l = 1, nsoil
1113  do j = jsta_2l, jend_2u
1114  do i = 1, im
1115 ! stc ( i, j, l ) = dum3d ( i, j, l )
1116  stc( i, j, l ) = dum3d( i, j, nsoil-l+1)
1117  end do
1118  end do
1119  end do
1120 
1121 ! bitmask out high, middle, and low cloud cover
1122  do j = jsta_2l, jend_2u
1123  do i = 1, im
1124  cfrach( i, j ) = spval/100.
1125  cfracl( i, j ) = spval/100.
1126  cfracm( i, j ) = spval/100.
1127  end do
1128  end do
1129 
1130  do l = 1, lm
1131  do j = jsta_2l, jend_2u
1132  do i = 1, im
1133  cfr( i, j, l ) = spval
1134  end do
1135  end do
1136  end do
1137 
1138  varname='SR'
1139  call getvariable(filename,datestr,datahandle,varname,dummy, &
1140  im,1,jm,1,im,js,je,1)
1141  do j = jsta_2l, jend_2u
1142  do i = 1, im
1143  sr( i, j ) = dummy( i, j )
1144  end do
1145  end do
1146 
1147 ! WRF EM outputs 3D cloud cover now
1148  IF(modelname == 'RAPR')THEN
1149  varname='CLDFRA_BL'
1150  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1151  im+1,1,jm+1,lm+1,im,js,je,lm)
1152  do l=1,lm
1153  do j = jsta_2l, jend_2u
1154  do i = 1, im
1155  cfr( i, j, l ) = dum3d( i, j, l )
1156  end do
1157  end do
1158  end do
1159  ELSE
1160  varname='CLDFRA'
1161  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1162  im+1,1,jm+1,lm+1,im,js,je,lm)
1163  do l=1,lm
1164  do j = jsta_2l, jend_2u
1165  do i = 1, im
1166  cfr( i, j, l ) = dum3d( i, j, l )
1167  end do
1168  end do
1169  end do
1170  END IF
1171 
1172  varname='QC_BL'
1173  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1174  im+1,1,jm+1,lm+1,im,js,je,lm)
1175 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1176 ! & IM,1,JM,1,IM,JS,JE,1)
1177  do l=1,lm
1178  do j = jsta_2l, jend_2u
1179  do i = 1, im
1180  qc_bl( i, j, l ) = dum3d( i, j, l )
1181  end do
1182  end do
1183  end do
1184 
1185  call ext_ncd_get_dom_ti_real(datahandle,'DX',tmp, &
1186  1,ioutcount,istatus)
1187  dxval=nint(tmp)
1188  write(6,*) 'dxval= ', dxval
1189 #ifdef COMMCODE
1190  IF(modelname == 'NCAR' .OR. modelname == 'RAPR')THEN
1191  if(imp_physics/=5 .and. imp_physics/=0)then
1192 #else
1193  IF(modelname == 'RAPR')THEN
1194 ! J. Kenyon - 4 Apr 2019: revised cloud-cover diagnostics for RAP/HRRR
1195 #endif
1196  ptop_low = 64200.
1197  ptop_mid = 35000.
1198  ptop_high = 15000.
1199 
1200  do j = jsta_2l, jend_2u
1201  do i = 1, im
1202  cfracl(i,j)=0.
1203  cfracm(i,j)=0.
1204  cfrach(i,j)=0.
1205 
1206  do k = 1,lm
1207  if (pmid(i,j,k) >= ptop_low) then
1208  cfracl(i,j)=max(cfracl(i,j),cfr(i,j,k))
1209  elseif (pmid(i,j,k) < ptop_low .and. pmid(i,j,k) >= ptop_mid) then
1210  cfracm(i,j)=max(cfracm(i,j),cfr(i,j,k))
1211  elseif (pmid(i,j,k) < ptop_mid .and. pmid(i,j,k) >= ptop_high) then
1212  cfrach(i,j)=max(cfrach(i,j),cfr(i,j,k))
1213  endif
1214  enddo
1215 
1216  enddo
1217  enddo
1218 #ifdef COMMCODE
1219  ENDIF ! Not Ferrier or null mp_physics
1220 #endif
1221  ENDIF ! NCAR or RAPR
1222 
1223 ! GRAUP
1224 ! E. James - 8 Dec 2017: SMOKE from WRF-CHEM
1225 !
1226  varname='smoke'
1227  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1228  im+1,1,jm+1,lm+1,im,js,je,lm)
1229  do l = 1,lm
1230  do j = jsta_2l, jend_2u
1231  do i = 1, im
1232  smoke( i, j, l, 1) = dum3d( i, j, l )
1233  end do
1234  end do
1235  end do
1236 !!
1237 !! E. James - 8 Dec 2017: tracer_2a is anthropogenic aerosol, not used
1238 !!
1239 ! VarName='tracer_2a'
1240 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D, &
1241 ! IM+1,1,JM+1,LM+1,IM,JS,JE,LM)
1242 ! do l=1,lm
1243 ! do j = jsta_2l, jend_2u
1244 ! do i = 1, im
1245 ! SMOKE ( i, j, l, 2) = dum3d ( i, j, l )
1246 ! end do
1247 ! end do
1248 ! end do
1249 
1250 ! CRA DUST FROM WRF-CHEM
1251 if(1==2) then
1252  varname='DUST_1'
1253  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1254  im+1,1,jm+1,lm+1,im,js,je,lm)
1255 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1256 ! & IM,1,JM,1,IM,JS,JE,1)
1257  do l=1,lm
1258  do j = jsta_2l, jend_2u
1259  do i = 1, im
1260 ! CLDFRA( i, j ) = dummy ( i, j )
1261  dust( i, j, l, 1) = dum3d( i, j, l )
1262  end do
1263  end do
1264  end do
1265 
1266  varname='DUST_2'
1267  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1268  im+1,1,jm+1,lm+1,im,js,je,lm)
1269 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1270 ! & IM,1,JM,1,IM,JS,JE,1)
1271  do l=1,lm
1272  do j = jsta_2l, jend_2u
1273  do i = 1, im
1274 ! CLDFRA( i, j ) = dummy ( i, j )
1275  dust( i, j, l, 2) = dum3d( i, j, l )
1276  end do
1277  end do
1278  end do
1279 
1280  varname='DUST_3'
1281  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1282  im+1,1,jm+1,lm+1,im,js,je,lm)
1283 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1284 ! & IM,1,JM,1,IM,JS,JE,1)
1285  do l=1,lm
1286  do j = jsta_2l, jend_2u
1287  do i = 1, im
1288 ! CLDFRA( i, j ) = dummy ( i, j )
1289  dust( i, j, l, 3) = dum3d( i, j, l )
1290  end do
1291  end do
1292  end do
1293 
1294  varname='DUST_4'
1295  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1296  im+1,1,jm+1,lm+1,im,js,je,lm)
1297 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1298 ! & IM,1,JM,1,IM,JS,JE,1)
1299  do l=1,lm
1300  do j = jsta_2l, jend_2u
1301  do i = 1, im
1302 ! CLDFRA( i, j ) = dummy ( i, j )
1303  dust( i, j, l, 4) = dum3d( i, j, l )
1304  end do
1305  end do
1306  end do
1307 
1308  varname='DUST_5'
1309  call getvariable(filename,datestr,datahandle,varname,dum3d, &
1310  im+1,1,jm+1,lm+1,im,js,je,lm)
1311 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1312 ! & IM,1,JM,1,IM,JS,JE,1)
1313  do l=1,lm
1314  do j = jsta_2l, jend_2u
1315  do i = 1, im
1316 ! CLDFRA( i, j ) = dummy ( i, j )
1317  dust( i, j, l, 5) = dum3d( i, j, l )
1318  end do
1319  end do
1320  end do
1321 ! CRA
1322 ! For HRRR-CHEM
1323  varname='PM2_5_DRY'
1324  call getvariable(filename,datestr,datahandle,varname,dum3d,&
1325  im+1,1,jm+1,lm+1,im,js,je,lm)
1326 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1327 ! & IM,1,JM,1,IM,JS,JE,1)
1328  do l=1,lm
1329  do j = jsta_2l, jend_2u
1330  do i = 1, im
1331  dust( i, j, l, 6) = dum3d( i, j, l )
1332  end do
1333  end do
1334  end do
1335  varname='PM10'
1336  call getvariable(filename,datestr,datahandle,varname,dum3d,&
1337  im+1,1,jm+1,lm+1,im,js,je,lm)
1338 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1339 ! & IM,1,JM,1,IM,JS,JE,1)
1340  do l=1,lm
1341  do j = jsta_2l, jend_2u
1342  do i = 1, im
1343  dust( i, j, l, 7) = dum3d( i, j, l )
1344  end do
1345  end do
1346  end do
1347 
1348  varname='so2'
1349  call getvariable(filename,datestr,datahandle,varname,dum3d,&
1350  im+1,1,jm+1,lm+1,im,js,je,lm)
1351 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUM3D,
1352 ! & IM,1,JM,1,IM,JS,JE,1)
1353  do l=1,lm
1354  do j = jsta_2l, jend_2u
1355  do i = 1, im
1356  dust( i, j, l, 8) = dum3d( i, j, l )
1357  end do
1358  end do
1359  end do
1360 endif ! 1==2
1361 ! END HRRR-CHEM
1362 
1363 ! CRA
1364 
1365 ! Soil layer/depth - extract thickness of soil layers from wrf output
1366  sldpth=0.0 !Assign bogus value
1367 
1368  ! RUC LSM - use depths of center of soil layer
1369  IF(isf_surface_physics==3)then ! RUC LSM
1370  call getvariable(filename,datestr,datahandle,'ZS',sllevel, &
1371  nsoil,1,1,1,nsoil,1,1,1)
1372  print*,'SLLEVEL= ',(sllevel(n),n=1,nsoil)
1373  ELSE
1374  call getvariable(filename,datestr,datahandle,'DZS',sldpth, &
1375  nsoil,1,1,1,nsoil,1,1,1)
1376  print*,'SLDPTH= ',(sldpth(n),n=1,nsoil)
1377  END IF
1378 
1379 ! SRD
1380 ! get 2-d variables
1381 
1382  varname='WSPD10MAX'
1383  call getvariable(filename,datestr,datahandle,varname,dummy, &
1384  im,1,jm,1,im,js,je,1)
1385  do j = jsta_2l, jend_2u
1386  do i = 1, im
1387  wspd10max( i, j ) = dummy( i, j )
1388  end do
1389  end do
1390 ! print*,'WSPD10MAX at ',ii,jj,' = ',WSPD10MAX(ii,jj)
1391 
1392  varname='WSPD10UMAX'
1393  call getvariable(filename,datestr,datahandle,varname,dummy, &
1394  im,1,jm,1,im,js,je,1)
1395  do j = jsta_2l, jend_2u
1396  do i = 1, im
1397  wspd10umax( i, j ) = dummy( i, j )
1398  end do
1399  end do
1400 ! print*,'WSPD10UMAX at ',ii,jj,' = ',WSPD10UMAX(ii,jj)
1401 
1402  varname='WSPD10VMAX'
1403  call getvariable(filename,datestr,datahandle,varname,dummy, &
1404  im,1,jm,1,im,js,je,1)
1405  do j = jsta_2l, jend_2u
1406  do i = 1, im
1407  wspd10vmax( i, j ) = dummy( i, j )
1408  end do
1409  end do
1410 ! print*,'WSPD10VMAX at ',ii,jj,' = ',WSPD10VMAX(ii,jj)
1411 
1412  varname='W_UP_MAX'
1413  call getvariable(filename,datestr,datahandle,varname,dummy, &
1414  im,1,jm,1,im,js,je,1)
1415  do j = jsta_2l, jend_2u
1416  do i = 1, im
1417  w_up_max( i, j ) = dummy( i, j )
1418 ! print *,' w_up_max, i,j, = ', w_up_max(i,j)
1419  end do
1420  end do
1421 ! print*,'W_UP_MAX at ',ii,jj,' = ',W_UP_MAX(ii,jj)
1422 
1423  varname='W_DN_MAX'
1424  call getvariable(filename,datestr,datahandle,varname,dummy, &
1425  im,1,jm,1,im,js,je,1)
1426  do j = jsta_2l, jend_2u
1427  do i = 1, im
1428  w_dn_max( i, j ) = dummy( i, j )
1429  end do
1430  end do
1431 ! print*,'W_DN_MAX at ',ii,jj,' = ',W_DN_MAX(ii,jj)
1432 
1433  varname='W_MEAN'
1434  call getvariable(filename,datestr,datahandle,varname,dummy, &
1435  im,1,jm,1,im,js,je,1)
1436  do j = jsta_2l, jend_2u
1437  do i = 1, im
1438  w_mean( i, j ) = dummy( i, j )
1439  end do
1440  end do
1441 ! print*,'W_MEAN at ',ii,jj,' = ',W_MEAN(ii,jj)
1442 
1443  varname='REFD_MAX'
1444  call getvariable(filename,datestr,datahandle,varname,dummy, &
1445  im,1,jm,1,im,js,je,1)
1446  do j = jsta_2l, jend_2u
1447  do i = 1, im
1448  refd_max( i, j ) = dummy( i, j )
1449  end do
1450  end do
1451 ! print*,'REFD_MAX at ',ii,jj,' = ',REFD_MAX(ii,jj)
1452 
1453  varname='REFDM10C_MAX'
1454  call getvariable(filename,datestr,datahandle,varname,dummy, &
1455  im,1,jm,1,im,js,je,1)
1456  do j = jsta_2l, jend_2u
1457  do i = 1, im
1458  refdm10c_max( i, j ) = dummy( i, j )
1459  end do
1460  end do
1461 ! print*,'REFDM10C_MAX at ',ii,jj,' = ',REFDM10C_MAX(ii,jj)
1462 
1463 
1464  varname='UP_HELI_MAX'
1465  call getvariable(filename,datestr,datahandle,varname,dummy, &
1466  im,1,jm,1,im,js,je,1)
1467  do j = jsta_2l, jend_2u
1468  do i = 1, im
1469  up_heli_max( i, j ) = dummy( i, j )
1470  end do
1471  end do
1472 ! print*,'UP_HELI_MAX at ',ii,jj,' = ',UP_HELI_MAX(ii,jj)
1473 
1474  varname='UP_HELI_MAX16'
1475  call getvariable(filename,datestr,datahandle,varname,dummy, &
1476  im,1,jm,1,im,js,je,1)
1477  do j = jsta_2l, jend_2u
1478  do i = 1, im
1479  up_heli_max16( i, j ) = dummy( i, j )
1480  end do
1481  end do
1482 ! print*,'UP_HELI_MAX16 at ',ii,jj,' = ',UP_HELI_MAX16(ii,jj)
1483 
1484  varname='UP_HELI_MIN'
1485  call getvariable(filename,datestr,datahandle,varname,dummy, &
1486  im,1,jm,1,im,js,je,1)
1487  do j = jsta_2l, jend_2u
1488  do i = 1, im
1489  up_heli_min( i, j ) = dummy( i, j )
1490  end do
1491  end do
1492 ! print*,'UP_HELI_MIN at ',ii,jj,' = ',UP_HELI_MIN(ii,jj)
1493 
1494  varname='UP_HELI_MIN16'
1495  call getvariable(filename,datestr,datahandle,varname,dummy, &
1496  im,1,jm,1,im,js,je,1)
1497  do j = jsta_2l, jend_2u
1498  do i = 1, im
1499  up_heli_min16( i, j ) = dummy( i, j )
1500  end do
1501  end do
1502 ! print*,'UP_HELI_MIN16 at ',ii,jj,' = ',UP_HELI_MIN16(ii,jj)
1503 
1504  varname='UP_HELI_MAX02'
1505  call getvariable(filename,datestr,datahandle,varname,dummy, &
1506  im,1,jm,1,im,js,je,1)
1507  do j = jsta_2l, jend_2u
1508  do i = 1, im
1509  up_heli_max02( i, j ) = dummy( i, j )
1510  end do
1511  end do
1512 ! print*,'UP_HELI_MAX02 at ',ii,jj,' = ',UP_HELI_MAX02(ii,jj)
1513 
1514  varname='UP_HELI_MIN02'
1515  call getvariable(filename,datestr,datahandle,varname,dummy, &
1516  im,1,jm,1,im,js,je,1)
1517  do j = jsta_2l, jend_2u
1518  do i = 1, im
1519  up_heli_min02( i, j ) = dummy( i, j )
1520  end do
1521  end do
1522 ! print*,'UP_HELI_MIN02 at ',ii,jj,' = ',UP_HELI_MIN02(ii,jj)
1523 
1524  varname='UP_HELI_MAX03'
1525  call getvariable(filename,datestr,datahandle,varname,dummy, &
1526  im,1,jm,1,im,js,je,1)
1527  do j = jsta_2l, jend_2u
1528  do i = 1, im
1529  up_heli_max03( i, j ) = dummy( i, j )
1530  end do
1531  end do
1532 ! print*,'UP_HELI_MAX03 at ',ii,jj,' = ',UP_HELI_MAX03(ii,jj)
1533 
1534  varname='UP_HELI_MIN03'
1535  call getvariable(filename,datestr,datahandle,varname,dummy, &
1536  im,1,jm,1,im,js,je,1)
1537  do j = jsta_2l, jend_2u
1538  do i = 1, im
1539  up_heli_min03( i, j ) = dummy( i, j )
1540  end do
1541  end do
1542 ! print*,'UP_HELI_MIN03 at ',ii,jj,' = ',UP_HELI_MIN03(ii,jj)
1543 
1544  varname='REL_VORT_MAX'
1545  call getvariable(filename,datestr,datahandle,varname,dummy, &
1546  im,1,jm,1,im,js,je,1)
1547  do j = jsta_2l, jend_2u
1548  do i = 1, im
1549  rel_vort_max( i, j ) = dummy( i, j )
1550  end do
1551  end do
1552 ! print*,'REL_VORT_MAX at ',ii,jj,' = ',REL_VORT_MAX(ii,jj)
1553 
1554  varname='REL_VORT_MAX01'
1555  call getvariable(filename,datestr,datahandle,varname,dummy, &
1556  im,1,jm,1,im,js,je,1)
1557  do j = jsta_2l, jend_2u
1558  do i = 1, im
1559  rel_vort_max01( i, j ) = dummy( i, j )
1560  end do
1561  end do
1562 ! print*,'REL_VORT_MAX01 at ',ii,jj,' = ',REL_VORT_MAX01(ii,jj)
1563 
1564  varname='GRPL_MAX'
1565  call getvariable(filename,datestr,datahandle,varname,dummy, &
1566  im,1,jm,1,im,js,je,1)
1567  do j = jsta_2l, jend_2u
1568  do i = 1, im
1569  grpl_max( i, j ) = dummy( i, j )
1570  end do
1571  end do
1572 ! print*,'GRPL_MAX at ',ii,jj,' = ',GRPL_MAX(ii,jj)
1573 
1574 
1575  varname='HAIL_MAXK1'
1576  call getvariable(filename,datestr,datahandle,varname,dummy, &
1577  im,1,jm,1,im,js,je,1)
1578  do j = jsta_2l, jend_2u
1579  do i = 1, im
1580  hail_maxk1( i, j ) = dummy( i, j )
1581  end do
1582  end do
1583 ! print*,'HAIL_MAXK1 at ',ii,jj,' = ',HAIL_MAXK1(ii,jj)
1584 
1585  varname='HAIL_MAX2D'
1586  call getvariable(filename,datestr,datahandle,varname,dummy, &
1587  im,1,jm,1,im,js,je,1)
1588  do j = jsta_2l, jend_2u
1589  do i = 1, im
1590  hail_max2d( i, j ) = dummy( i, j )
1591  end do
1592  end do
1593 ! print*,'HAIL_MAX2D at ',ii,jj,' = ',HAIL_MAX2D(ii,jj)
1594 
1595  varname='HAILCAST_DIAM_MAX'
1596  call getvariable(filename,datestr,datahandle,varname,dummy, &
1597  im,1,jm,1,im,js,je,1)
1598  do j = jsta_2l, jend_2u
1599  do i = 1, im
1600  hail_maxhailcast( i, j ) = dummy( i, j )
1601  end do
1602  end do
1603 ! print*,'HAIL_MAXHAILCAST at ',ii,jj,' = ',HAIL_MAXHAILCAST(ii,jj)
1604 
1605  varname='UH'
1606  call getvariable(filename,datestr,datahandle,varname,dummy, &
1607  im,1,jm,1,im,js,je,1)
1608  do j = jsta_2l, jend_2u
1609  do i = 1, im
1610  up_heli( i, j ) = dummy( i, j )
1611  end do
1612  end do
1613 ! print*,'UP_HELI at ',ii,jj,' = ',UP_HELI(ii,jj)
1614 
1615  varname='UH16'
1616  call getvariable(filename,datestr,datahandle,varname,dummy, &
1617  im,1,jm,1,im,js,je,1)
1618  do j = jsta_2l, jend_2u
1619  do i = 1, im
1620  up_heli16( i, j ) = dummy( i, j )
1621  end do
1622  end do
1623 ! print*,'UP_HELI16 at ',ii,jj,' = ',UP_HELI16(ii,jj)
1624 
1625  varname='LTG1_MAX'
1626  call getvariable(filename,datestr,datahandle,varname,dummy, &
1627  im,1,jm,1,im,js,je,1)
1628  do j = jsta_2l, jend_2u
1629  do i = 1, im
1630  ltg1_max( i, j ) = dummy( i, j )
1631  end do
1632  end do
1633 ! print*,'LTG1_MAX at ',ii,jj,' = ',LTG1_MAX(ii,jj)
1634 
1635  varname='LTG2_MAX'
1636  call getvariable(filename,datestr,datahandle,varname,dummy, &
1637  im,1,jm,1,im,js,je,1)
1638  do j = jsta_2l, jend_2u
1639  do i = 1, im
1640  ltg2_max( i, j ) = dummy( i, j )
1641  end do
1642  end do
1643 ! print*,'LTG2_MAX at ',ii,jj,' = ',LTG2_MAX(ii,jj)
1644 
1645  varname='LTG3_MAX'
1646  call getvariable(filename,datestr,datahandle,varname,dummy, &
1647  im,1,jm,1,im,js,je,1)
1648  do j = jsta_2l, jend_2u
1649  do i = 1, im
1650  ltg3_max( i, j ) = dummy( i, j )
1651  end do
1652  end do
1653 ! print*,'LTG3_MAX at ',ii,jj,' = ',LTG3_MAX(ii,jj)
1654 
1655  varname='NCI_LTG'
1656  call getvariable(filename,datestr,datahandle,varname,dummy, &
1657  im,1,jm,1,im,js,je,1)
1658  do j = jsta_2l, jend_2u
1659  do i = 1, im
1660  nci_ltg( i, j ) = dummy( i, j )
1661  end do
1662  end do
1663 ! print*,'NCI_LTG at ',ii,jj,' = ',NCI_LTG(ii,jj)
1664 
1665  varname='NCA_LTG'
1666  call getvariable(filename,datestr,datahandle,varname,dummy, &
1667  im,1,jm,1,im,js,je,1)
1668  do j = jsta_2l, jend_2u
1669  do i = 1, im
1670  nca_ltg( i, j ) = dummy( i, j )
1671  end do
1672  end do
1673 ! print*,'NCA_LTG at ',ii,jj,' = ',NCA_LTG(ii,jj)
1674 
1675  varname='NCI_WQ'
1676  call getvariable(filename,datestr,datahandle,varname,dummy, &
1677  im,1,jm,1,im,js,je,1)
1678  do j = jsta_2l, jend_2u
1679  do i = 1, im
1680  nci_wq( i, j ) = dummy( i, j )
1681  end do
1682  end do
1683 ! print*,'NCI_WQ at ',ii,jj,' = ',NCI_WQ(ii,jj)
1684 
1685  varname='NCA_WQ'
1686  call getvariable(filename,datestr,datahandle,varname,dummy, &
1687  im,1,jm,1,im,js,je,1)
1688  do j = jsta_2l, jend_2u
1689  do i = 1, im
1690  nca_wq( i, j ) = dummy( i, j )
1691  end do
1692  end do
1693 ! print*,'NCA_WQ at ',ii,jj,' = ',NCA_WQ(ii,jj)
1694 
1695  varname='NCI_REFD'
1696  call getvariable(filename,datestr,datahandle,varname,dummy, &
1697  im,1,jm,1,im,js,je,1)
1698  do j = jsta_2l, jend_2u
1699  do i = 1, im
1700  nci_refd( i, j ) = dummy( i, j )
1701  end do
1702  end do
1703 ! print*,'NCI_REFD at ',ii,jj,' = ',NCI_REFD(ii,jj)
1704 
1705  varname='NCA_REFD'
1706  call getvariable(filename,datestr,datahandle,varname,dummy, &
1707  im,1,jm,1,im,js,je,1)
1708  do j = jsta_2l, jend_2u
1709  do i = 1, im
1710  nca_refd( i, j ) = dummy( i, j )
1711  end do
1712  end do
1713 ! print*,'NCA_REFD at ',ii,jj,' = ',NCA_REFD(ii,jj)
1714 !
1715 ! SRD
1716 !
1717 ! CRA REFLECTIVITY VARIABLES
1718  varname='REFL_10CM'
1719  call getvariable(filename,datestr,datahandle,varname,dum3d,&
1720  im+1,1,jm+1,lm+1,im,js,je,lm)
1721  do l=1,lm
1722  do j = jsta_2l, jend_2u
1723  do i = 1, im
1724  ref_10cm( i, j, l) = dum3d( i, j, l )
1725  end do
1726  end do
1727  end do
1728  deallocate(dum3d)
1729 
1730  varname='COMPOSITE_REFL_10CM'
1731  call getvariable(filename,datestr,datahandle,varname,dummy, &
1732  im,1,jm,1,im,js,je,1)
1733  do j = jsta_2l, jend_2u
1734  do i = 1, im
1735  refc_10cm( i, j ) = dummy( i, j )
1736  end do
1737  end do
1738 
1739  varname='REFL_10CM_1KM'
1740  call getvariable(filename,datestr,datahandle,varname,dummy, &
1741  im,1,jm,1,im,js,je,1)
1742  do j = jsta_2l, jend_2u
1743  do i = 1, im
1744  ref1km_10cm( i, j ) = dummy( i, j )
1745  end do
1746  end do
1747 
1748  varname='REFL_10CM_4KM'
1749  call getvariable(filename,datestr,datahandle,varname,dummy, &
1750  im,1,jm,1,im,js,je,1)
1751  do j = jsta_2l, jend_2u
1752  do i = 1, im
1753  ref4km_10cm( i, j ) = dummy( i, j )
1754  end do
1755  end do
1756 ! CRA
1757 ! get 2-d variables
1758 
1759  varname='U10'
1760  call getvariable(filename,datestr,datahandle,varname,dummy, &
1761  im,1,jm,1,im,js,je,1)
1762  do j = jsta_2l, jend_2u
1763  do i = 1, im
1764  IF(submodelname == 'RTMA' .and. modelname == 'RAPR')THEN !use 1st level of unstaggered U for U10
1765  u10( i, j ) = uh( i, j, lm )
1766  ELSE
1767  u10( i, j ) = dummy( i, j )
1768  ENDIF
1769  end do
1770  end do
1771  varname='V10'
1772  call getvariable(filename,datestr,datahandle,varname,dummy, &
1773  im,1,jm,1,im,js,je,1)
1774  do j = jsta_2l, jend_2u
1775  do i = 1, im
1776  IF( submodelname == 'RTMA' .and. modelname == 'RAPR')then!use 1st level of unstaggered V for V10
1777  v10( i, j ) = vh( i, j, lm )
1778  ELSE
1779  v10( i, j ) = dummy( i, j )
1780  ENDIF
1781  end do
1782  end do
1783 ! print*,'V10 at ',ii,jj,' = ',V10(ii,jj)
1784 
1785 ! RAP/HRRR time-averaged wind
1786  varname='U10MEAN'
1787  call getvariable(filename,datestr,datahandle,varname,dummy, &
1788  im,1,jm,1,im,js,je,1)
1789  do j = jsta_2l, jend_2u
1790  do i = 1, im
1791  u10mean( i, j ) = dummy( i, j )
1792  end do
1793  end do
1794 ! print*,'U10mean at ',ii,jj,' = ',U10mean(ii,jj)
1795 !
1796  varname='V10MEAN'
1797  call getvariable(filename,datestr,datahandle,varname,dummy, &
1798  im,1,jm,1,im,js,je,1)
1799  do j = jsta_2l, jend_2u
1800  do i = 1, im
1801  v10mean( i, j ) = dummy( i, j )
1802  end do
1803  end do
1804 ! print*,'V10mean at ',ii,jj,' = ',V10mean(ii,jj)
1805 !
1806  varname='SPDUV10MEAN'
1807  call getvariable(filename,datestr,datahandle,varname,dummy, &
1808  im,1,jm,1,im,js,je,1)
1809  do j = jsta_2l, jend_2u
1810  do i = 1, im
1811  spduv10mean( i, j ) = dummy( i, j )
1812  end do
1813  end do
1814 ! print*,'SPDUV10mean at ',ii,jj,' = ',SPDUV10mean(ii,jj)
1815 !
1816 
1817  do j = jsta_2l, jend_2u
1818  do i = 1, im
1819  th10( i, j ) = spval
1820  q10( i, j ) = spval
1821  end do
1822  end do
1823 
1824 ! get 2-m theta
1825  varname='TH2'
1826  call getvariable(filename,datestr,datahandle,varname,dummy, &
1827  im,1,jm,1,im,js,je,1)
1828  do j = jsta_2l, jend_2u
1829  do i = 1, im
1830  tshltr( i, j ) = dummy( i, j )
1831  end do
1832  end do
1833 ! print*,'TSHLTR at ',ii,jj,' = ',TSHLTR(ii,jj)
1834 ! get 2-m mixing ratio
1835  varname='Q2'
1836  call getvariable(filename,datestr,datahandle,varname,dummy, &
1837  im,1,jm,1,im,js,je,1)
1838  do j = jsta_2l, jend_2u
1839  do i = 1, im
1840  mrshltr( i, j ) = dummy(i, j ) ! Shelter Mixing ratio
1841  IF(modelname == 'RAPR')THEN
1842 ! QV2M = first level QV
1843 ! QV2M ( i, j ) = q ( i, j, lm )/(1.-q ( i, j, lm )) ! 1st level mix. ratio
1844 ! QSHLTR ( i, j ) = q ( i, j, lm ) ! 1st level spec. humidity
1845 ! QV2M = diagnosed in WRF 2-m QV
1846  qv2m( i, j ) = dummy( i, j ) ! 2-m mix. ratio
1847  qshltr( i, j ) = dummy( i, j )/(1.0+dummy( i, j )) ! 2-m spec. hum.
1848  qvl1( i, j ) = q( i, j, lm ) ! spec. humidity at lev. 1
1849  ELSE
1850 !HC QSHLTR ( i, j ) = dummy ( i, j )
1851 !HC CONVERT FROM MIXING RATIO TO SPECIFIC HUMIDITY
1852  qv2m( i, j ) = dummy( i, j )
1853  qshltr( i, j ) = dummy( i, j )/(1.0+dummy( i, j ))
1854  ENDIF
1855  end do
1856  end do
1857 ! print*,'QSHLTR at ',ii,jj,' = ',QSHLTR(ii,jj)
1858 
1859  IF(modelname == 'RAPR')THEN
1860  varname='MAVAIL'
1861  ELSE
1862  varname='SMSTAV'
1863  END IF
1864 
1865  call getvariable(filename,datestr,datahandle,varname,dummy, &
1866  im,1,jm,1,im,js,je,1)
1867  do j = jsta_2l, jend_2u
1868  do i = 1, im
1869  smstav( i, j ) = dummy( i, j )
1870  end do
1871  end do
1872 
1873 ! VarName='SMSTOT'
1874 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
1875 ! IM,1,JM,1,IM,JS,JE,1)
1876 ! do j = jsta_2l, jend_2u
1877 ! do i = 1, im
1878 ! SMSTOT ( i, j ) = dummy ( i, j )
1879 ! end do
1880 ! end do
1881 
1882 !mhu VarName='SSROFF'
1883  varname='SFROFF'
1884  call getvariable(filename,datestr,datahandle,varname,dummy, &
1885  im,1,jm,1,im,js,je,1)
1886  do j = jsta_2l, jend_2u
1887  do i = 1, im
1888  ssroff( i, j ) = dummy( i, j )
1889  end do
1890  end do
1891  varname='UDROFF'
1892  call getvariable(filename,datestr,datahandle,varname,dummy, &
1893  im,1,jm,1,im,js,je,1)
1894  do j = jsta_2l, jend_2u
1895  do i = 1, im
1896  bgroff( i, j ) = dummy( i, j )
1897  end do
1898  end do
1899 
1900 ! VarName='SFCEVP'
1901 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
1902 ! IM,1,JM,1,IM,JS,JE,1)
1903 ! do j = jsta_2l, jend_2u
1904 ! do i = 1, im
1905 ! SFCEVP( i, j ) = dummy ( i, j )
1906 ! end do
1907 ! end do
1908 ! print*,'SFCEVP at ',ii,jj,' = ',SFCEVP(ii,jj)
1909 
1910 ! VarName='SFCEXC'
1911 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
1912 ! IM,1,JM,1,IM,JS,JE,1)
1913 ! do j = jsta_2l, jend_2u
1914 ! do i = 1, im
1915 ! SFCEXC ( i, j ) = dummy ( i, j )
1916 ! end do
1917 ! end do
1918 
1919  varname='VEGFRA'
1920  call getvariable(filename,datestr,datahandle,varname,dummy, &
1921  im,1,jm,1,im,js,je,1)
1922  do j = jsta_2l, jend_2u
1923  do i = 1, im
1924  vegfrc( i, j ) = dummy( i, j )/100.
1925  end do
1926  end do
1927  print*,'VEGFRC at ',ii,jj,' = ',vegfrc(ii,jj)
1928 
1929  varname='SHDMIN'
1930  call getvariable(filename,datestr,datahandle,varname,dummy, &
1931  im,1,jm,1,im,js,je,1)
1932  do j = jsta_2l, jend_2u
1933  do i = 1, im
1934  shdmin( i, j ) = dummy( i, j )/100.
1935  end do
1936  end do
1937  print*,'SHDMIN at ',ii,jj,' = ',shdmin(ii,jj)
1938 
1939  varname='SHDMAX'
1940  call getvariable(filename,datestr,datahandle,varname,dummy, &
1941  im,1,jm,1,im,js,je,1)
1942  do j = jsta_2l, jend_2u
1943  do i = 1, im
1944  shdmax( i, j ) = dummy( i, j )/100.
1945  end do
1946  end do
1947  print*,'SHDMAX at ',ii,jj,' = ',shdmax(ii,jj)
1948 
1949  varname='LAI'
1950  call getvariable(filename,datestr,datahandle,varname,dummy, &
1951  im,1,jm,1,im,js,je,1)
1952  do j = jsta_2l, jend_2u
1953  do i = 1, im
1954  lai( i, j ) = dummy( i, j )
1955  end do
1956  end do
1957  print*,'LAI at ',ii,jj,' = ',lai(ii,jj)
1958 
1959  varname='ACSNOW'
1960  call getvariable(filename,datestr,datahandle,varname,dummy, &
1961  im,1,jm,1,im,js,je,1)
1962  do j = jsta_2l, jend_2u
1963  do i = 1, im
1964  acsnow( i, j ) = dummy( i, j )
1965  end do
1966  end do
1967  print*,'maxval ACSNOW: ', maxval(acsnow)
1968  varname='ACSNOM'
1969  call getvariable(filename,datestr,datahandle,varname,dummy, &
1970  im,1,jm,1,im,js,je,1)
1971  do j = jsta_2l, jend_2u
1972  do i = 1, im
1973  acsnom( i, j ) = dummy( i, j )
1974  end do
1975  end do
1976  varname='CANWAT'
1977  call getvariable(filename,datestr,datahandle,varname,dummy, &
1978  im,1,jm,1,im,js,je,1)
1979  do j = jsta_2l, jend_2u
1980  do i = 1, im
1981  cmc( i, j ) = dummy( i, j )
1982  end do
1983  end do
1984  varname='SST'
1985  call getvariable(filename,datestr,datahandle,varname,dummy, &
1986  im,1,jm,1,im,js,je,1)
1987  do j = jsta_2l, jend_2u
1988  do i = 1, im
1989  sst( i, j ) = dummy( i, j )
1990  end do
1991  end do
1992 ! print*,'SST at ',ii,jj,' = ',sst(ii,jj)
1993  varname='THZ0'
1994  call getvariable(filename,datestr,datahandle,varname,dummy, &
1995  im,1,jm,1,im,js,je,1)
1996  do j = jsta_2l, jend_2u
1997  do i = 1, im
1998  thz0( i, j ) = dummy( i, j )
1999  end do
2000  end do
2001 ! print*,'THZ0 at ',ii,jj,' = ',THZ0(ii,jj)
2002 ! VarName='QZ0'
2003 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2004 ! IM,1,JM,1,IM,JS,JE,1)
2005 ! do j = jsta_2l, jend_2u
2006 ! do i = 1, im
2007 ! QZ0 ( i, j ) = dummy ( i, j )
2008 ! end do
2009 ! end do
2010 ! print*,'QZ0 at ',ii,jj,' = ',QZ0(ii,jj)
2011 ! VarName='UZ0'
2012 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2013 ! IM,1,JM,1,IM,JS,JE,1)
2014 ! do j = jsta_2l, jend_2u
2015 ! do i = 1, im
2016 ! UZ0 ( i, j ) = dummy ( i, j )
2017 ! end do
2018 ! end do
2019 ! print*,'UZ0 at ',ii,jj,' = ',UZ0(ii,jj)
2020 ! VarName='VZ0'
2021 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2022 ! IM,1,JM,1,IM,JS,JE,1)
2023 ! do j = jsta_2l, jend_2u
2024 ! do i = 1, im
2025 ! VZ0 ( i, j ) = dummy ( i, j )
2026 ! end do
2027 ! end do
2028 ! print*,'VZ0 at ',ii,jj,' = ',VZ0(ii,jj)
2029 ! VarName='QSFC'
2030 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2031 ! IM,1,JM,1,IM,JS,JE,1)
2032 ! do j = jsta_2l, jend_2u
2033 ! do i = 1, im
2034 ! QS ( i, j ) = dummy ( i, j )
2035 ! QVG ( i, j ) = dummy ( i, j )/(1.-dummy ( i, j ))
2036 ! end do
2037 ! end do
2038 ! print*,'QS at ',ii,jj,' = ',QS(ii,jj)
2039 
2040  IF(modelname == 'RAPR')THEN
2041  varname='ZNT'
2042  call getvariable(filename,datestr,datahandle,varname,dummy, &
2043  im,1,jm,1,im,js,je,1)
2044  do j = jsta_2l, jend_2u
2045  do i = 1, im
2046  z0( i, j ) = dummy( i, j )
2047  end do
2048  end do
2049  ELSE
2050  varname='Z0'
2051  call getvariable(filename,datestr,datahandle,varname,dummy, &
2052  im,1,jm,1,im,js,je,1)
2053  do j = jsta_2l, jend_2u
2054  do i = 1, im
2055  z0( i, j ) = dummy( i, j )
2056  end do
2057  end do
2058  END IF
2059 
2060 ! VarName='USTAR'
2061  varname='UST'
2062  call getvariable(filename,datestr,datahandle,varname,dummy, &
2063  im,1,jm,1,im,js,je,1)
2064  do j = jsta_2l, jend_2u
2065  do i = 1, im
2066  ustar( i, j ) = dummy( i, j )
2067  end do
2068  end do
2069 ! print*,'USTAR at ',ii,jj,' = ',USTAR(ii,jj)
2070 
2071 ! IF(MODELNAME == 'RAPR')THEN
2072 ! VarName='FLHC'
2073 ! ELSE
2074 ! VarName='AKHS'
2075 ! ENDIF
2076 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2077 ! IM,1,JM,1,IM,JS,JE,1)
2078 ! do j = jsta_2l, jend_2u
2079 ! do i = 1, im
2080 ! AKHS ( i, j ) = dummy ( i, j )
2081 ! end do
2082 ! end do
2083 ! VarName='AKMS'
2084 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY, &
2085 ! IM,1,JM,1,IM,JS,JE,1)
2086 ! do j = jsta_2l, jend_2u
2087 ! do i = 1, im
2088 ! AKMS ( i, j ) = dummy ( i, j )
2089 ! end do
2090 ! end do
2091 
2092 !
2093 ! In my version, variable is TSK (skin temp, not skin pot temp)
2094 !
2095 !mp call getVariable(fileName,DateStr,DataHandle,'THSK',DUMMY,
2096  varname='TSK'
2097  call getvariable(filename,datestr,datahandle,varname,dummy, &
2098  im,1,jm,1,im,js,je,1)
2099  do j = jsta_2l, jend_2u
2100  do i = 1, im
2101 !HC THS ( i, j ) = dummy ( i, j ) ! this is WRONG (should be theta)
2102 !HC CONVERT SKIN TEMPERATURE TO SKIN POTENTIAL TEMPERATURE
2103 ! CHC: deriving outgoing longwave fluxes by assuming emmissitivity=1
2104  radot( i, j ) = dummy(i,j)**4.0/stbol
2105  ths( i, j ) = dummy( i, j ) &
2106  *(p1000/pint(i,j,nint(lmh(i,j))+1))**capa
2107  end do
2108  end do
2109 ! print*,'THS at ',ii,jj,' = ',THS(ii,jj)
2110 
2111  varname='EMISS'
2112  IF(modelname == 'RAPR')THEN
2113 ! Update "RADOT" variable (calculated above) using model emissivity
2114  call getvariable(filename,datestr,datahandle,varname,dummy, &
2115  im,1,jm,1,im,js,je,1)
2116  do j = jsta_2l, jend_2u
2117  do i = 1, im
2118  radot( i, j ) = radot(i, j) * dummy( i, j )
2119  end do
2120  end do
2121  END IF
2122 
2123 !C
2124 !CMP
2125 !C
2126 !C RAINC is "ACCUMULATED TOTAL CUMULUS PRECIPITATION"
2127 !C RAINNC is "ACCUMULATED TOTAL GRID SCALE PRECIPITATION"
2128 
2129  write(6,*) 'getting RAINC'
2130  varname='RAINC'
2131  call getvariable(filename,datestr,datahandle,varname,dummy, &
2132  im,1,jm,1,im,js,je,1)
2133  do j = jsta_2l, jend_2u
2134  do i = 1, im
2135  cuprec( i, j ) = dummy( i, j ) * 0.001
2136  end do
2137  end do
2138 ! print*,'CUPREC at ',ii,jj,' = ',CUPREC(ii,jj)
2139  write(6,*) 'getting RAINNC'
2140  varname='RAINNC'
2141  call getvariable(filename,datestr,datahandle,varname,dummy, &
2142  im,1,jm,1,im,js,je,1)
2143  do j = jsta_2l, jend_2u
2144  do i = 1, im
2145  ancprc( i, j ) = dummy( i, j )* 0.001
2146  end do
2147  end do
2148 ! print*,'ANCPRC at ',ii,jj,' = ',ANCPRC(ii,jj)
2149  write(6,*) 'past getting RAINNC'
2150 
2151  do j = jsta_2l, jend_2u
2152  do i = 1, im
2153  acprec(i,j)=ancprc(i,j)+cuprec(i,j)
2154  end do
2155  end do
2156 
2157 !-- RAINC_bucket is "ACCUMULATED CUMULUS PRECIPITATION OVER BUCKET_DT PERIODS OF TIME"
2158 
2159  write(6,*) 'getting PREC_ACC_C, [mm] '
2160 ! VarName='RAINC_BUCKET'
2161  varname='PREC_ACC_C'
2162  call getvariable(filename,datestr,datahandle,varname,dummy, &
2163  im,1,jm,1,im,js,je,1)
2164  do j = jsta_2l, jend_2u
2165  do i = 1, im
2166  rainc_bucket( i, j ) = dummy( i, j )
2167  end do
2168  end do
2169 
2170 !-- RAINC_bucket1 is "ACCUMULATED CUMULUS PRECIPITATION OVER BUCKET_DT1 PERIODS OF TIME"
2171 
2172  write(6,*) 'getting PREC_ACC_C1, [mm] '
2173  varname='PREC_ACC_C1'
2174  call getvariable(filename,datestr,datahandle,varname,dummy, &
2175  im,1,jm,1,im,js,je,1)
2176  do j = jsta_2l, jend_2u
2177  do i = 1, im
2178  rainc_bucket1( i, j ) = dummy( i, j )
2179  end do
2180  end do
2181 
2182 !-- RAINNC_bucket is "ACCUMULATED GRID SCALE PRECIPITATION OVER BUCKET_DT PERIODS OF TIME"
2183 
2184  write(6,*) 'getting PREC_ACC_NC, [mm]'
2185 ! VarName='RAINNC_BUCKET'
2186  varname='PREC_ACC_NC'
2187  call getvariable(filename,datestr,datahandle,varname,dummy, &
2188  im,1,jm,1,im,js,je,1)
2189  do j = jsta_2l, jend_2u
2190  do i = 1, im
2191  rainnc_bucket( i, j ) = dummy( i, j )
2192  end do
2193  end do
2194 
2195 !-- RAINNC_bucket1 is "ACCUMULATED GRID SCALE PRECIPITATION OVER BUCKET_DT1 PERIODS OF TIME"
2196 
2197  write(6,*) 'getting PREC_ACC_NC1, [mm]'
2198  varname='PREC_ACC_NC1'
2199  call getvariable(filename,datestr,datahandle,varname,dummy, &
2200  im,1,jm,1,im,js,je,1)
2201  do j = jsta_2l, jend_2u
2202  do i = 1, im
2203  rainnc_bucket1( i, j ) = dummy( i, j )
2204  end do
2205  end do
2206 
2207  do j = jsta_2l, jend_2u
2208  do i = 1, im
2209  pcp_bucket(i,j)=rainc_bucket(i,j)+rainnc_bucket(i,j)
2210  pcp_bucket1(i,j)=rainc_bucket1(i,j)+rainnc_bucket1(i,j)
2211  end do
2212  end do
2213 
2214  varname='RAINCV'
2215  dummy=0.0
2216  call getvariable(filename,datestr,datahandle,varname,dummy, &
2217  im,1,jm,1,im,js,je,1)
2218  do j = jsta_2l, jend_2u
2219  do i = 1, im
2220 !-- CPRATE is in [m] per time step
2221  cprate( i, j ) = dummy( i, j )* 0.001
2222  end do
2223  end do
2224 
2225 
2226  varname='RAINNCV'
2227  dummy2=0.0
2228  call getvariable(filename,datestr,datahandle,varname,dummy2, &
2229  im,1,jm,1,im,js,je,1)
2230  do j = jsta_2l, jend_2u
2231  do i = 1, im
2232 !-- PREC is in [m] per time step
2233  prec( i, j ) = (dummy( i, j )+dummy2(i,j))* 0.001
2234  end do
2235  end do
2236 
2237  varname='SNOWNCV'
2238  call getvariable(filename,datestr,datahandle,varname,dummy, &
2239  im,1,jm,1,im,js,je,1)
2240  do j = jsta_2l, jend_2u
2241  do i = 1, im
2242 !-- SNOW is in [m] per time sep
2243  snownc( i, j ) = dummy( i, j ) * 0.001
2244  end do
2245  end do
2246 
2247 !-- SNOW_bucket is "ACCUMULATED GRID SCALE SNOW OVER BUCKET_DT PERIODS OF TIME"
2248 
2249  write(6,*) 'getting SNOW_ACC_NC, [mm] '
2250 ! VarName='SNOW_BUCKET'
2251  varname='SNOW_ACC_NC'
2252  call getvariable(filename,datestr,datahandle,varname,dummy, &
2253  im,1,jm,1,im,js,je,1)
2254  do j = jsta_2l, jend_2u
2255  do i = 1, im
2256  snow_bucket( i, j ) = dummy( i, j )
2257  end do
2258  end do
2259 
2260 !-- SNOW_bucket1 is "ACCUMULATED GRID SCALE SNOW OVER BUCKET_DT1 PERIODS OF TIME"
2261 
2262  write(6,*) 'getting SNOW_ACC_NC1, [mm] '
2263  varname='SNOW_ACC_NC1'
2264  call getvariable(filename,datestr,datahandle,varname,dummy, &
2265  im,1,jm,1,im,js,je,1)
2266  do j = jsta_2l, jend_2u
2267  do i = 1, im
2268  snow_bucket1( i, j ) = dummy( i, j )
2269  end do
2270  end do
2271 
2272 !-- GRAUP_bucket is "ACCUMULATED GRID SCALE GRAUPEL OVER BUCKET_DT PERIODS OF TIME"
2273 
2274  write(6,*) 'getting GRAUP_ACC_NC, [mm] '
2275  varname='GRAUP_ACC_NC'
2276  call getvariable(filename,datestr,datahandle,varname,dummy, &
2277  im,1,jm,1,im,js,je,1)
2278  do j = jsta_2l, jend_2u
2279  do i = 1, im
2280  graup_bucket( i, j ) = dummy( i, j )
2281  end do
2282  end do
2283 
2284 !-- GRAUP_bucket1 is "ACCUMULATED GRID SCALE GRAUPEL OVER BUCKET_DT1 PERIODS OF TIME"
2285 
2286  write(6,*) 'getting GRAUP_ACC_NC1, [mm] '
2287  varname='GRAUP_ACC_NC1'
2288  call getvariable(filename,datestr,datahandle,varname,dummy, &
2289  im,1,jm,1,im,js,je,1)
2290  do j = jsta_2l, jend_2u
2291  do i = 1, im
2292  graup_bucket1( i, j ) = dummy( i, j )
2293  end do
2294  end do
2295 
2296  varname='ACGRAUP'
2297  call getvariable(filename,datestr,datahandle,varname,dummy, &
2298  im,1,jm,1,im,js,je,1)
2299  do j = jsta_2l, jend_2u
2300  do i = 1, im
2301  acgraup( i, j ) = dummy( i, j )
2302  end do
2303  end do
2304 
2305  varname='ACFRAIN'
2306  call getvariable(filename,datestr,datahandle,varname,dummy, &
2307  im,1,jm,1,im,js,je,1)
2308  do j = jsta_2l, jend_2u
2309  do i = 1, im
2310  acfrain( i, j ) = dummy( i, j )
2311  end do
2312  end do
2313 
2314  varname='GRAUPELNCV'
2315  call getvariable(filename,datestr,datahandle,varname,dummy, &
2316  im,1,jm,1,im,js,je,1)
2317  do j = jsta_2l, jend_2u
2318  do i = 1, im
2319 !-- GRAUPEL in in [m] per time step
2320  graupelnc( i, j ) = dummy( i, j ) * 0.001
2321  end do
2322  end do
2323 
2324 
2325  varname='ALBSOL'
2326  call getvariable(filename,datestr,datahandle,varname,dummy, &
2327  im,1,jm,1,im,js,je,1)
2328  do j = jsta_2l, jend_2u
2329  do i = 1, im
2330  albedo( i, j ) = dummy( i, j )
2331  end do
2332  end do
2333 !
2334 ! GSD trunk as 2014.4.14
2335 ! VarName='GSW'
2336 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY,&
2337 ! IM,1,JM,1,IM,JS,JE,1)
2338 ! do j = jsta_2l, jend_2u
2339 ! do i = 1, im
2340 ! RSWIN ( i, j ) = dummy ( i, j )
2341 ! HCHUANG: GSW is actually net downward shortwave in ncar wrf
2342 ! RSWIN ( i, j ) = dummy ( i, j )/(1.0-albedo(i,j))
2343 ! RSWOUT ( i, j ) = RSWIN ( i, j ) - dummy ( i, j )
2344 ! end do
2345 ! end do
2346 
2347  varname='SWDOWN'
2348  call getvariable(filename,datestr,datahandle,varname,dummy, &
2349  im,1,jm,1,im,js,je,1)
2350  do j = jsta_2l, jend_2u
2351  do i = 1, im
2352 ! HCHUANG: SWDOWN is actually net downward shortwave in ncar wrf
2353  rswin( i, j ) = dummy( i, j )
2354  rswout( i, j ) = rswin( i, j ) * albedo( i, j )
2355  end do
2356  end do
2357 
2358  varname='SWDDNI'
2359 ! Shortwave surface downward direct normal irradiance
2360  call getvariable(filename,datestr,datahandle,varname,dummy, &
2361  im,1,jm,1,im,js,je,1)
2362  do j = jsta_2l, jend_2u
2363  do i = 1, im
2364  swddni( i, j ) = dummy( i, j )
2365  end do
2366  end do
2367 
2368  varname='SWDDIF'
2369 ! Shortwave surface downward diffuse horizontal irradiance
2370  call getvariable(filename,datestr,datahandle,varname,dummy, &
2371  im,1,jm,1,im,js,je,1)
2372  do j = jsta_2l, jend_2u
2373  do i = 1, im
2374  swddif( i, j ) = dummy( i, j )
2375  end do
2376  end do
2377 
2378  varname='SWDNBC'
2379 ! Shortwave surface downward clear-sky shortwave irradiance
2380  call getvariable(filename,datestr,datahandle,varname,dummy, &
2381  im,1,jm,1,im,js,je,1)
2382  do j = jsta_2l, jend_2u
2383  do i = 1, im
2384  swdnbc( i, j ) = dummy( i, j )
2385  end do
2386  end do
2387 
2388  varname='SWDDNIC'
2389 ! Clear-sky shortwave surface downward direct normal irradiance
2390  call getvariable(filename,datestr,datahandle,varname,dummy, &
2391  im,1,jm,1,im,js,je,1)
2392  do j = jsta_2l, jend_2u
2393  do i = 1, im
2394  swddnic( i, j ) = dummy( i, j )
2395  end do
2396  end do
2397 
2398  varname='SWDDIFC'
2399 ! Clear-sky shortwave surface downward diffuse horizontal irradiance
2400  call getvariable(filename,datestr,datahandle,varname,dummy, &
2401  im,1,jm,1,im,js,je,1)
2402  do j = jsta_2l, jend_2u
2403  do i = 1, im
2404  swddifc( i, j ) = dummy( i, j )
2405  end do
2406  end do
2407 
2408  varname='SWUPBC'
2409 ! Clear-sky surface upwelling shortwave flux
2410  call getvariable(filename,datestr,datahandle,varname,dummy, &
2411  im,1,jm,1,im,js,je,1)
2412  do j = jsta_2l, jend_2u
2413  do i = 1, im
2414  swupbc( i, j ) = dummy( i, j )
2415  end do
2416  end do
2417 
2418  varname='SWUPT'
2419 ! Upward shortwave flux at top of atmosphere
2420  call getvariable(filename,datestr,datahandle,varname,dummy, &
2421  im,1,jm,1,im,js,je,1)
2422  do j = jsta_2l, jend_2u
2423  do i = 1, im
2424  swupt( i, j ) = dummy( i, j )
2425  end do
2426  end do
2427 
2428 !
2429 ! E. James - 8 Dec 2017: Fire Radiative Power from HRRR-smoke
2430 !
2431  varname='MEAN_FRP'
2432 ! Mean fire radiative power
2433  call getvariable(filename,datestr,datahandle,varname,dummy, &
2434  im,1,jm,1,im,js,je,1)
2435  do j = jsta_2l, jend_2u
2436  do i = 1, im
2437  mean_frp( i, j ) = dummy( i, j )
2438  end do
2439  end do
2440 
2441  varname='TAOD5502D'
2442 ! Total aerosol optical depth
2443  call getvariable(filename,datestr,datahandle,varname,dummy, &
2444  im,1,jm,1,im,js,je,1)
2445  do j = jsta_2l, jend_2u
2446  do i = 1, im
2447  taod5502d( i, j ) = dummy( i, j )
2448  end do
2449  end do
2450 
2451  varname='AERASY2D'
2452 ! Aerosol asymmetry parameter
2453  call getvariable(filename,datestr,datahandle,varname,dummy, &
2454  im,1,jm,1,im,js,je,1)
2455  do j = jsta_2l, jend_2u
2456  do i = 1, im
2457  aerasy2d( i, j ) = dummy( i, j )
2458  end do
2459  end do
2460 
2461  varname='AERSSA2D'
2462 ! Aerosol single-scattering albedo
2463  call getvariable(filename,datestr,datahandle,varname,dummy, &
2464  im,1,jm,1,im,js,je,1)
2465  do j = jsta_2l, jend_2u
2466  do i = 1, im
2467  aerssa2d( i, j ) = dummy( i, j )
2468  end do
2469  end do
2470 
2471  varname='LWP'
2472 ! Liquid water path
2473  call getvariable(filename,datestr,datahandle,varname,dummy, &
2474  im,1,jm,1,im,js,je,1)
2475  do j = jsta_2l, jend_2u
2476  do i = 1, im
2477  lwp( i, j ) = dummy( i, j )
2478  end do
2479  end do
2480 
2481  varname='IWP'
2482 ! Ice water path
2483  call getvariable(filename,datestr,datahandle,varname,dummy, &
2484  im,1,jm,1,im,js,je,1)
2485  do j = jsta_2l, jend_2u
2486  do i = 1, im
2487  iwp( i, j ) = dummy( i, j )
2488  end do
2489  end do
2490 
2491 ! time_averaged SWDOWN
2492  varname='SWRADMEAN'
2493  call getvariable(filename,datestr,datahandle,varname,dummy, &
2494  im,1,jm,1,im,js,je,1)
2495  do j = jsta_2l, jend_2u
2496  do i = 1, im
2497 ! averaged incoming solar radiation
2498  swradmean( i, j ) = dummy( i, j )
2499  end do
2500  end do
2501  print*,'SWRADmean at ',ii,jj,' = ',swradmean(ii,jj)
2502 
2503 ! time_averaged SWNORM
2504  varname='SWNORMMEAN'
2505  call getvariable(filename,datestr,datahandle,varname,dummy, &
2506  im,1,jm,1,im,js,je,1)
2507  do j = jsta_2l, jend_2u
2508  do i = 1, im
2509 ! averaged incoming solar radiation
2510  swnormmean( i, j ) = dummy( i, j )
2511  end do
2512  end do
2513  print*,'SWNORMmean at ',ii,jj,' = ',swnormmean(ii,jj)
2514 
2515  varname='GLW'
2516 ! Downwelling longwave flux at surface
2517  call getvariable(filename,datestr,datahandle,varname,dummy, &
2518  im,1,jm,1,im,js,je,1)
2519  do j = jsta_2l, jend_2u
2520  do i = 1, im
2521  rlwin( i, j ) = dummy( i, j )
2522  end do
2523  end do
2524 ! ncar wrf does not output sigt4 so make sig4=sigma*tlmh**4
2525  do j = jsta_2l, jend_2u
2526  do i = 1, im
2527  tlmh=t(i,j,nint(lmh(i,j)))
2528  sigt4( i, j ) = 5.67e-8*tlmh*tlmh*tlmh*tlmh
2529  end do
2530  end do
2531 
2532  varname='LWDNBC'
2533 ! Clear-sky downwelling longwave flux at surface
2534  call getvariable(filename,datestr,datahandle,varname,dummy, &
2535  im,1,jm,1,im,js,je,1)
2536  do j = jsta_2l, jend_2u
2537  do i = 1, im
2538  lwdnbc( i, j ) = dummy( i, j )
2539  end do
2540  end do
2541 
2542  varname='LWUPBC'
2543 ! Clear-sky upwelling longwave flux at surface
2544  call getvariable(filename,datestr,datahandle,varname,dummy, &
2545  im,1,jm,1,im,js,je,1)
2546  do j = jsta_2l, jend_2u
2547  do i = 1, im
2548  lwupbc( i, j ) = dummy( i, j )
2549  end do
2550  end do
2551 
2552 ! Top of the atmosphere outgoing LW radiation
2553  varname='OLR'
2554  call getvariable(filename,datestr,datahandle,varname,dummy, &
2555  im,1,jm,1,im,js,je,1)
2556  do j = jsta_2l, jend_2u
2557  do i = 1, im
2558  rlwtoa( i, j ) = dummy( i, j )
2559  end do
2560  end do
2561 
2562 
2563 ! NCAR WRF does not output accumulated fluxes so set the bitmap of these fluxes to 0
2564  do j = jsta_2l, jend_2u
2565  do i = 1, im
2566 ! RLWTOA(I,J)=SPVAL
2567  rswinc(i,j)=spval
2568  aswin(i,j)=spval
2569  aswout(i,j)=spval
2570  alwin(i,j)=spval
2571  alwout(i,j)=spval
2572  alwtoa(i,j)=spval
2573  aswtoa(i,j)=spval
2574  ardlw=1.0
2575  ardsw=1.0
2576  nrdlw=1
2577  nrdsw=1
2578  end do
2579  end do
2580 
2581  varname='TMN'
2582  call getvariable(filename,datestr,datahandle,varname,dummy, &
2583  im,1,jm,1,im,js,je,1)
2584  do j = jsta_2l, jend_2u
2585  do i = 1, im
2586  tg( i, j ) = dummy( i, j )
2587  soiltb( i, j ) = dummy( i, j )
2588  end do
2589  end do
2590 
2591  varname='HFX'
2592  call getvariable(filename,datestr,datahandle,varname,dummy, &
2593  im,1,jm,1,im,js,je,1)
2594  do j = jsta_2l, jend_2u
2595  do i = 1, im
2596  twbs(i,j)= dummy( i, j )
2597 ! SFCSHX ( i, j ) = dummy ( i, j )
2598 ! ASRFC=1.0
2599  end do
2600  end do
2601 
2602 ! latent heat flux
2603  IF(isf_surface_physics/=3) then
2604  varname='LH'
2605  call getvariable(filename,datestr,datahandle,varname,dummy, &
2606  im,1,jm,1,im,js,je,1)
2607  do j = jsta_2l, jend_2u
2608  do i = 1, im
2609  qwbs(i,j) = dummy( i, j )
2610 ! SFCLHX ( i, j ) = dummy ( i, j )
2611  end do
2612  end do
2613  else
2614  varname='QFX'
2615  call getvariable(filename,datestr,datahandle,varname,dummy, &
2616  im,1,jm,1,im,js,je,1)
2617  do j = jsta_2l, jend_2u
2618  do i = 1, im
2619  qwbs(i,j) = dummy( i, j ) * lheat
2620  end do
2621  end do
2622  ENDIF
2623 
2624 ! ground heat fluxes
2625  varname='GRDFLX'
2626  call getvariable(filename,datestr,datahandle,varname,dummy, &
2627  im,1,jm,1,im,js,je,1)
2628  do j = jsta_2l, jend_2u
2629  do i = 1, im
2630  grnflx(i,j) = dummy( i, j )
2631  end do
2632  end do
2633 
2634 ! NCAR WRF does not output accumulated fluxes so bitmask out these fields
2635  do j = jsta_2l, jend_2u
2636  do i = 1, im
2637  sfcshx(i,j)=spval
2638  sfclhx(i,j)=spval
2639  subshx(i,j)=spval
2640  snopcx(i,j)=spval
2641  sfcuvx(i,j)=spval
2642  potevp(i,j)=spval
2643  ncfrcv(i,j)=spval
2644  ncfrst(i,j)=spval
2645  asrfc=1.0
2646  nsrfc=1
2647  end do
2648  end do
2649 
2650 ! VarName='WEASD'
2651 ! Snow water equivalent
2652  varname='SNOW' ! WRF V2 replace WEASD with SNOW
2653  call getvariable(filename,datestr,datahandle,varname,dummy, &
2654  im,1,jm,1,im,js,je,1)
2655  do j = jsta_2l, jend_2u
2656  do i = 1, im
2657  if( dummy( i, j ) <= 5000.0 .and. dummy( i, j ) >=0.0) then
2658  sno( i, j ) = dummy( i, j )
2659  elseif( dummy( i, j ) > 5000.0) then
2660  sno( i, j ) = 5000.0
2661  write(*,*) 'too large SNOW=',i,j,dummy( i, j )
2662  elseif( dummy( i, j ) < 0.0 ) then
2663  sno( i, j ) = 0.0
2664  write(*,*) 'negative SNOW=',i,j,dummy( i, j )
2665  else
2666  sno( i, j ) = 0.0
2667  write(*,*) 'strange SNOW=',i,j,dummy( i, j )
2668  endif
2669  end do
2670  end do
2671 ! Snow depth
2672  varname='SNOWH'
2673  call getvariable(filename,datestr,datahandle,varname,dummy, &
2674  im,1,jm,1,im,js,je,1)
2675  do j = jsta_2l, jend_2u
2676  do i = 1, im
2677  if( dummy( i, j ) <= 50.0 .and. dummy( i, j ) >=0.0) then
2678  si( i, j ) = dummy( i, j ) * 1000.
2679  elseif( dummy( i, j ) > 50.0) then
2680  si( i, j ) = 50.0 * 1000.
2681  write(*,*) 'too large SNOWH=',i,j,dummy( i, j )
2682  elseif( dummy( i, j ) < 0.0 ) then
2683  si( i, j ) = 0.0
2684  write(*,*) 'negative SNOWH=',i,j,dummy( i, j )
2685  else
2686  si( i, j ) = 0.0
2687  write(*,*) 'strange SNOWH=',i,j,dummy( i, j )
2688  endif
2689  end do
2690  end do
2691 
2692 ! snow cover
2693  varname='SNOWC'
2694  call getvariable(filename,datestr,datahandle,varname,dummy, &
2695  im,1,jm,1,im,js,je,1)
2696  do j = jsta_2l, jend_2u
2697  do i = 1, im
2698  pctsno( i, j ) = dummy( i, j )
2699  end do
2700  end do
2701 
2702 ! Accumulated grid-scale snow and ice precipitation
2703  varname='SNOWNC'
2704  call getvariable(filename,datestr,datahandle,varname,dummy, &
2705  im,1,jm,1,im,js,je,1)
2706  do j = jsta_2l, jend_2u
2707  do i = 1, im
2708  snonc( i, j ) = dummy( i, j )
2709  end do
2710  end do
2711 
2712 ! snowfall density
2713  varname='RHOSNF'
2714  call getvariable(filename,datestr,datahandle,varname,dummy, &
2715  im,1,jm,1,im,js,je,1)
2716  do j = jsta_2l, jend_2u
2717  do i = 1, im
2718  snfden( i, j ) = max(0.,dummy( i, j ))
2719  end do
2720  end do
2721  print *,' MIN/MAX SNFDEN ',minval(snfden),maxval(snfden)
2722 
2723 ! snowfall accumulation
2724  varname='SNOWFALLAC'
2725  call getvariable(filename,datestr,datahandle,varname,dummy, &
2726  im,1,jm,1,im,js,je,1)
2727  do j = jsta_2l, jend_2u
2728  do i = 1, im
2729  sndepac( i, j ) = dummy( i, j )
2730  end do
2731  end do
2732  print *,' MIN/MAX SNDEPAC ',minval(sndepac),maxval(sndepac)
2733 
2734 ! snow temperature at the interface of 2 snow layers
2735  varname='SOILT1'
2736  call getvariable(filename,datestr,datahandle,varname,dummy, &
2737  im,1,jm,1,im,js,je,1)
2738  do j = jsta_2l, jend_2u
2739  do i = 1, im
2740  tsnow( i, j ) = dummy( i, j )
2741  end do
2742  end do
2743 
2744 ! GET VEGETATION TYPE
2745 
2746  call getivariablen(filename,datestr,datahandle,'IVGTYP',idummy, &
2747  im,1,jm,1,im,js,je,1)
2748 ! print*,'sample VEG TYPE',IDUMMY(20,20)
2749  do j = jsta_2l, jend_2u
2750  do i = 1, im
2751  ivgtyp( i, j ) = idummy( i, j )
2752  end do
2753  end do
2754 
2755  varname='ISLTYP'
2756  call getivariablen(filename,datestr,datahandle,varname,idummy, &
2757  im,1,jm,1,im,js,je,1)
2758  do j = jsta_2l, jend_2u
2759  do i = 1, im
2760  isltyp( i, j ) = idummy( i, j )
2761  end do
2762  end do
2763  print*,'MAX ISLTYP=', maxval(idummy)
2764 
2765 ! VarName='ISLOPE'
2766 ! call getIVariableN(fileName,DateStr,DataHandle,VarName,IDUMMY, &
2767 ! IM,1,JM,1,IM,JS,JE,1)
2768 ! do j = jsta_2l, jend_2u
2769 ! do i = 1, im
2770 ! ISLOPE( i, j ) = idummy ( i, j )
2771 ! end do
2772 ! end do
2773 
2774 
2775 ! XLAND 1 land 2 sea
2776  varname='XLAND'
2777  call getvariable(filename,datestr,datahandle,varname,dummy, &
2778  im,1,jm,1,im,js,je,1)
2779  do j = jsta_2l, jend_2u
2780  do i = 1, im
2781  sm( i, j ) = dummy( i, j ) - 1.0
2782  end do
2783  end do
2784 
2785 ! PBL depth
2786  varname='PBLH'
2787  call getvariable(filename,datestr,datahandle,varname,dummy, &
2788  im,1,jm,1,im,js,je,1)
2789  do j = jsta_2l, jend_2u
2790  do i = 1, im
2791  pblh( i, j ) = dummy( i, j )
2792  end do
2793  end do
2794  IF(modelname == 'RAPR')THEN
2795 ! PBL depth from GSD
2796  delta_theta4gust=0.5
2797  do j = jsta_2l, jend_2u
2798  do i = 1, im
2799 !! Is there any mixed layer at all?
2800  if (thv(i,j,lm-1) < (thv(i,j,lm) + delta_theta4gust)) then
2801  zsf=zint(i,j,nint(lmh(i,j))+1)
2802 !! Calculate k1 level as first above PBL top
2803  do 34 k=3,lm
2804  k1 = k
2805 !! - give theta-v at the sfc a 0.5K boost in
2806 !! the PBL height definition
2807  if (thv(i,j,lm-k+1)>(thv(i,j,lm) + delta_theta4gust)) &
2808 ! go to 341
2809  exit
2810  34 continue
2811  continue
2812  !341 continue
2813  zpbltop = zmid(i,j,lm-k1+1) + &
2814  ((thv(i,j,lm)+delta_theta4gust)-thv(i,j,lm-k1+1)) &
2815  * (zmid(i,j,lm-k1+2)-zmid(i,j,lm-k1+1)) &
2816  / (thv(i,j,lm-k1+2) - thv(i,j,lm-k1+1))
2817  pblhgust( i, j ) = zpbltop - zsf
2818  else
2819  pblhgust( i, j ) = 0.
2820  endif
2821  end do
2822  end do
2823  ENDIF
2824 
2825  varname='XLAT'
2826  call getvariable(filename,datestr,datahandle,varname,dummy, &
2827  im,1,jm,1,im,js,je,1)
2828  do j = jsta_2l, jend_2u
2829  do i = 1, im
2830  gdlat( i, j ) = dummy( i, j )
2831 ! compute F = 2*omg*sin(xlat)
2832  f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr)
2833  end do
2834  end do
2835 ! pos north
2836 ! print*,'GDLAT at ',ii,jj,' = ',GDLAT(ii,jj)
2837  print*,'read past GDLAT'
2838  varname='XLONG'
2839  call getvariable(filename,datestr,datahandle,varname,dummy, &
2840  im,1,jm,1,im,js,je,1)
2841  do j = jsta_2l, jend_2u
2842  do i = 1, im
2843  gdlon( i, j ) = dummy( i, j )
2844 ! if(abs(GDLAT(i,j)-20.0)<0.5 .and. abs(GDLON(I,J)
2845 ! 1 +157.0)<5.)print*
2846 ! 2 ,'Debug:I,J,GDLON,GDLAT,SM,HGT,psfc= ',i,j,GDLON(i,j)
2847 ! 3 ,GDLAT(i,j),SM(i,j),FIS(i,j)/G,PINT(I,j,lm+1)
2848  end do
2849  end do
2850 ! print*,'GDLON at ',ii,jj,' = ',GDLON(ii,jj)
2851  print*,'read past GDLON'
2852 ! pos east
2853  call collect_loc(gdlat,dummy)
2854  if(me==0)then
2855  latstart=nint(dummy(1,1)*gdsdegr)
2856  latlast=nint(dummy(im,jm)*gdsdegr)
2857 ! print*,'LL corner from model output= ',dummy(1,1)
2858 ! print*,'LR corner from model output= ',dummy(im,1)
2859 ! print*,'UL corner from model output= ',dummy(1,jm)
2860 ! print*,'UR corner from model output= ',dummy(im,jm)
2861  end if
2862  write(6,*) 'laststart,latlast B calling bcast= ',latstart,latlast
2863  call mpi_bcast(latstart,1,mpi_integer,0,mpi_comm_comp,irtn)
2864  call mpi_bcast(latlast,1,mpi_integer,0,mpi_comm_comp,irtn)
2865  write(6,*) 'laststart,latlast A calling bcast= ',latstart,latlast
2866  call collect_loc(gdlon,dummy)
2867  if(me==0)then
2868  if(dummy(1,1)<0.0) dummy(1,1)=360.0+dummy(1,1)
2869  if(dummy(im,jm)<0.0) dummy(im,jm)=360.0+dummy(im,jm)
2870  lonstart=nint(dummy(1,1)*gdsdegr)
2871  lonlast=nint(dummy(im,jm)*gdsdegr)
2872 ! print*,'LL corner from model output= ',dummy(1,1)
2873 ! print*,'LR corner from model output= ',dummy(im,1)
2874 ! print*,'UL corner from model output= ',dummy(1,jm)
2875 ! print*,'UR corner from model output= ',dummy(im,jm)
2876  end if
2877  write(6,*)'lonstart,lonlast B calling bcast=',lonstart,lonlast
2878  call mpi_bcast(lonstart,1,mpi_integer,0,mpi_comm_comp,irtn)
2879  call mpi_bcast(lonlast,1,mpi_integer,0,mpi_comm_comp,irtn)
2880  write(6,*)'lonstart,lonlast A calling bcast= ',lonstart,lonlast
2881 !
2882 ! obtain map scale factor
2883 ! VarName='msft'
2884  varname='MAPFAC_M'
2885 
2886  allocate(msft(im,jsta_2l:jend_2u))
2887 
2888  call getvariable(filename,datestr,datahandle,varname,dummy, &
2889  im,1,jm,1,im,js,je,1)
2890  do j = jsta_2l, jend_2u
2891  do i = 1, im
2892  msft( i, j ) = dummy( i, j )
2893  end do
2894  end do
2895 
2896 ! physics calling frequency
2897  varname='STEPBL'
2898  call getivariablen(filename,datestr,datahandle,varname,nphs, &
2899  1,1,1,1,1,1,1,1)
2900 
2901 
2902 ! ncdump -h
2903 
2904 ! ncar wrf does not output zenith angle so make czen=czmean so that
2905 ! RSWIN can be output normally in SURFCE
2906  IF(modelname /= 'RAPR')THEN
2907  do j = jsta_2l, jend_2u
2908  do i = 1, im
2909  czen( i, j ) = 1.0
2910  czmean( i, j ) = czen( i, j )
2911  end do
2912  end do
2913  ELSE
2914 
2915  jdn=iw3jdn(idat(3),idat(1),idat(2))
2916  do j=jsta,jend
2917  do i=1,im
2918  call zensun(jdn,float(idat(4)),gdlat(i,j),gdlon(i,j) &
2919  ,pi,sun_zenith,sun_azimuth)
2920  temp=sun_zenith/rtd
2921  czen(i,j)=cos(temp)
2922  czmean( i, j ) = czen( i, j )
2923  end do
2924  end do
2925  print*,'sample RAPR zenith angle=',acos(czen(ii,jj))*rtd
2926  ENDIF
2927 
2928 
2929 !!
2930 !!
2931 !!
2932  write(6,*) 'filename in INITPOST=', filename,' is'
2933 
2934 ! status=nf_open(filename,NF_NOWRITE,ncid)
2935 ! write(6,*) 'returned ncid= ', ncid
2936 ! status=nf_get_att_real(ncid,varid,'DX',tmp)
2937 ! dxval=int(tmp)
2938 ! status=nf_get_att_real(ncid,varid,'DY',tmp)
2939 ! dyval=int(tmp)
2940 ! status=nf_get_att_real(ncid,varid,'CEN_LAT',tmp)
2941 ! cenlat=int(1000.*tmp)
2942 ! status=nf_get_att_real(ncid,varid,'CEN_LON',tmp)
2943 ! cenlon=int(1000.*tmp)
2944 ! status=nf_get_att_real(ncid,varid,'TRUELAT1',tmp)
2945 ! truelat1=int(1000.*tmp)
2946 ! status=nf_get_att_real(ncid,varid,'TRUELAT2',tmp)
2947 ! truelat2=int(1000.*tmp)
2948 ! status=nf_get_att_real(ncid,varid,'MAP_PROJ',tmp)
2949 ! maptype=int(tmp)
2950 ! status=nf_close(ncid)
2951 
2952 ! dxval=30000.
2953 ! dyval=30000.
2954 !
2955 ! write(6,*) 'dxval= ', dxval
2956 ! write(6,*) 'dyval= ', dyval
2957 ! write(6,*) 'cenlat= ', cenlat
2958 ! write(6,*) 'cenlon= ', cenlon
2959 ! write(6,*) 'truelat1= ', truelat1
2960 ! write(6,*) 'truelat2= ', truelat2
2961 ! write(6,*) 'maptype is ', maptype
2962 !
2963 !tgs call ext_ncd_get_dom_ti_real(DataHandle,'DX',tmp, &
2964 ! 1,ioutcount,istatus)
2965 ! dxval=nint(tmp)
2966 ! write(6,*) 'dxval= ', dxval
2967  call ext_ncd_get_dom_ti_real(datahandle,'DY',tmp, &
2968  1,ioutcount,istatus)
2969  dyval=nint(tmp)
2970  write(6,*) 'dyval= ', dyval
2971  call ext_ncd_get_dom_ti_real(datahandle,'CEN_LAT',tmp, &
2972  1,ioutcount,istatus)
2973  cenlat=nint(gdsdegr*tmp)
2974  write(6,*) 'cenlat= ', cenlat
2975  call ext_ncd_get_dom_ti_real(datahandle,'CEN_LON',tmp, &
2976  1,ioutcount,istatus)
2977  if(tmp < 0) tmp=360.0 + tmp
2978  cenlon=nint(gdsdegr*tmp)
2979  write(6,*) 'cenlon= ', cenlon
2980  call ext_ncd_get_dom_ti_integer(datahandle,'MAP_PROJ',itmp, &
2981  1,ioutcount,istatus)
2982  maptype=itmp
2983  write(6,*) 'maptype is ', maptype
2984  if(maptype/=6)then
2985  call ext_ncd_get_dom_ti_real(datahandle,'TRUELAT1',tmp, &
2986  1,ioutcount,istatus)
2987  truelat1=nint(gdsdegr*tmp)
2988  write(6,*) 'truelat1= ', truelat1
2989  if(maptype/=2)then !PS projection excluded
2990  call ext_ncd_get_dom_ti_real(datahandle,'TRUELAT2',tmp, &
2991  1,ioutcount,istatus)
2992  truelat2=nint(gdsdegr*tmp)
2993  write(6,*) 'truelat2= ', truelat2
2994  endif
2995  endif
2996  call ext_ncd_get_dom_ti_real(datahandle,'STAND_LON',tmp, &
2997  1,ioutcount,istatus)
2998  if(tmp < 0) tmp=360.0 + tmp
2999  standlon=nint(gdsdegr*tmp)
3000  write(6,*) 'STANDLON= ', standlon
3001 
3002 !MEB not sure how to get these
3003  do j = jsta_2l, jend_2u
3004  do i = 1, im
3005  dx( i, j ) = dxval/msft(i,j)
3006  dy( i, j ) = dyval/msft(i,j)
3007  end do
3008  end do
3009  ii=im/2
3010  jj=(jend+jsta)/2
3011  print*,'sample dx,dy,msft=',ii,jj,dx(ii,jj),dy(ii,jj) &
3012  ,msft(ii,jj)
3013 
3014 
3015 ! Convert DXVAL and DYVAL for ARW rotated latlon from meters to radian
3016  if(maptype==6)then
3017  dxval=(dxval * 360.)/(erad*2.*pi)*gdsdegr
3018  dyval=(dyval * 360.)/(erad*2.*pi)*gdsdegr
3019 
3020  print*,'dx and dy for arw rotated latlon= ', &
3021  dxval,dyval
3022  end if
3023 
3024 !tgs Define smoothing flag for isobaric output
3025  IF(modelname == 'RAPR')THEN
3026  smflag=.true.
3027  ELSE
3028  smflag=.false.
3029  ENDIF
3030 
3031 ! generate look up table for lifted parcel calculations
3032 
3033  thl=210.
3034  plq=70000.
3035 
3036  CALL table(ptbl,ttbl,pt, &
3037  rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
3038 
3039  CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
3040 
3041 !
3042 !
3043  IF(me==0)THEN
3044  WRITE(6,*)' SPL (POSTED PRESSURE LEVELS) BELOW: '
3045  WRITE(6,51) (spl(l),l=1,lsm)
3046  50 FORMAT(14(f4.1,1x))
3047  51 FORMAT(8(f8.1,1x))
3048  ENDIF
3049 !
3050 ! COMPUTE DERIVED TIME STEPPING CONSTANTS.
3051 !
3052 !need to get DT
3053  call ext_ncd_get_dom_ti_real(datahandle,'DT',tmp,1,ioutcount,istatus)
3054  dt=abs(tmp)
3055  print*,'DT= ',dt
3056 
3057 !need to get period of time for precipitation buckets
3058  call ext_ncd_get_dom_ti_real(datahandle,'PREC_ACC_DT',tmp,1,ioutcount,istatus)
3059  prec_acc_dt=abs(tmp)
3060  print*,'PREC_ACC_DT= ',prec_acc_dt
3061 
3062 !need to get period of time for precipitation bucket 1 (15-min precip)
3063 !talk to Tanya about getting this output in wrfout file
3064  prec_acc_dt1=15.0
3065  print*,'PREC_ACC_DT1= ',prec_acc_dt1
3066 
3067 ! DT = 120. !MEB need to get DT
3068  nphs = 1 !CHUANG SET IT TO 1 BECAUSE ALL THE INST PRECIP ARE ACCUMULATED 1 TIME STEP
3069  dtq2 = dt * nphs !MEB need to get physics DT
3070  tsph = 3600./dt !MEB need to get DT
3071 ! Randomly specify accumulation period because WRF EM does not
3072 ! output accumulation fluxes yet and accumulated fluxes are bit
3073 ! masked out
3074 
3075  tsrfc=1.0
3076  trdlw=1.0
3077  trdsw=1.0
3078  theat=1.0
3079  tclod=1.0
3080 
3081  tprec=float(nprec)/tsph
3082  IF(nprec==0)tprec=float(ifhr) !in case buket does not get emptied
3083  print*,'NPREC,TPREC = ',nprec,tprec
3084 
3085 !tgs TPREC=float(ifhr) ! WRF EM does not empty precip buket at all
3086 
3087 ! TSRFC=float(NSRFC)/TSPH
3088 ! TRDLW=float(NRDLW)/TSPH
3089 ! TRDSW=float(NRDSW)/TSPH
3090 ! THEAT=float(NHEAT)/TSPH
3091 ! TCLOD=float(NCLOD)/TSPH
3092 ! TPREC=float(NPREC)/TSPH
3093  print*,'TSRFC TRDLW TRDSW= ',tsrfc, trdlw, trdsw
3094 !MEB need to get DT
3095 
3096 !how am i going to get this information?
3097 ! NPREC = INT(TPREC *TSPH+D50)
3098 ! NHEAT = INT(THEAT *TSPH+D50)
3099 ! NCLOD = INT(TCLOD *TSPH+D50)
3100 ! NRDSW = INT(TRDSW *TSPH+D50)
3101 ! NRDLW = INT(TRDLW *TSPH+D50)
3102 ! NSRFC = INT(TSRFC *TSPH+D50)
3103 !how am i going to get this information?
3104 !
3105 ! IF(ME==0)THEN
3106 ! WRITE(6,*)' '
3107 ! WRITE(6,*)'DERIVED TIME STEPPING CONSTANTS'
3108 ! WRITE(6,*)' NPREC,NHEAT,NSRFC : ',NPREC,NHEAT,NSRFC
3109 ! WRITE(6,*)' NCLOD,NRDSW,NRDLW : ',NCLOD,NRDSW,NRDLW
3110 ! ENDIF
3111 !
3112 
3113 ! VarName='RAINCV'
3114 ! call getVariable(fileName,DateStr,DataHandle,VarName,DUMMY,
3115 ! & IM,1,JM,1,IM,JS,JE,1)
3116 ! do j = jsta_2l, jend_2u
3117 ! do i = 1, im
3118 ! CUPPT ( i, j ) = dummy ( i, j )* 0.001*(TRDLW*3600.)
3119 ! end do
3120 ! end do
3121 
3122 
3123 
3124 ! COMPUTE DERIVED MAP OUTPUT CONSTANTS.
3125  DO l = 1,lsm
3126  alsl(l) = alog(spl(l))
3127  END DO
3128 ! close up shop
3129  call ext_ncd_ioclose( datahandle, status )
3130 !
3131 !HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
3132  if(me==0)then
3133  print*,'writing out igds'
3134  igdout=110
3135 ! open(igdout,file='griddef.out',form='unformatted'
3136 ! + ,status='unknown')
3137  if(maptype == 1)THEN ! Lambert conformal
3138  WRITE(igdout)3
3139  WRITE(6,*)'igd(1)=',3
3140  WRITE(igdout)im
3141  WRITE(igdout)jm
3142  WRITE(igdout)latstart
3143  WRITE(igdout)lonstart
3144  WRITE(igdout)8
3145 ! WRITE(igdout)CENLON
3146  WRITE(igdout)standlon
3147  WRITE(igdout)dxval
3148  WRITE(igdout)dyval
3149  WRITE(igdout)0
3150  WRITE(igdout)64
3151  WRITE(igdout)truelat2
3152  WRITE(igdout)truelat1
3153  WRITE(igdout)255
3154  ELSE IF(maptype == 2)THEN !Polar stereographic
3155  WRITE(igdout)5
3156  WRITE(igdout)im
3157  WRITE(igdout)jm
3158  WRITE(igdout)latstart
3159  WRITE(igdout)lonstart
3160  WRITE(igdout)8
3161  WRITE(igdout)cenlon
3162  WRITE(igdout)dxval
3163  WRITE(igdout)dyval
3164  WRITE(igdout)0
3165  WRITE(igdout)64
3166  WRITE(igdout)truelat2 !Assume projection at +-90
3167  WRITE(igdout)truelat1
3168  WRITE(igdout)255
3169  ! Note: The calculation of the map scale factor at the standard
3170  ! lat/lon and the PSMAPF
3171  ! Get map factor at 60 degrees (N or S) for PS projection, which will
3172  ! be needed to correctly define the DX and DY values in the GRIB GDS
3173  if (truelat1 < 0.) THEN
3174  lat = -60.
3175  else
3176  lat = 60.
3177  end if
3178 
3179  CALL msfps(lat,truelat1*0.001,psmapf)
3180 
3181  ELSE IF(maptype == 3)THEN !Mercator
3182  WRITE(igdout)1
3183  WRITE(igdout)im
3184  WRITE(igdout)jm
3185  WRITE(igdout)latstart
3186  WRITE(igdout)lonstart
3187  WRITE(igdout)8
3188  WRITE(igdout)latlast
3189  WRITE(igdout)lonlast
3190  WRITE(igdout)truelat1
3191  WRITE(igdout)0
3192  WRITE(igdout)64
3193  WRITE(igdout)dxval
3194  WRITE(igdout)dyval
3195  WRITE(igdout)255
3196  ELSE IF(maptype==6 )THEN ! ARW rotated lat/lon grid
3197  WRITE(igdout)205
3198  WRITE(igdout)im
3199  WRITE(igdout)jm
3200  WRITE(igdout)latstart
3201  WRITE(igdout)lonstart
3202  WRITE(igdout)136
3203  WRITE(igdout)cenlat
3204  WRITE(igdout)cenlon
3205  WRITE(igdout)dxval
3206  WRITE(igdout)dyval
3207  WRITE(igdout)64
3208  WRITE(igdout)latlast
3209  WRITE(igdout)lonlast
3210  WRITE(igdout)0
3211 
3212  END IF
3213 
3214 ! following for hurricane wrf post
3215 
3216  open(10,file='copygb_hwrf.txt',form='formatted',status='unknown')
3217  idxvald = abs(lonlast-lonstart)/(im-2)
3218  idyvald = abs(latlast-latstart)/(jm-2)
3219  print*,'dxval,dyval in degree',dxval/107000.,dyval/107000.
3220  print*,'idxvald,idyvald,LATSTART,LONSTART,LATLAST,LONLAST= ', &
3221  idxvald,idyvald,latstart,lonstart,latlast,lonlast
3222  write(10,1010) im-1,jm-1,latstart,lonstart,latlast,lonlast, &
3223  idxvald,idyvald
3224 
3225 1010 format('255 0 ',2(i4,x),i8,x,i9,x,'136 ',i8,x,i9,x, &
3226  2(i8,x),'0')
3227  close (10)
3228  end if
3229 !
3230  DEALLOCATE (thv)
3231  deallocate (msft)
3232 
3233 !
3234 ! convert dxval, dyval from mtere to mm
3235 !
3236  if (grib=="grib2" )then
3237  dxval=dxval*1000.
3238  dyval=dyval*1000.
3239  endif
3240 !
3241 
3242  RETURN
3243  END
wrf_io_flags_mod
Definition: wrf_io_flags.f:1
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
zensun
subroutine zensun(day, time, lat, lon, pi, sun_zenith, sun_azimuth)
This subroutine computes solar position information as a function of geographic coordinates, date and time.
Definition: ZENSUN.f:61
soil
Definition: SOIL_mod.f:1
lookup_mod
Definition: LOOKUP.f:1