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