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