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