UPP (develop)
Loading...
Searching...
No Matches
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
1252if(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
1361endif ! 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
32261010 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
collect_loc
subroutine collect_loc(a, b)
COLLECT_LOC.
Definition COLLECT_LOC.f:23
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27
initpost
subroutine initpost
This routine initializes constants and variables at the start of an ETA model or post processor run.
Definition INITPOST.F:26
msfps
subroutine msfps(lat, truelat1, msf)
msfps() computes the map scale factor for a polar stereographic grid at a give latitude.
Definition MSFPS.f:27
zensun
subroutine zensun(day, time, lat, lon, pi, sun_zenith, sun_azimuth)
Makes sun zenith and sun azimuth angle.
Definition ZENSUN.f:62