UPP  V11.0.0
 All Data Structures Files Functions Pages
INITPOST_NETCDF.f
Go to the documentation of this file.
1 
5 
22  SUBROUTINE initpost_netcdf(ncid2d,ncid3d)
23 
24 
25  use netcdf
26  use vrbls4d, only: dust, salt, suso, soot, waso
27  use vrbls3d, only: t, q, uh, vh, pmid, pint, alpint, dpres, zint, zmid, o3, &
28  qqr, qqs, cwm, qqi, qqw, omga, rhomid, q2, cfr, rlwtt, rswtt, tcucn, &
29  tcucns, train, el_pbl, exch_h, vdifftt, vdiffmois, dconvmois, nradtt, &
30  o3vdiff, o3prod, o3tndy, mwpv, unknown, vdiffzacce, zgdrag,cnvctummixing, &
31  vdiffmacce, mgdrag, cnvctvmmixing, ncnvctcfrac, cnvctumflx, cnvctdmflx, &
32  cnvctzgdrag, sconvmois, cnvctmgdrag, cnvctdetmflx, duwt, duem, dusd, dudp, &
33  wh, qqg, ref_10cm, qqnifa, qqnwfa, pmtf, ozcon
34 
35  use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate, &
36  cprate, avgprec, prec, lspa, sno, si, cldefi, th10, q10, tshltr, pshltr, &
37  tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, landfrac, radot, sigt4, &
38  cfrach, cfracl, cfracm, avgcfrach, qshltr, avgcfracl, avgcfracm, cnvcfr, &
39  islope, cmc, grnflx, vegfrc, acfrcv, ncfrcv, acfrst, ncfrst, ssroff, &
40  bgroff, rlwin, rlwtoa, cldwork, alwin, alwout, alwtoa, rswin, rswinc, &
41  rswout, aswin, auvbin, auvbinc, aswout, aswtoa, sfcshx, sfclhx, subshx, &
42  snopcx, sfcux, sfcvx, sfcuxi, sfcvxi, sfcuvx, gtaux, gtauy, potevp, u10, v10, smstav, &
43  smstot, ivgtyp, isltyp, sfcevp, sfcexc, acsnow, acsnom, sst, thz0, qz0, &
44  uz0, vz0, ptop, htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm, pbotm, ttopm, &
45  ptoph, pboth, pblcfr, ttoph, runoff, tecan, tetran, tedir, twa, maxtshltr, &
46  mintshltr, maxrhshltr, fdnsst, &
47  minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots, &
48  cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa,rel_vort_maxhy1, &
49  maxqshltr, minqshltr, acond, sr, u10h, v10h,refd_max, w_up_max, w_dn_max, &
50  up_heli_max,up_heli_min,up_heli_max03,up_heli_min03,rel_vort_max01,u10max, v10max, &
51  avgedir,avgecan,paha,pahi,avgetrans,avgesnow,avgprec_cont,avgcprate_cont,rel_vort_max, &
52  avisbeamswin,avisdiffswin,airbeamswin,airdiffswin,refdm10c_max,wspd10max, &
53  alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, &
54  ti,aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550,prate_max, &
55  pwat
56  use soil, only: sldpth, sllevel, sh2o, smc, stc
57  use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
58  use physcons_post, only: grav => con_g, fv => con_fvirt, rgas => con_rd, &
59  eps => con_eps, epsm1 => con_epsm1
60  use params_mod, only: erad, dtr, tfrz, h1, d608, rd, p1000, capa,pi
61  use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe, &
62  ttblq, rdpq, rdtheq, stheq, the0q, the0
63  use ctlblk_mod, only: me, mpi_comm_comp, icnt, idsp, jsta, jend, ihrst, idat, sdat, ifhr, &
64  ifmin, filename, tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, td3d, restrt, sdat, &
65  jend_m, imin, imp_physics, dt, spval, pdtop, pt, qmin, nbin_du, nphs, dtq2, ardlw,&
66  ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm, &
67  jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc, &
68  nbin_oc, nbin_su, gocart_on, pt_tbl, hyb_sigp, filenameflux, filenameaer, &
69  isf_surface_physics,rdaod, aqfcmaq_on, modelname, &
70  ista, iend, ista_2l, iend_2u,iend_m
71  use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, &
72  dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, &
73  latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, standlon
74  use upp_physics, only: fpvsnew
75 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76  implicit none
77 !
78 ! type(nemsio_gfile) :: nfile,ffile,rfile
79  integer,parameter :: nvar2d=48
80 ! character(nemsio_charkind) :: name2d(nvar2d)
81  integer :: nvar3d, numdims
82 ! character(nemsio_charkind), allocatable :: name3din(:), name3dout(:)
83 ! character(nemsio_charkind) :: varname,levtype
84 !
85 ! INCLUDE/SET PARAMETERS.
86 !
87  include "mpif.h"
88 
89 ! integer,parameter:: MAXPTS=1000000 ! max im*jm points
90 !
91 ! real,parameter:: con_g =9.80665e+0! gravity
92 ! real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
93 ! real,parameter:: con_rd =2.8705e+2 ! gas constant air
94 ! real,parameter:: con_fvirt =con_rv/con_rd-1.
95 ! real,parameter:: con_eps =con_rd/con_rv
96 ! real,parameter:: con_epsm1 =con_rd/con_rv-1
97 !
98 ! This version of INITPOST shows how to initialize, open, read from, and
99 ! close a NetCDF dataset. In order to change it to read an internal (binary)
100 ! dataset, do a global replacement of _ncd_ with _int_.
101 
102  real, parameter :: gravi = 1.0/grav
103  character(len=20) :: varname, vcoordname
104  integer :: status, fldsize, fldst, recn, recn_vvel
105  character startdate*19,sysdepinfo*80,cgar*1
106  character startdate2(19)*4
107  logical :: read_lonlat=.true.
108 !
109 ! NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
110 ! AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
111 !
112 ! ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
113 ! INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
114  LOGICAL runb,singlrst,subpost,nest,hydro,ioomg,ioall
115  logical, parameter :: debugprint = .false., zerout = .false.
116 ! logical, parameter :: debugprint = .true., zerout = .false.
117  logical :: convert_rad_to_deg=.false.
118  CHARACTER*32 varcharval
119 ! CHARACTER*40 CONTRL,FILALL,FILMST,FILTMP,FILTKE,FILUNV,FILCLD,FILRAD,FILSFC
120  CHARACTER*4 resthr
121  CHARACTER fname*255,envar*50
122  INTEGER idate(8),jdate(8),jpds(200),jgds(200),kpds(200),kgds(200)
123 ! LOGICAL*1 LB(IM,JM)
124 !
125 ! INCLUDE COMMON BLOCKS.
126 !
127 ! DECLARE VARIABLES.
128 !
129 ! REAL fhour
130  integer nfhour ! forecast hour from nems io file
131  integer fhzero !bucket
132  real dtp !physics time step
133  REAL rinc(5)
134 
135  REAL dummy(im,jm)
136 !jw
137  integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus, &
138  i,j,l,ll,k,kf,irtn,igdout,n,index,nframe, &
139  nframed2,iunitd3d,ierr,idum,iret,nrec,idrt
140  integer ncid3d,ncid2d,varid,nhcas
141  real tstart,tlmh,tsph,es,fact,soilayert,soilayerb,zhour,dum, &
142  tvll,pmll,tv, tx1, tx2
143 
144  character*20,allocatable :: recname(:)
145  integer, allocatable :: reclev(:), kmsk(:,:)
146  real, allocatable :: glat1d(:), glon1d(:), qstl(:)
147  real, allocatable :: wrk1(:,:), wrk2(:,:)
148  real, allocatable :: p2d(:,:), t2d(:,:), q2d(:,:), &
149  qs2d(:,:), cw2d(:,:), cfr2d(:,:)
150  real, dimension(lm+1) :: ak5, bk5
151  real*8, allocatable :: pm2d(:,:), pi2d(:,:)
152  real, allocatable :: tmp(:)
153  real :: buf(ista_2l:iend_2u,jsta_2l:jend_2u)
154  real :: buf3d(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
155 
156 ! real buf(ista_2l:iend_2u,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u) &
157 ! ,buf3d(ista_2l:iend_2u,jsta_2l:jend_2u,lm),buf3d2(im,lp1,jsta_2l:jend_2u)
158 
159  real lat
160  integer isa, jsa, latghf, jtem, idvc, idsl, nvcoord, ip1, nn, npass
161 
162  integer, parameter :: npass2=5, npass3=30
163  real, parameter :: third=1.0/3.0
164  INTEGER, DIMENSION(2) :: ij4min, ij4max
165  REAL :: omgmin, omgmax
166  real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:)
167  REAL, ALLOCATABLE :: ps2d(:,:),psx2d(:,:),psy2d(:,:)
168  real, allocatable :: div3d(:,:,:)
169  real(kind=4),allocatable :: vcrd(:,:)
170  real :: dum_const
171 
172 ! AQF
173 
174  real, allocatable :: aacd(:,:,:), aalj(:,:,:) &
175  ,aalk1j(:,:,:), aalk2j(:,:,:) &
176  ,abnz1j(:,:,:), abnz2j(:,:,:), abnz3j(:,:,:) &
177  ,acaj(:,:,:), acet(:,:,:) &
178  ,acli(:,:,:), aclj(:,:,:), aclk(:,:,:) &
179  ,acors(:,:,:), acro_primary(:,:,:) &
180  ,acrolein(:,:,:), aeci(:,:,:) &
181  ,aecj(:,:,:), afej(:,:,:) &
182  ,aglyj(:,:,:) &
183  ,ah2oi(:,:,:), ah2oj(:,:,:), ah2ok(:,:,:) &
184  ,ah3opi(:,:,:), ah3opj(:,:,:), ah3opk(:,:,:) &
185  ,aiso1j(:,:,:), aiso2j(:,:,:), aiso3j(:,:,:) &
186  ,aivpo1j(:,:,:), akj(:,:,:) &
187  ,ald2(:,:,:), ald2_primary(:,:,:) &
188  ,aldx(:,:,:) &
189  ,alvoo1i(:,:,:), alvoo1j(:,:,:) &
190  ,alvoo2i(:,:,:), alvoo2j(:,:,:) &
191  ,alvpo1i(:,:,:), alvpo1j(:,:,:) &
192  ,amgj(:,:,:), amnj(:,:,:) &
193  ,amgk(:,:,:), akk(:,:,:), acak(:,:,:) &
194  ,anai(:,:,:), anaj(:,:,:), anak(:,:,:) &
195  ,anh4i(:,:,:), anh4j(:,:,:), anh4k(:,:,:) &
196  ,ano3i(:,:,:), ano3j(:,:,:), ano3k(:,:,:) &
197  ,aolgaj(:,:,:), aolgbj(:,:,:), aorgcj(:,:,:) &
198  ,aomi(:,:,:), aomj(:,:,:) &
199  ,aothri(:,:,:), aothrj(:,:,:) &
200  ,apah1j(:,:,:), apah2j(:,:,:), apah3j(:,:,:) &
201  ,apomi(:,:,:), apomj(:,:,:) &
202  ,apcsoj(:,:,:), aseacat(:,:,:), asij(:,:,:) &
203  ,aso4i(:,:,:), aso4j(:,:,:), aso4k(:,:,:) &
204  ,asoil(:,:,:), asqtj(:,:,:) &
205  ,asomi(:,:,:), asomj(:,:,:) &
206  ,asvoo1i(:,:,:), asvoo1j(:,:,:) &
207  ,asvoo2i(:,:,:), asvoo2j(:,:,:) &
208  ,asvoo3j(:,:,:) &
209  ,asvpo1i(:,:,:), asvpo1j(:,:,:) &
210  ,asvpo2i(:,:,:), asvpo2j(:,:,:) &
211  ,asvpo3j(:,:,:) &
212  ,atij(:,:,:) &
213  ,atol1j(:,:,:), atol2j(:,:,:), atol3j(:,:,:) &
214  ,atoti(:,:,:), atotj(:,:,:), atotk(:,:,:) &
215  ,atrp1j(:,:,:), atrp2j(:,:,:) &
216  ,axyl1j(:,:,:), axyl2j(:,:,:), axyl3j(:,:,:) &
217  ,pm25ac(:,:,:), pm25at(:,:,:), pm25co(:,:,:)
218 
219 
220  if (me == 0) print *,' aqfcmaq_on=', aqfcmaq_on
221 
222  if (aqfcmaq_on) then
223 
224  allocate(aacd(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
225  allocate(aalj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
226  allocate(aalk1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
227  allocate(aalk2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
228 
229  allocate(abnz1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
230  allocate(abnz2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
231  allocate(abnz3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
232 
233  allocate(acaj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
234  allocate(acet(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
235 
236  allocate(acli(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
237  allocate(aclj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
238  allocate(aclk(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
239 
240  allocate(acors(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
241  allocate(acro_primary(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
242  allocate(acrolein(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
243  allocate(aeci(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
244  allocate(aecj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
245  allocate(afej(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
246  allocate(aglyj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
247 
248  allocate(ah2oi(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
249  allocate(ah2oj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
250  allocate(ah2ok(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
251 
252  allocate(ah3opi(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
253  allocate(ah3opj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
254  allocate(ah3opk(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
255 
256  allocate(aiso1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
257  allocate(aiso2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
258  allocate(aiso3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
259 
260  allocate(aivpo1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
261  allocate(akj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
262 
263  allocate(ald2(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
264  allocate(ald2_primary(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
265 
266  allocate(aldx(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
267 
268  allocate(alvoo1i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
269  allocate(alvoo1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
270  allocate(alvoo2i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
271  allocate(alvoo2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
272  allocate(alvpo1i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
273  allocate(alvpo1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
274 
275  allocate(amgj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
276  allocate(amnj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
277 
278  allocate(anai(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
279  allocate(anaj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
280  allocate(anak(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
281 
282  allocate(anh4i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
283  allocate(anh4j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
284  allocate(anh4k(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
285 
286  allocate(ano3i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
287  allocate(ano3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
288  allocate(ano3k(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
289 
290  allocate(aolgaj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
291  allocate(aolgbj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
292 
293  allocate(aomi(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
294  allocate(aomj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
295 
296  allocate(aorgcj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
297 
298  allocate(aothri(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
299  allocate(aothrj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
300 
301  allocate(apah1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
302  allocate(apah2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
303  allocate(apah3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
304 
305  allocate(apcsoj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
306 
307  allocate(apomi(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
308  allocate(apomj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
309 
310  allocate(aseacat(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
311  allocate(asij(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
312 
313  allocate(aso4i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
314  allocate(aso4j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
315  allocate(aso4k(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
316  allocate(asoil(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
317 
318  allocate(asomi(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
319  allocate(asomj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
320 
321  allocate(asqtj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
322 
323  allocate(asvoo1i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
324  allocate(asvoo1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
325  allocate(asvoo2i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
326  allocate(asvoo2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
327  allocate(asvoo3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
328 
329  allocate(asvpo1i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
330  allocate(asvpo1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
331  allocate(asvpo2i(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
332  allocate(asvpo2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
333  allocate(asvpo3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
334 
335  allocate(atij(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
336 
337  allocate(atol1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
338  allocate(atol2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
339  allocate(atol3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
340 
341  allocate(atoti(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
342  allocate(atotj(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
343  allocate(atotk(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
344 
345  allocate(atrp1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
346  allocate(atrp2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
347 
348  allocate(axyl1j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
349  allocate(axyl2j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
350  allocate(axyl3j(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
351 
352  allocate(pm25ac(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
353  allocate(pm25at(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
354  allocate(pm25co(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
355 
356  endif
357 
358 !***********************************************************************
359 ! START INIT HERE.
360 !
361  WRITE(6,*)'INITPOST: ENTER INITPOST_NETCDF'
362  WRITE(6,*)'me=',me, &
363  'jsta_2l=',jsta_2l,'jend_2u=', &
364  jend_2u,'im=',im, &
365  'ista_2l=',ista_2l,'iend_2u=',iend_2u, &
366  'ista=',ista,'iend=',iend, &
367  'iend_m=',iend_m
368 !
369  isa = (ista+iend) / 2
370  jsa = (jsta+jend) / 2
371 
372 !$omp parallel do private(i,j)
373  do j = jsta_2l, jend_2u
374  do i= ista_2l, iend_2u
375  buf(i,j) = spval
376  enddo
377  enddo
378 
379  status=nf90_get_att(ncid3d,nf90_global,'ak',ak5)
380  if(status/=0)then
381  print*,'ak not found; assigning missing value'
382  ak5=spval
383  else
384  if(me==0)print*,'ak5= ',ak5
385  end if
386  status=nf90_get_att(ncid3d,nf90_global,'idrt',idrt)
387  if(status/=0)then
388  print*,'idrt not in netcdf file,reading grid'
389  status=nf90_get_att(ncid3d,nf90_global,'grid',varcharval)
390  if(status/=0)then
391  print*,'idrt and grid not in netcdf file, set default to latlon'
392  idrt=0
393  maptype=0
394  else
395  if(trim(varcharval)=='rotated_latlon')then
396  maptype=207
397  idrt=207
398  status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const)
399  if(status/=0)then
400  print*,'cen_lon not found; assigning missing value'
401  cenlon=spval
402  else
403  if(dum_const<0.)then
404  cenlon=nint((dum_const+360.)*gdsdegr)
405  else
406  cenlon=dum_const*gdsdegr
407  end if
408  end if
409  status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const)
410  if(status/=0)then
411  print*,'cen_lat not found; assigning missing value'
412  cenlat=spval
413  else
414  cenlat=dum_const*gdsdegr
415  end if
416 
417  status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
418  if(status/=0)then
419  print*,'lonstart_r not found; assigning missing value'
420  lonstart_r=spval
421  else
422  if(dum_const<0.)then
423  lonstart_r=nint((dum_const+360.)*gdsdegr)
424  else
425  lonstart_r=dum_const*gdsdegr
426  end if
427  end if
428  status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
429  if(status/=0)then
430  print*,'latstart_r not found; assigning missing value'
431  latstart_r=spval
432  else
433  latstart_r=dum_const*gdsdegr
434  end if
435 
436  status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const)
437  if(status/=0)then
438  print*,'lonlast_r not found; assigning missing value'
439  lonlast_r=spval
440  else
441  if(dum_const<0.)then
442  lonlast_r=nint((dum_const+360.)*gdsdegr)
443  else
444  lonlast_r=dum_const*gdsdegr
445  end if
446  end if
447  status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const)
448  if(status/=0)then
449  print*,'latlast_r not found; assigning missing value'
450  latlast_r=spval
451  else
452  latlast_r=dum_const*gdsdegr
453  end if
454 
455  status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const)
456  if(status/=0)then
457  print*,'dlmd not found; assigning missing value'
458  dxval=spval
459  else
460  dxval=dum_const*gdsdegr
461  end if
462  status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const)
463  if(status/=0)then
464  print*,'dphd not found; assigning missing value'
465  dyval=spval
466  else
467  dyval=dum_const*gdsdegr
468  end if
469 
470  print*,'lonstart,latstart,cenlon,cenlat,dyval,dxval', &
471  lonstart,latstart,cenlon,cenlat,dyval,dxval
472 
473 ! Jili Dong add support for regular lat lon (2019/03/22) start
474  else if(trim(varcharval)=='latlon')then
475  maptype=0
476  idrt=0
477 
478  status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
479  if(status/=0)then
480  print*,'lonstart not found; assigning missing value'
481  lonstart=spval
482  else
483  if(dum_const<0.)then
484  lonstart=nint((dum_const+360.)*gdsdegr)
485  else
486  lonstart=dum_const*gdsdegr
487  end if
488  end if
489  status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
490  if(status/=0)then
491  print*,'latstart not found; assigning missing value'
492  latstart=spval
493  else
494  latstart=dum_const*gdsdegr
495  end if
496 
497  status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const)
498  if(status/=0)then
499  print*,'lonlast not found; assigning missing value'
500  lonlast=spval
501  else
502  if(dum_const<0.)then
503  lonlast=nint((dum_const+360.)*gdsdegr)
504  else
505  lonlast=dum_const*gdsdegr
506  end if
507  end if
508  status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const)
509  if(status/=0)then
510  print*,'latlast not found; assigning missing value'
511  latlast=spval
512  else
513  latlast=dum_const*gdsdegr
514  end if
515 
516  status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const)
517  if(status/=0)then
518  print*,'dlmd not found; assigning missing value'
519  dxval=spval
520  else
521  dxval=dum_const*gdsdegr
522  end if
523  status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const)
524  if(status/=0)then
525  print*,'dphd not found; assigning missing value'
526  dyval=spval
527  else
528  dyval=dum_const*gdsdegr
529  end if
530 
531  print*,'lonstart,latstart,dyval,dxval', &
532  lonstart,lonlast,latstart,latlast,dyval,dxval
533 
534 ! Jili Dong add support for regular lat lon (2019/03/22) end
535 
536  ELSE IF (trim(varcharval)=='lambert_conformal')then
537 
538  maptype=1
539  idrt=1
540  status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const)
541  if(status/=0)then
542  print*,'cen_lon not found; assigning missing value'
543  cenlon=spval
544  else
545  if(dum_const<0.)then
546  cenlon=nint((dum_const+360.)*gdsdegr)
547  else
548  cenlon=dum_const*gdsdegr
549  end if
550  end if
551  status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const)
552  if(status/=0)then
553  print*,'cen_lat not found; assigning missing value'
554  cenlat=spval
555  else
556  cenlat=dum_const*gdsdegr
557  end if
558 
559  status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const)
560  if(status/=0)then
561  print*,'lonstart not found; assigning missing value'
562  lonstart=spval
563  else
564  if(dum_const<0.)then
565  lonstart=nint((dum_const+360.)*gdsdegr)
566  else
567  lonstart=dum_const*gdsdegr
568  end if
569  end if
570  status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const)
571  if(status/=0)then
572  print*,'latstart not found; assigning missing value'
573  latstart=spval
574  else
575  latstart=dum_const*gdsdegr
576  end if
577 
578  status=nf90_get_att(ncid3d,nf90_global,'stdlat1',dum_const)
579  if(status/=0)then
580  print*,'stdlat1 not found; assigning missing value'
581  truelat1=spval
582  else
583  truelat1=dum_const*gdsdegr
584  end if
585  status=nf90_get_att(ncid3d,nf90_global,'stdlat2',dum_const)
586  if(status/=0)then
587  print*,'stdlat2 not found; assigning missing value'
588  truelat2=spval
589  else
590  truelat2=dum_const*gdsdegr
591  end if
592 
593  status=nf90_get_att(ncid3d,nf90_global,'dx',dum_const)
594  if(status/=0)then
595  print*,'dx not found; assigning missing value'
596  dxval=spval
597  else
598  dxval=dum_const*1.e3
599  end if
600  status=nf90_get_att(ncid3d,nf90_global,'dy',dum_const)
601  if(status/=0)then
602  print*,'dphd not found; assigning missing value'
603  dyval=spval
604  else
605  dyval=dum_const*1.e3
606  end if
607 
608  standlon = cenlon
609  print*,'lonstart,latstart,cenlon,cenlat,truelat1,truelat2, &
610  stadlon,dyval,dxval', &
611  lonstart,latstart,cenlon,cenlat,truelat1,truelat2,standlon,dyval,dxval
612 
613  else if(trim(varcharval)=='gaussian')then
614  maptype=4
615  idrt=4
616  else ! setting default maptype
617  maptype=0
618  idrt=0
619  end if
620  end if
621  end if
622  if(me==0)print*,'idrt MAPTYPE= ',idrt,maptype
623 ! STEP 1. READ MODEL OUTPUT FILE
624 !
625 !
626 !***
627 !
628 ! LMH and LMV always = LM for sigma-type vert coord
629 
630 !$omp parallel do private(i,j)
631  do j = jsta_2l, jend_2u
632  do i = ista_2l, iend_2u
633  lmv(i,j) = lm
634  lmh(i,j) = lm
635  end do
636  end do
637 
638 ! HTM VTM all 1 for sigma-type vert coord
639 
640 !$omp parallel do private(i,j,l)
641  do l = 1, lm
642  do j = jsta_2l, jend_2u
643  do i = ista_2l, iend_2u
644  htm(i,j,l) = 1.0
645  vtm(i,j,l) = 1.0
646  end do
647  end do
648  end do
649 
650  status=nf90_get_att(ncid3d,nf90_global,'nhcas',nhcas)
651  if(status/=0)then
652  print*,'nhcas not in netcdf file, set default to nonhydro'
653  nhcas=0
654  end if
655  if(me==0)print*,'nhcas= ',nhcas
656  if (nhcas == 0 ) then !non-hydrostatic case
657  nrec=14
658  allocate (recname(nrec))
659  recname=[character(len=20) :: 'ugrd','vgrd','spfh','tmp','o3mr', &
660  'presnh','dzdt', 'clwmr','dpres', &
661  'delz','icmr','rwmr', &
662  'snmr','grle']
663  else
664  nrec=8
665  allocate (recname(nrec))
666  recname=[character(len=20) :: 'ugrd','vgrd','tmp','spfh','o3mr', &
667  'hypres', 'clwmr','dpres']
668  endif
669 
670 ! write(0,*)'nrec=',nrec
671  !allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
672  allocate(glat1d(jm),glon1d(im))
673 
674 ! hardwire idate for now
675 ! idate=(/2017,08,07,00,0,0,0,0/)
676 ! get cycle start time
677  status=nf90_inq_varid(ncid3d,'time',varid)
678  if(status/=0)then
679  print*,'time not in netcdf file, stopping'
680  stop 1
681  else
682  status=nf90_get_att(ncid3d,varid,'units',varcharval)
683  if(status/=0)then
684  print*,'time unit not available'
685  else
686  print*,'time unit read from netcdf file= ',varcharval
687 ! assume use hours as unit
688 ! idate_loc=index(varcharval,'since')+6
689  read(varcharval,101)idate(1),idate(2),idate(3),idate(4),idate(5)
690  end if
691 ! Status=nf90_inquire_dimension(ncid3d,varid,len=ntimes)
692 ! allocate(fhours(ntimes))
693 ! status = nf90_inq_varid(ncid3d,varid,fhours)
694 ! Status=nf90_get_var(ncid3d,varid,nfhour,start=(/1/), &
695 ! count=(/1/))
696 ! if(Status/=0)then
697 ! print*,'forecast hour not in netcdf file, stopping'
698 ! STOP 1
699 ! end if
700  end if
701  101 format(t13,i4,1x,i2,1x,i2,1x,i2,1x,i2)
702  print*,'idate= ',idate(1:5)
703 
704 ! Jili Dong check output format for coordinate reading
705  status=nf90_inq_varid(ncid3d,'grid_xt',varid)
706  status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
707  if(numdims==1.and.modelname=="FV3R") then
708  read_lonlat=.true.
709  else
710  read_lonlat=.false.
711  end if
712 
713 
714 ! Jili Dong add support for new write component output
715 ! get longitude
716  if (read_lonlat) then
717  status=nf90_inq_varid(ncid3d,'lon',varid)
718  status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
719  if(debugprint)print*,'number of dim for gdlon ',numdims
720  else
721  status=nf90_inq_varid(ncid3d,'grid_xt',varid)
722  status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
723  if(debugprint)print*,'number of dim for gdlon ',numdims
724  end if
725  if(numdims==1)then
726  status=nf90_get_var(ncid3d,varid,glon1d)
727  do j=jsta,jend
728  do i=ista,iend
729  gdlon(i,j) = real(glon1d(i),kind=4)
730  end do
731  end do
732  lonstart = nint(glon1d(1)*gdsdegr)
733  lonlast = nint(glon1d(im)*gdsdegr)
734 
735 ! Jili Dong add support for regular lat lon (2019/03/22) start
736  if (maptype == 0) then
737  if(lonstart<0.)then
738  lonstart=lonstart+360.*gdsdegr
739  end if
740  if(lonlast<0.)then
741  lonlast=lonlast+360.*gdsdegr
742  end if
743  end if
744 ! Jili Dong add support for regular lat lon (2019/03/22) end
745 
746  else if(numdims==2)then
747  status=nf90_get_var(ncid3d,varid,dummy)
748  if(maxval(abs(dummy))<2.0*pi)convert_rad_to_deg=.true.
749  if(convert_rad_to_deg)then
750  do j=jsta,jend
751  do i=ista,iend
752  gdlon(i,j) = real(dummy(i,j),kind=4)*180./pi
753  end do
754  end do
755  else
756  do j=jsta,jend
757  do i=ista,iend
758  gdlon(i,j) = real(dummy(i,j),kind=4)
759  end do
760  end do
761  end if
762  if(convert_rad_to_deg)then
763  lonstart = nint(dummy(1,1)*gdsdegr)*180./pi
764  lonlast = nint(dummy(im,jm)*gdsdegr)*180./pi
765  else
766  lonstart = nint(dummy(1,1)*gdsdegr)
767  lonlast = nint(dummy(im,jm)*gdsdegr)
768  end if
769 
770 ! Jili Dong add support for regular lat lon (2019/03/22) start
771  if (maptype == 0) then
772  if(lonstart<0.)then
773  lonstart=lonstart+360.*gdsdegr
774  end if
775  if(lonlast<0.)then
776  lonlast=lonlast+360.*gdsdegr
777  end if
778  end if
779 ! Jili Dong add support for regular lat lon (2019/03/22) end
780 
781  end if
782  print*,'lonstart,lonlast ',lonstart,lonlast
783 ! Jili Dong add support for new write component output
784 ! get latitude
785  if (read_lonlat) then
786  status=nf90_inq_varid(ncid3d,'lat',varid)
787  status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
788  if(debugprint)print*,'number of dim for gdlat ',numdims
789  else
790  status=nf90_inq_varid(ncid3d,'grid_yt',varid)
791  status=nf90_inquire_variable(ncid3d,varid,ndims = numdims)
792  if(debugprint)print*,'number of dim for gdlat ',numdims
793  end if
794  if(numdims==1)then
795  status=nf90_get_var(ncid3d,varid,glat1d)
796  do j=jsta,jend
797  do i=ista,iend
798  gdlat(i,j) = real(glat1d(j),kind=4)
799  end do
800  end do
801  latstart = nint(glat1d(1)*gdsdegr)
802  latlast = nint(glat1d(jm)*gdsdegr)
803  else if(numdims==2)then
804  status=nf90_get_var(ncid3d,varid,dummy)
805  if(maxval(abs(dummy))<pi)convert_rad_to_deg=.true.
806  if(convert_rad_to_deg)then
807  do j=jsta,jend
808  do i=ista,iend
809  gdlat(i,j) = real(dummy(i,j),kind=4)*180./pi
810  end do
811  end do
812  else
813  do j=jsta,jend
814  do i=ista,iend
815  gdlat(i,j) = real(dummy(i,j),kind=4)
816  end do
817  end do
818  end if
819  if(convert_rad_to_deg)then
820  latstart = nint(dummy(1,1)*gdsdegr)*180./pi
821  latlast = nint(dummy(im,jm)*gdsdegr)*180./pi
822  else
823  latstart = nint(dummy(1,1)*gdsdegr)
824  latlast = nint(dummy(im,jm)*gdsdegr)
825  end if
826  end if
827  print*,'laststart,latlast = ',latstart,latlast
828  if(debugprint)print*,'me sample gdlon gdlat= ' &
829  ,me,gdlon(isa,jsa),gdlat(isa,jsa)
830 
831 ! Specify grid staggering type
832  gridtype = 'A'
833  maptype=idrt
834  if (me == 0) print *, 'maptype and gridtype is ', &
835  maptype,gridtype
836 
837  if(gridtype == 'A')then
838  lonstartv=lonstart
839  lonlastv=lonlast
840  latstartv=latstart
841  latlastv=latlast
842  cenlatv=cenlat
843  cenlonv=cenlon
844  end if
845 
846  if(debugprint)then
847  if (me == 0)then
848  do i=1,nrec
849  print *,'recname=',trim(recname(i))
850  end do
851  end if
852  end if
853 
854 
855  deallocate(glat1d,glon1d)
856 
857  print*,'idate = ',(idate(i),i=1,7)
858 ! print*,'nfhour = ',nfhour
859 
860 ! sample print point
861  ii = im/2
862  jj = jm/2
863 
864  print *,me,'max(gdlat)=', maxval(gdlat), &
865  'max(gdlon)=', maxval(gdlon)
866  CALL exch(gdlat(ista_2l,jsta_2l))
867  CALL exch(gdlon(ista_2l,jsta_2l))
868  print *,'after call EXCH,me=',me
869 
870 !$omp parallel do private(i,j,ip1)
871  do j = jsta, jend_m
872  do i = ista, iend_m
873  ip1 = i + 1
874 ! if (ip1 > im) ip1 = ip1 - im
875  dx(i,j) = erad*cos(gdlat(i,j)*dtr) *(gdlon(ip1,j)-gdlon(i,j))*dtr
876  dy(i,j) = erad*(gdlat(i,j+1)-gdlat(i,j))*dtr ! like A*DPH
877 ! F(I,J)=1.454441e-4*sin(gdlat(i,j)*DTR) ! 2*omeg*sin(phi)
878 ! if (i == ii .and. j == jj) print*,'sample LATLON, DY, DY=' &
879 ! ,i,j,GDLAT(I,J),GDLON(I,J),DX(I,J),DY(I,J)
880  end do
881  end do
882  if(debugprint)print*,'me sample dx dy= ' &
883  ,me,dx(isa,jsa),dy(isa,jsa)
884 !$omp parallel do private(i,j)
885  do j=jsta,jend
886  do i=ista,iend
887  f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr) ! 2*omeg*sin(phi)
888  end do
889  end do
890 
891  iyear = idate(1)
892  imn = idate(2)
893  iday = idate(3)
894  ihrst = idate(4)
895  imin = idate(5)
896  jdate = 0
897  idate = 0
898 !
899  print*,'start yr mo day hr min =',iyear,imn,iday,ihrst,imin
900  print*,'processing yr mo day hr min=' &
901  ,idat(3),idat(1),idat(2),idat(4),idat(5)
902 !
903  idate(1) = iyear
904  idate(2) = imn
905  idate(3) = iday
906  idate(5) = ihrst
907  idate(6) = imin
908  sdat(1) = imn
909  sdat(2) = iday
910  sdat(3) = iyear
911  jdate(1) = idat(3)
912  jdate(2) = idat(1)
913  jdate(3) = idat(2)
914  jdate(5) = idat(4)
915  jdate(6) = idat(5)
916 !
917  print *,' idate=',idate
918  print *,' jdate=',jdate
919 !
920  CALL w3difdat(jdate,idate,0,rinc)
921 !
922  print *,' rinc=',rinc
923  ifhr = nint(rinc(2)+rinc(1)*24.)
924  print *,' ifhr=',ifhr
925  ifmin = nint(rinc(3))
926 ! if(ifhr /= nint(fhour))print*,'find wrong Grib file';stop
927  print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,filename
928 
929 ! Getting tstart
930  tstart = 0.
931  print*,'tstart= ',tstart
932 
933 ! Getiing restart
934 
935  restrt = .true. ! set RESTRT as default
936 
937  IF(tstart > 1.0e-2)THEN
938  ifhr = ifhr+nint(tstart)
939  rinc = 0
940  idate = 0
941  rinc(2) = -1.0*ifhr
942  call w3movdat(rinc,jdate,idate)
943  sdat(1) = idate(2)
944  sdat(2) = idate(3)
945  sdat(3) = idate(1)
946  ihrst = idate(5)
947  print*,'new forecast hours for restrt run= ',ifhr
948  print*,'new start yr mo day hr min =',sdat(3),sdat(1) &
949  ,sdat(2),ihrst,imin
950  END IF
951 
952 ! GFS does not need DT to compute accumulated fields, set it to one
953 ! VarName='DT'
954  dt = 1
955 
956  hbm2 = 1.0
957 
958 ! start reading 3d netcdf output
959  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
960  spval,recname(1),uh(ista_2l,jsta_2l,1),lm)
961  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
962  spval,recname(2),vh(ista_2l,jsta_2l,1),lm)
963  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
964  spval,recname(3),q(ista_2l,jsta_2l,1),lm)
965  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
966  spval,recname(4),t(ista_2l,jsta_2l,1),lm)
967  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
968  spval,recname(5),o3(ista_2l,jsta_2l,1),lm)
969  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
970  spval,recname(7),wh(ista_2l,jsta_2l,1),lm)
971  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
972  spval,recname(8),qqw(ista_2l,jsta_2l,1),lm)
973  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
974  spval,recname(9),dpres(ista_2l,jsta_2l,1),lm)
975  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
976  spval,recname(10),buf3d(ista_2l,jsta_2l,1),lm)
977  do l=1,lm
978  do j=jsta,jend
979  do i=ista,iend
980  cwm(i,j,l)=spval
981 ! dong add missing value
982  if (wh(i,j,l) < spval) then
983  omga(i,j,l)=(-1.)*wh(i,j,l)*dpres(i,j,l)/abs(buf3d(i,j,l))
984  else
985  omga(i,j,l) = spval
986  end if
987 ! if(t(i,j,l)>1000.)print*,'bad T ',t(i,j,l)
988  enddo
989  enddo
990  enddo
991  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
992  spval,recname(11),qqi(ista_2l,jsta_2l,1),lm)
993  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
994  spval,recname(12),qqr(ista_2l,jsta_2l,1),lm)
995  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
996  spval,recname(13),qqs(ista_2l,jsta_2l,1),lm)
997  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
998  spval,recname(14),qqg(ista_2l,jsta_2l,1),lm)
999 
1000 ! calculate CWM from FV3 output
1001  do l=1,lm
1002  do j=jsta,jend
1003  do i=ista,iend
1004  cwm(i,j,l)=qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l)
1005  enddo
1006  enddo
1007  if(debugprint)print*,'sample l,t,q,u,v,w= ',isa,jsa,l &
1008  ,t(isa,jsa,l),q(isa,jsa,l),uh(isa,jsa,l),vh(isa,jsa,l) &
1009  ,wh(isa,jsa,l)
1010  if(debugprint)print*,'sample l cwm for FV3',l, &
1011  cwm(isa,jsa,l)
1012  end do
1013 
1014 ! instantaneous 3D cloud fraction
1015  if ( imp_physics==11) then !GFDL MP
1016  varname='cld_amt'
1017  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1018  spval,varname,cfr(ista_2l,jsta_2l,1),lm)
1019  else
1020  varname='cldfra'
1021  call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1022  spval,varname,cfr(ista_2l,jsta_2l,1),lm)
1023  endif
1024 ! do l=1,lm
1025 ! if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' &
1026 ! ,cfr(isa,jsa,l),isa,jsa,l
1027 ! enddo
1028 
1029 !=============================
1030 ! For AQF Chemical species
1031 !=============================
1032 
1033  if (aqfcmaq_on) then
1034 
1035  ! *********** VarName need to be in lower case ************
1036  ! === It will cause problem if not use the lower case =====
1037  ! *********************************************************
1038 
1039  !--------------------------------------------------------------
1040  !-- rename input o3 to NCO grib2 name ozcon -------------------
1041 
1042  varname='o3'
1043  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1044  spval,varname,ozcon(ista_2l,jsta_2l,1),lm)
1045 
1046  !--------------------------------------------------------------
1047 
1048  varname='aacd'
1049  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1050  spval,varname,aacd(ista_2l,jsta_2l,1),lm)
1051 
1052  varname='aalj'
1053  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1054  spval,varname,aalj(ista_2l,jsta_2l,1),lm)
1055 
1056  varname='aalk1j'
1057  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1058  spval,varname,aalk1j(ista_2l,jsta_2l,1),lm)
1059 
1060  varname='aalk2j'
1061  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1062  spval,varname,aalk2j(ista_2l,jsta_2l,1),lm)
1063 
1064  varname='abnz1j'
1065  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1066  spval,varname,abnz1j(ista_2l,jsta_2l,1),lm)
1067 
1068  varname='abnz2j'
1069  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1070  spval,varname,abnz2j(ista_2l,jsta_2l,1),lm)
1071 
1072  varname='abnz3j'
1073  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1074  spval,varname,abnz3j(ista_2l,jsta_2l,1),lm)
1075 
1076  varname='acaj'
1077  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1078  spval,varname,acaj(ista_2l,jsta_2l,1),lm)
1079 
1080  varname='acet'
1081  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1082  spval,varname,acet(ista_2l,jsta_2l,1),lm)
1083 
1084  varname='acli'
1085  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1086  spval,varname,acli(ista_2l,jsta_2l,1),lm)
1087 
1088  varname='aclj'
1089  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1090  spval,varname,aclj(ista_2l,jsta_2l,1),lm)
1091 
1092  varname='aclk'
1093  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1094  spval,varname,aclk(ista_2l,jsta_2l,1),lm)
1095 
1096  varname='acors'
1097  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1098  spval,varname,acors(ista_2l,jsta_2l,1),lm)
1099 
1100  varname='acro_primary'
1101  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1102  spval,varname,acro_primary(ista_2l,jsta_2l,1),lm)
1103 
1104  varname='acrolein'
1105  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1106  spval,varname,acrolein(ista_2l,jsta_2l,1),lm)
1107 
1108  varname='aeci'
1109  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1110  spval,varname,aeci(ista_2l,jsta_2l,1),lm)
1111 
1112  varname='aecj'
1113  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1114  spval,varname,aecj(ista_2l,jsta_2l,1),lm)
1115 
1116  varname='afej'
1117  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1118  spval,varname,afej(ista_2l,jsta_2l,1),lm)
1119 
1120  varname='aglyj'
1121  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1122  spval,varname,aglyj(ista_2l,jsta_2l,1),lm)
1123 
1124  varname='ah2oi'
1125  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1126  spval,varname,ah2oi(ista_2l,jsta_2l,1),lm)
1127 
1128  varname='ah2oj'
1129  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1130  spval,varname,ah2oj(ista_2l,jsta_2l,1),lm)
1131 
1132  varname='ah2ok'
1133  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1134  spval,varname,ah2ok(ista_2l,jsta_2l,1),lm)
1135 
1136  varname='ah3opi'
1137  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1138  spval,varname,ah3opi(ista_2l,jsta_2l,1),lm)
1139 
1140  varname='ah3opj'
1141  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1142  spval,varname,ah3opj(ista_2l,jsta_2l,1),lm)
1143 
1144  varname='ah3opk'
1145  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1146  spval,varname,ah3opk(ista_2l,jsta_2l,1),lm)
1147 
1148  varname='aiso1j'
1149  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1150  spval,varname,aiso1j(ista_2l,jsta_2l,1),lm)
1151 
1152  varname='aiso2j'
1153  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1154  spval,varname,aiso2j(ista_2l,jsta_2l,1),lm)
1155 
1156  varname='aiso3j'
1157  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1158  spval,varname,aiso3j(ista_2l,jsta_2l,1),lm)
1159 
1160  varname='aivpo1j'
1161  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1162  spval,varname,aivpo1j(ista_2l,jsta_2l,1),lm)
1163 
1164  varname='akj'
1165  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1166  spval,varname,akj(ista_2l,jsta_2l,1),lm)
1167 
1168  varname='ald2'
1169  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1170  spval,varname,ald2(ista_2l,jsta_2l,1),lm)
1171 
1172  varname='ald2_primary'
1173  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1174  spval,varname,ald2_primary(ista_2l,jsta_2l,1),lm)
1175 
1176  varname='aldx'
1177  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1178  spval,varname,aldx(ista_2l,jsta_2l,1),lm)
1179 
1180  varname='alvoo1i'
1181  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1182  spval,varname,alvoo1i(ista_2l,jsta_2l,1),lm)
1183 
1184  varname='alvoo1j'
1185  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1186  spval,varname,alvoo1j(ista_2l,jsta_2l,1),lm)
1187 
1188  varname='alvoo2i'
1189  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1190  spval,varname,alvoo2i(ista_2l,jsta_2l,1),lm)
1191 
1192  varname='alvoo2j'
1193  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1194  spval,varname,alvoo2j(ista_2l,jsta_2l,1),lm)
1195 
1196  varname='alvpo1i'
1197  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1198  spval,varname,alvpo1i(ista_2l,jsta_2l,1),lm)
1199 
1200  varname='alvpo1j'
1201  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1202  spval,varname,alvpo1j(ista_2l,jsta_2l,1),lm)
1203 
1204  varname='amgj'
1205  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1206  spval,varname,amgj(ista_2l,jsta_2l,1),lm)
1207 
1208  varname='amnj'
1209  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1210  spval,varname,amnj(ista_2l,jsta_2l,1),lm)
1211 
1212  varname='anai'
1213  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1214  spval,varname,anai(ista_2l,jsta_2l,1),lm)
1215 
1216  varname='anaj'
1217  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1218  spval,varname,anaj(ista_2l,jsta_2l,1),lm)
1219 
1220  varname='anh4i'
1221  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1222  spval,varname,anh4i(ista_2l,jsta_2l,1),lm)
1223 
1224  varname='anh4j'
1225  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1226  spval,varname,anh4j(ista_2l,jsta_2l,1),lm)
1227 
1228  varname='anh4k'
1229  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1230  spval,varname,anh4k(ista_2l,jsta_2l,1),lm)
1231 
1232  varname='ano3i'
1233  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1234  spval,varname,ano3i(ista_2l,jsta_2l,1),lm)
1235 
1236  varname='ano3j'
1237  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1238  spval,varname,ano3j(ista_2l,jsta_2l,1),lm)
1239 
1240  varname='ano3k'
1241  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1242  spval,varname,ano3k(ista_2l,jsta_2l,1),lm)
1243 
1244  varname='aolgaj'
1245  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1246  spval,varname,aolgaj(ista_2l,jsta_2l,1),lm)
1247 
1248  varname='aolgbj'
1249  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1250  spval,varname,aolgbj(ista_2l,jsta_2l,1),lm)
1251 
1252  varname='aorgcj'
1253  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1254  spval,varname,aorgcj(ista_2l,jsta_2l,1),lm)
1255 
1256  varname='aothri'
1257  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1258  spval,varname,aothri(ista_2l,jsta_2l,1),lm)
1259 
1260  varname='aothrj'
1261  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1262  spval,varname,aothrj(ista_2l,jsta_2l,1),lm)
1263 
1264  varname='apah1j'
1265  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1266  spval,varname,apah1j(ista_2l,jsta_2l,1),lm)
1267 
1268  varname='apah2j'
1269  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1270  spval,varname,apah2j(ista_2l,jsta_2l,1),lm)
1271 
1272  varname='apah3j'
1273  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1274  spval,varname,apah3j(ista_2l,jsta_2l,1),lm)
1275 
1276  varname='apcsoj'
1277  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1278  spval,varname,apcsoj(ista_2l,jsta_2l,1),lm)
1279 
1280  varname='aseacat'
1281  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1282  spval,varname,aseacat(ista_2l,jsta_2l,1),lm)
1283 
1284  varname='asij'
1285  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1286  spval,varname,asij(ista_2l,jsta_2l,1),lm)
1287 
1288  varname='aso4i'
1289  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1290  spval,varname,aso4i(ista_2l,jsta_2l,1),lm)
1291 
1292  varname='aso4j'
1293  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1294  spval,varname,aso4j(ista_2l,jsta_2l,1),lm)
1295 
1296  varname='aso4k'
1297  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1298  spval,varname,aso4k(ista_2l,jsta_2l,1),lm)
1299 
1300  varname='asoil'
1301  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1302  spval,varname,asoil(ista_2l,jsta_2l,1),lm)
1303 
1304  varname='asqtj'
1305  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1306  spval,varname,asqtj(ista_2l,jsta_2l,1),lm)
1307 
1308  varname='asvoo1i'
1309  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1310  spval,varname,asvoo1i(ista_2l,jsta_2l,1),lm)
1311 
1312  varname='asvoo1j'
1313  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1314  spval,varname,asvoo1j(ista_2l,jsta_2l,1),lm)
1315 
1316  varname='asvoo2i'
1317  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1318  spval,varname,asvoo2i(ista_2l,jsta_2l,1),lm)
1319 
1320  varname='asvoo2j'
1321  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1322  spval,varname,asvoo2j(ista_2l,jsta_2l,1),lm)
1323 
1324  varname='asvoo3j'
1325  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1326  spval,varname,asvoo3j(ista_2l,jsta_2l,1),lm)
1327 
1328  varname='asvpo1i'
1329  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1330  spval,varname,asvpo1i(ista_2l,jsta_2l,1),lm)
1331 
1332  varname='asvpo1j'
1333  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1334  spval,varname,asvpo1j(ista_2l,jsta_2l,1),lm)
1335 
1336  varname='asvpo2i'
1337  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1338  spval,varname,asvpo2i(ista_2l,jsta_2l,1),lm)
1339 
1340  varname='asvpo2j'
1341  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1342  spval,varname,asvpo2j(ista_2l,jsta_2l,1),lm)
1343 
1344  varname='asvpo3j'
1345  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1346  spval,varname,asvpo3j(ista_2l,jsta_2l,1),lm)
1347 
1348  varname='atij'
1349  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1350  spval,varname,atij(ista_2l,jsta_2l,1),lm)
1351 
1352  varname='atol1j'
1353  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1354  spval,varname,atol1j(ista_2l,jsta_2l,1),lm)
1355 
1356  varname='atol2j'
1357  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1358  spval,varname,atol2j(ista_2l,jsta_2l,1),lm)
1359 
1360  varname='atol3j'
1361  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1362  spval,varname,atol3j(ista_2l,jsta_2l,1),lm)
1363 
1364  varname='atrp1j'
1365  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1366  spval,varname,atrp1j(ista_2l,jsta_2l,1),lm)
1367 
1368  varname='atrp2j'
1369  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1370  spval,varname,atrp2j(ista_2l,jsta_2l,1),lm)
1371 
1372  varname='axyl1j'
1373  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1374  spval,varname,axyl1j(ista_2l,jsta_2l,1),lm)
1375 
1376  varname='axyl2j'
1377  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1378  spval,varname,axyl2j(ista_2l,jsta_2l,1),lm)
1379 
1380  varname='axyl3j'
1381  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1382  spval,varname,axyl3j(ista_2l,jsta_2l,1),lm)
1383 
1384  varname='pm25ac'
1385  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1386  spval,varname,pm25ac(ista_2l,jsta_2l,1),lm)
1387 
1388  varname='pm25at'
1389  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1390  spval,varname,pm25at(ista_2l,jsta_2l,1),lm)
1391 
1392  varname='pm25co'
1393  call read_netcdf_3d_para(ncid3d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1394  spval,varname,pm25co(ista_2l,jsta_2l,1),lm)
1395 
1396 !=========================
1397 ! PM2.5 SPECIES
1398 !=========================
1399 
1400  ! do l=1,lm
1401  ! do j=jsta,jend
1402  ! do i=ista,iend
1403  ! pm25hp(i,j,l) = ( ah3opi(i,j,l)*pm25at(i,j,l) &
1404  ! + ah3opj(i,j,l)*pm25ac(i,j,l) &
1405  ! + ah3opk(i,j,l)*pm25co(i,j,l) ) / 19.0
1406 
1407  ! pm25cl(i,j,l) = acli(i,j,l)*pm25at(i,j,l) &
1408  ! + aclj(i,j,l)*pm25ac(i,j,l) &
1409  ! + aclk(i,j,l)*pm25co(i,j,l)
1410 
1411  ! pm25ec(i,j,l) = aeci(i,j,l)*pm25at(i,j,l) &
1412  ! + aecj(i,j,l)*pm25ac(i,j,l)
1413  ! enddo
1414  ! enddo
1415  ! enddo
1416 
1417 
1418  ! do l=1,lm
1419  ! do j=jsta,jend
1420  ! do i=ista,iend
1421 
1422  ! anak(i,j,l) = 0.8373 * aseacat(i,j,l) &
1423  ! + 0.0626 * asoil(i,j,l) &
1424  ! + 0.0023 * acors(i,j,l)
1425 
1426  ! pm25na(i,j,l) = anai(i,j,l)*pm25at(i,j,l) &
1427  ! + anaj(i,j,l)*pm25ac(i,j,l) &
1428  ! + anak(i,j,l)*pm25co(i,j,l)
1429  ! enddo
1430  ! enddo
1431  ! enddo
1432 
1433  do l=1,lm
1434  do j=jsta,jend
1435  do i=ista,iend
1436 
1437  apomi(i,j,l) = alvpo1i(i,j,l) &
1438  +asvpo1i(i,j,l) + asvpo2i(i,j,l)
1439 
1440  apomj(i,j,l) = alvpo1j(i,j,l) &
1441  +asvpo1j(i,j,l) + asvpo2j(i,j,l) + asvpo3j(i,j,l) &
1442  +aivpo1j(i,j,l)
1443 
1444  asomi(i,j,l) = alvoo1i(i,j,l) + alvoo2i(i,j,l) &
1445  +asvoo1i(i,j,l) + asvoo2i(i,j,l)
1446 
1447  asomj(i,j,l) = axyl1j(i,j,l) + axyl2j(i,j,l) + axyl3j(i,j,l) &
1448  +atol1j(i,j,l) + atol2j(i,j,l) + atol3j(i,j,l) &
1449  +abnz1j(i,j,l) + abnz2j(i,j,l) + abnz3j(i,j,l) &
1450  +aiso1j(i,j,l) + aiso2j(i,j,l) + aiso3j(i,j,l) &
1451  +atrp1j(i,j,l) + atrp2j(i,j,l) + asqtj(i,j,l) &
1452  +aalk1j(i,j,l) + aalk2j(i,j,l) &
1453  +apah1j(i,j,l) + apah2j(i,j,l) + apah3j(i,j,l) &
1454  +aorgcj(i,j,l) + aolgbj(i,j,l) + aolgaj(i,j,l) &
1455  +alvoo1j(i,j,l) + alvoo2j(i,j,l) &
1456  +asvoo1j(i,j,l) + asvoo2j(i,j,l) + asvoo3j(i,j,l) &
1457  +apcsoj(i,j,l)
1458 
1459  aomi(i,j,l) = apomi(i,j,l) + asomi(i,j,l)
1460  aomj(i,j,l) = apomj(i,j,l) + asomj(i,j,l)
1461 
1462  atoti(i,j,l) = aso4i(i,j,l) + ano3i(i,j,l) + anh4i(i,j,l) &
1463  + anai(i,j,l) + acli(i,j,l) + aeci(i,j,l) &
1464  + aomi(i,j,l) +aothri(i,j,l)
1465 
1466  atotj(i,j,l) = aso4j(i,j,l) + ano3j(i,j,l) + anh4j(i,j,l) &
1467  + anaj(i,j,l) + aclj(i,j,l) + aecj(i,j,l) &
1468  + aomj(i,j,l) +aothrj(i,j,l) &
1469  + afej(i,j,l) + asij(i,j,l) + atij(i,j,l) &
1470  + acaj(i,j,l) + amgj(i,j,l) + amnj(i,j,l) &
1471  + aalj(i,j,l) + akj(i,j,l)
1472 
1473  atotk(i,j,l) = asoil(i,j,l) + acors(i,j,l) + aseacat(i,j,l)&
1474  + aclk(i,j,l) &
1475  +aso4k(i,j,l) + ano3k(i,j,l) + anh4k(i,j,l)
1476 
1477  pmtf(i,j,l) = atoti(i,j,l)*pm25at(i,j,l) &
1478  + atotj(i,j,l)*pm25ac(i,j,l) &
1479  + atotk(i,j,l)*pm25co(i,j,l)
1480  enddo
1481  enddo
1482  enddo
1483 
1484  deallocate (aacd, aalj, aalk1j, aalk2j, abnz1j, abnz2j, abnz3j)
1485  deallocate (acaj, acet, acli, aclj, aclk)
1486  deallocate (acors, acro_primary, acrolein)
1487 
1488  deallocate (aeci, aecj, afej, aglyj, ah2oi, ah2oj, ah2ok)
1489  deallocate (ah3opi, ah3opj, ah3opk, aiso1j, aiso2j, aiso3j)
1490 
1491  deallocate (aivpo1j, akj, ald2, ald2_primary, aldx)
1492  deallocate (alvoo1i, alvoo1j, alvoo2i, alvoo2j, alvpo1i, alvpo1j)
1493 
1494  deallocate (amgj, amnj, anai, anaj, anak)
1495  deallocate (anh4i, anh4j, anh4k, ano3i, ano3j, ano3k)
1496 
1497  deallocate (aolgaj, aolgbj, aomi, aomj)
1498  deallocate (aorgcj, aothri, aothrj, apah1j, apah2j, apah3j)
1499 
1500  deallocate (apcsoj, apomi, apomj, aseacat, asij)
1501  deallocate (aso4i, aso4j, aso4k, asoil, asomi, asomj, asqtj)
1502 
1503  deallocate (asvoo1i, asvoo1j, asvoo2i, asvoo2j, asvoo3j)
1504  deallocate (asvpo1i, asvpo1j, asvpo2i, asvpo2j, asvpo3j)
1505 
1506  deallocate (atij, atol1j, atol2j, atol3j, atrp1j, atrp2j)
1507  deallocate (atoti, atotj, atotk, axyl1j, axyl2j, axyl3j)
1508 
1509  deallocate (pm25ac, pm25at, pm25co)
1510 
1511  endif ! -- aqfcmaq_on
1512 !============================
1513 
1514 ! read for regional FV3
1515  if (modelname == 'FV3R') then
1516 ! max hourly updraft velocity
1517  varname='upvvelmax'
1518  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1519  spval,varname,w_up_max(ista_2l,jsta_2l))
1520  if(debugprint)print*,'sample ',varname,' = ',w_up_max(isa,jsa)
1521 ! max hourly downdraft velocity
1522  varname='dnvvelmax'
1523  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1524  spval,varname,w_dn_max(ista_2l,jsta_2l))
1525  if(debugprint)print*,'sample ',varname,' = ',w_dn_max(isa,jsa)
1526 ! max hourly updraft helicity
1527  varname='uhmax25'
1528  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1529  spval,varname,up_heli_max(ista_2l,jsta_2l))
1530  if(debugprint)print*,'sample ',varname,' = ',up_heli_max(isa,jsa)
1531 ! min hourly updraft helicity
1532  varname='uhmin25'
1533  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1534  spval,varname,up_heli_min(ista_2l,jsta_2l))
1535  if(debugprint)print*,'sample ',varname,' = ',up_heli_min(isa,jsa)
1536 ! max hourly 0-3km updraft helicity
1537  varname='uhmax03'
1538  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1539  spval,varname,up_heli_max03(ista_2l,jsta_2l))
1540  if(debugprint)print*,'sample ',varname,' = ',up_heli_max03(isa,jsa)
1541 ! min hourly 0-3km updraft helicity
1542  varname='uhmin03'
1543  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1544  spval,varname,up_heli_min03(ista_2l,jsta_2l))
1545  if(debugprint)print*,'sample ',varname,' = ',up_heli_min03(isa,jsa)
1546 
1547 ! max 0-1km relative vorticity max
1548  varname='maxvort01'
1549  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1550  spval,varname,rel_vort_max01(ista_2l,jsta_2l))
1551  if(debugprint)print*,'sample ',varname,' = ',rel_vort_max01(isa,jsa)
1552 ! max 0-2km relative vorticity max
1553  varname='maxvort02'
1554  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1555  spval,varname,rel_vort_max(ista_2l,jsta_2l))
1556  if(debugprint)print*,'sample ',varname,' =',rel_vort_max(isa,jsa)
1557 ! max hybrid lev 1 relative vorticity max
1558  varname='maxvorthy1'
1559  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1560  spval,varname,rel_vort_maxhy1(ista_2l,jsta_2l))
1561  if(debugprint)print*,'sample ',varname,' =',rel_vort_maxhy1(isa,jsa)
1562  endif
1563 
1564 ! surface pressure
1565  varname='pressfc'
1566  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1567  spval,varname,pint(ista_2l,jsta_2l,lp1))
1568  do j=jsta,jend
1569  do i=ista,iend
1570 ! if(pint(i,j,lp1)>1.0E6 .or. pint(ista_2l,jsta_2l,lp1)<50000.) &
1571 ! print*,'bad psfc ',i,j,pint(i,j,lp1)
1572  end do
1573  end do
1574  if(debugprint)print*,'sample ',varname,' =',pint(isa,jsa,lp1)
1575 
1576  pt = ak5(1)
1577 
1578  do j=jsta,jend
1579  do i=ista,iend
1580  pint(i,j,1)= pt
1581  end do
1582  end do
1583 
1584  do l=2,lp1
1585  do j=jsta,jend
1586  do i=ista,iend
1587  if (dpres(i,j,l-1)<spval .and. pint(i,j,l-1)<spval) then
1588  pint(i,j,l)= pint(i,j,l-1) + dpres(i,j,l-1)
1589  else
1590  pint(i,j,l)=spval
1591  endif
1592  enddo
1593  enddo
1594 ! if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
1595 ! ,pint(ii,jj,l),pmid(ii,jj,l)
1596  end do
1597 
1598 !compute pmid from averaged two layer pint
1599  do l=lm,1,-1
1600  do j=jsta,jend
1601  do i=ista,iend
1602  if (pint(i,j,l)<spval .and. pint(i,j,l+1)<spval) then
1603  pmid(i,j,l)=0.5*(pint(i,j,l)+pint(i,j,l+1))
1604  else
1605  pmid(i,j,l)=spval
1606  endif
1607  enddo
1608  enddo
1609  enddo
1610 
1611 ! surface height from FV3
1612 ! dong set missing value for zint
1613 ! zint=spval
1614  varname='hgtsfc'
1615  call read_netcdf_2d_para(ncid3d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1616  spval,varname,zint(ista_2l,jsta_2l,lp1))
1617  if(debugprint)print*,'sample ',varname,' =',zint(isa,jsa,lp1)
1618  do j=jsta,jend
1619  do i=ista,iend
1620  if (zint(i,j,lp1) /= spval) then
1621  fis(i,j) = zint(i,j,lp1) * grav
1622  else
1623  fis(i,j) = spval
1624  endif
1625  enddo
1626  enddo
1627 
1628  do l=lm,1,-1
1629  do j=jsta,jend
1630  do i=ista,iend
1631  if(zint(i,j,l+1)/=spval .and. buf3d(i,j,l)/=spval)then
1632 !make sure delz is positive
1633  zint(i,j,l)=zint(i,j,l+1)+abs(buf3d(i,j,l))
1634 ! if(zint(i,j,l)>1.0E6)print*,'bad H ',i,j,l,zint(i,j,l)
1635  else
1636  zint(i,j,l)=spval
1637  end if
1638  end do
1639  end do
1640  if(debugprint)print*,'sample zint= ',isa,jsa,l,zint(isa,jsa,l)
1641  end do
1642 
1643  do l=lp1,1,-1
1644  do j=jsta,jend
1645  do i=ista,iend
1646  alpint(i,j,l)=log(pint(i,j,l))
1647  end do
1648  end do
1649  end do
1650 
1651  do l=lm,1,-1
1652  do j=jsta,jend
1653  do i=ista,iend
1654  if(zint(i,j,l+1)/=spval .and. zint(i,j,l)/=spval &
1655  .and. pmid(i,j,l)/=spval)then
1656  zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* &
1657  (log(pmid(i,j,l))-alpint(i,j,l+1))/ &
1658  (alpint(i,j,l)-alpint(i,j,l+1))
1659  if(zmid(i,j,l)>1.0e6)print*,'bad Hmid ',i,j,l,zmid(i,j,l)
1660  else
1661  zmid(i,j,l)=spval
1662  endif
1663  end do
1664  end do
1665  end do
1666 
1667 
1668 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1669 
1670 !
1671 
1672 ! done with 3d file, close it for now
1673  status=nf90_close(ncid3d)
1674  deallocate(recname)
1675 
1676 ! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD
1677  varname='IVEGSRC'
1678  status=nf90_get_att(ncid2d,nf90_global,'IVEGSRC',ivegsrc)
1679  if (status /= 0) then
1680  print*,varname,' not found-Assigned 1 for IGBP as default'
1681  ivegsrc=1
1682  end if
1683  if (me == 0) print*,'IVEGSRC= ',ivegsrc
1684 
1685 ! set novegtype based on vegetation classification
1686  if(ivegsrc==2)then
1687  novegtype=13
1688  else if(ivegsrc==1)then
1689  novegtype=20
1690  else if(ivegsrc==0)then
1691  novegtype=24
1692  end if
1693  if (me == 0) print*,'novegtype= ',novegtype
1694 
1695  status=nf90_get_att(ncid2d,nf90_global,'fhzero',fhzero)
1696  if (status /= 0) then
1697  print*,'fhzero not found-Assigned 3 hours as default'
1698  fhzero=3
1699  end if
1700  if (me == 0) print*,'fhzero= ',fhzero
1701 !
1702  status=nf90_get_att(ncid2d,nf90_global,'dtp',dtp)
1703  if (status /= 0) then
1704  print*,'dtp not found-Assigned 90s as default'
1705  dtp=90
1706  end if
1707  if (me == 0) print*,'dtp= ',dtp
1708 ! Initializes constants for Ferrier microphysics
1709  if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then
1710  CALL microinit(imp_physics)
1711  end if
1712 
1713  tprec = float(fhzero)
1714  if(ifhr>240)tprec=12.
1715  tclod = tprec
1716  trdlw = tprec
1717  trdsw = tprec
1718  tsrfc = tprec
1719  tmaxmin = tprec
1720  td3d = tprec
1721  print*,'tprec = ',tprec
1722 
1723 
1724  varname='refl_10cm'
1725  call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1726  spval,varname,ref_10cm(ista_2l,jsta_2l,1),lm)
1727 ! do l=1,lm
1728 ! if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' &
1729 ! ,REF_10CM(isa,jsa,l),isa,jsa,l
1730 ! enddo
1731 
1732 ! turbulence kinetic energy (QKE = 2*TKE)
1733  varname='qke'
1734  call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1735  spval,varname,q2(ista_2l,jsta_2l,1),lm)
1736  do l=1,lm
1737  do j=jsta,jend
1738  do i=ista,iend
1739  q2(i,j,l)=q2(i,j,l)/2.0
1740  enddo
1741  enddo
1742  enddo
1743 
1744 ! ice-friendly aerosol number concentration
1745  varname='nifa'
1746  call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1747  spval,varname,qqnifa(ista_2l,jsta_2l,1),lm)
1748 
1749 ! water-friendly aerosol number concentration
1750  varname='nwfa'
1751  call read_netcdf_3d_para(ncid2d,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1752  spval,varname,qqnwfa(ista_2l,jsta_2l,1),lm)
1753 
1754  varname='land'
1755  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1756  spval,varname,sm)
1757  if(debugprint)print*,'sample ',varname,' =',sm((ista+iend)/2,(jsta+jend)/2)
1758 
1759 !$omp parallel do private(i,j)
1760  do j=jsta,jend
1761  do i=ista,iend
1762  if (sm(i,j) /= spval) sm(i,j) = 1.0 - sm(i,j)
1763  enddo
1764  enddo
1765 
1766 ! sea ice mask
1767 
1768  varname = 'icec'
1769  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1770  spval,varname,sice)
1771  if(debugprint)print*,'sample ',varname,' = ',sice(isa,jsa)
1772 
1773 ! where(sice /=spval .and. sice >=1.0)sm=0.0 !sea ice has sea
1774 ! mask=0
1775 ! GFS flux files have land points with non-zero sea ice, per Iredell,
1776 ! these
1777 ! points have sea ice changed to zero, i.e., trust land mask more than
1778 ! sea ice
1779 ! where(sm/=spval .and. sm==0.0)sice=0.0 !specify sea ice=0 at land
1780 
1781 !$omp parallel do private(i,j)
1782  do j=jsta,jend
1783  do i=ista,iend
1784  if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0
1785  enddo
1786  enddo
1787 
1788 
1789 ! PBL height using nemsio
1790  varname = 'hpbl'
1791  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1792  spval,varname,pblh)
1793  if(debugprint)print*,'sample ',varname,' = ',pblh(isa,jsa)
1794 
1795 ! frictional velocity using nemsio
1796  varname='fricv'
1797  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1798  spval,varname,ustar)
1799 ! if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa)
1800 
1801 ! roughness length using getgb
1802  varname='sfcr'
1803  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1804  spval,varname,z0)
1805 ! if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa)
1806 
1807 ! sfc exchange coeff
1808  varname='sfexc'
1809  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1810  spval,varname,sfcexc)
1811 
1812 ! aerodynamic conductance
1813  varname='acond'
1814  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1815  spval,varname,acond)
1816  if(debugprint)print*,'sample ',varname,' = ',acond(isa,jsa)
1817 ! mid day avg albedo
1818  varname='albdo_ave'
1819  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1820  spval,varname,avgalbedo)
1821 !$omp parallel do private(i,j)
1822  do j=jsta,jend
1823  do i=ista,iend
1824  if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01
1825  enddo
1826  enddo
1827  if(debugprint)print*,'sample ',varname,' = ',avgalbedo(isa,jsa)
1828 
1829 ! surface potential T using getgb
1830  varname='tmpsfc'
1831  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1832  spval,varname,ths)
1833 
1834 ! where(ths/=spval)ths=ths*(p1000/pint(:,:,lp1))**CAPA ! convert to THS
1835 
1836 !$omp parallel do private(i,j)
1837  do j=jsta,jend
1838  do i=ista,iend
1839  if (ths(i,j) /= spval) then
1840 ! write(0,*)' i=',i,' j=',j,' ths=',ths(i,j),' pint=',pint(i,j,lp1)
1841  ths(i,j) = ths(i,j) * (p1000/pint(i,j,lp1))**capa
1842  endif
1843  qs(i,j) = spval ! GFS does not have surface specific humidity
1844  twbs(i,j) = spval ! GFS does not have inst sensible heat flux
1845  qwbs(i,j) = spval ! GFS does not have inst latent heat flux
1846 !assign sst
1847  if (sm(i,j) /= 0.0 .and. ths(i,j) < spval ) then
1848  if (sice(i,j) >= 0.15) then
1849  sst(i,j) = 271.4
1850  else
1851  sst(i,j) = ths(i,j) * (pint(i,j,lp1)/p1000)**capa
1852  endif
1853  else
1854  sst(i,j) = spval
1855  endif
1856  enddo
1857  enddo
1858  if(debugprint)print*,'sample ',varname,' = ',ths(isa,jsa)
1859 
1860 ! foundation temperature
1861  varname='tref'
1862  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1863  spval,varname,fdnsst)
1864  if(debugprint)print*,'sample ',varname,' = ',fdnsst(isa,jsa)
1865 
1866 
1867 ! GFS does not have time step and physics time step, make up ones since they
1868 ! are not really used anyway
1869 ! NPHS=1.
1870 ! DT=90.
1871 ! DTQ2 = DT * NPHS !MEB need to get physics DT
1872  dtq2 = dtp !MEB need to get physics DT
1873  nphs=1
1874  dt = dtq2/nphs !MEB need to get DT
1875  tsph = 3600./dt
1876 
1877 ! convective precip in m per physics time step using getgb
1878 ! read 3 hour bucket
1879  varname='cpratb_ave'
1880  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1881  spval,varname,avgcprate)
1882 ! where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m
1883 !$omp parallel do private(i,j)
1884  do j=jsta,jend
1885  do i=ista,iend
1886  if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001)
1887  enddo
1888  enddo
1889 ! if(debugprint)print*,'sample ',VarName,' = ',avgcprate(isa,jsa)
1890 
1891 ! print*,'maxval CPRATE: ', maxval(CPRATE)
1892 
1893 ! read continuous bucket
1894  varname='cprat_ave'
1895  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1896  spval,varname,avgcprate_cont)
1897 !$omp parallel do private(i,j)
1898  do j=jsta,jend
1899  do i=ista,iend
1900  if (avgcprate_cont(i,j) /= spval) avgcprate_cont(i,j) = &
1901  avgcprate_cont(i,j) * (dtq2*0.001)
1902  enddo
1903  enddo
1904 ! if(debugprint)print*,'sample ',VarName,' = ',avgcprate_cont(isa,jsa)
1905 
1906 ! print*,'maxval CPRATE: ', maxval(CPRATE)
1907 
1908 ! precip rate in m per physics time step using getgb
1909  varname='prateb_ave'
1910  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1911  spval,varname,avgprec)
1912 !$omp parallel do private(i,j)
1913  do j=jsta,jend
1914  do i=ista,iend
1915  if(avgprec(i,j) /= spval)avgprec(i,j)=avgprec(i,j)*(dtq2*0.001)
1916  enddo
1917  enddo
1918 
1919  if(debugprint)print*,'sample ',varname,' = ',avgprec(isa,jsa)
1920 
1921 ! prec = avgprec !set avg cprate to inst one to derive other fields
1922 
1923  varname='prate_ave'
1924  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1925  spval,varname,avgprec_cont)
1926 ! where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
1927 !$omp parallel do private(i,j)
1928  do j=jsta,jend
1929  do i=ista,iend
1930  if (avgprec_cont(i,j) /=spval)avgprec_cont(i,j)=avgprec_cont(i,j) &
1931  * (dtq2*0.001)
1932  enddo
1933  enddo
1934 
1935  if(debugprint)print*,'sample ',varname,' = ',avgprec_cont(isa,jsa)
1936 ! precip rate in m per physics time step
1937  varname='tprcp'
1938  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1939  spval,varname,prec)
1940 !$omp parallel do private(i,j)
1941  do j=jsta,jend
1942  do i=ista,iend
1943  if (prec(i,j) /= spval) prec(i,j)=prec(i,j)* (dtq2*0.001) &
1944  * 1000. / dtp
1945  enddo
1946  enddo
1947 
1948 ! convective precip rate in m per physics time step
1949  varname='cnvprcp'
1950  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1951  spval,varname,cprate)
1952 !set cprate as 0.
1953  do j=jsta,jend
1954  do i=ista,iend
1955  if (cprate(i,j) /= spval) then
1956  cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) &
1957  * 1000. / dtp
1958  else
1959  cprate(i,j) = 0.
1960  endif
1961  enddo
1962  enddo
1963 
1964 ! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f
1965 
1966 ! max hourly surface precipitation rate
1967  varname='pratemax'
1968  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1969  spval,varname,prate_max)
1970  if(debugprint)print*,'sample ',varname,' = ',prate_max(isa,jsa)
1971 ! max hourly 1-km agl reflectivity
1972  varname='refdmax'
1973  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1974  spval,varname,refd_max)
1975  if(debugprint)print*,'sample ',varname,' = ',refd_max(isa,jsa)
1976 ! max hourly -10C reflectivity
1977  varname='refdmax263k'
1978  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1979  spval,varname,refdm10c_max)
1980  if(debugprint)print*,'sample ',varname,' = ',refdm10c_max(isa,jsa)
1981 
1982 ! max hourly u comp of 10m agl wind
1983  varname='u10max'
1984  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1985  spval,varname,u10max)
1986  if(debugprint)print*,'sample ',varname,' = ',u10max(isa,jsa)
1987 ! max hourly v comp of 10m agl wind
1988  varname='v10max'
1989  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1990  spval,varname,v10max)
1991  if(debugprint)print*,'sample ',varname,' = ',v10max(isa,jsa)
1992 ! max hourly 10m agl wind speed
1993  varname='spd10max'
1994  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
1995  spval,varname,wspd10max)
1996  if(debugprint)print*,'sample ',varname,' = ',wspd10max(isa,jsa)
1997 
1998 ! inst snow water eqivalent using nemsio
1999  varname='weasd'
2000  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2001  spval,varname,sno)
2002 ! mask water areas
2003 !$omp parallel do private(i,j)
2004  do j=jsta,jend
2005  do i=ista,iend
2006  if (sm(i,j) == 1.0 .and. sice(i,j)==0.) sno(i,j) = spval
2007  enddo
2008  enddo
2009  if(debugprint)print*,'sample ',varname,' = ',sno(isa,jsa)
2010 
2011 ! ave snow cover
2012  varname='snowc_ave'
2013  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2014  spval,varname,snoavg)
2015 ! snow cover is multipled by 100 in SURFCE before writing it out
2016  do j=jsta,jend
2017  do i=ista,iend
2018  if (sm(i,j)==1.0 .and. sice(i,j)==0.) snoavg(i,j)=spval
2019  if(snoavg(i,j)/=spval)snoavg(i,j)=snoavg(i,j)/100.
2020  end do
2021  end do
2022 
2023 ! snow depth in mm using nemsio
2024  varname='snod'
2025  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2026  spval,varname,si)
2027 !$omp parallel do private(i,j)
2028  do j=jsta,jend
2029  do i=ista,iend
2030  if (sm(i,j)==1.0 .and. sice(i,j)==0.) si(i,j)=spval
2031  if (si(i,j) /= spval) si(i,j) = si(i,j) * 1000.0
2032  cldefi(i,j) = spval ! GFS does not have convective cloud efficiency
2033  lspa(i,j) = spval ! GFS does not have similated precip
2034  th10(i,j) = spval ! GFS does not have 10 m theta
2035  th10(i,j) = spval ! GFS does not have 10 m theta
2036  q10(i,j) = spval ! GFS does not have 10 m humidity
2037  albase(i,j) = spval ! GFS does not have snow free albedo
2038  enddo
2039  enddo
2040  if(debugprint)print*,'sample ',varname,' = ',si(isa,jsa)
2041 
2042 ! 2m T using nemsio
2043  varname='tmp2m'
2044  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2045  spval,varname,tshltr)
2046  if(debugprint)print*,'sample ',varname,' = ',tshltr(isa,jsa)
2047 
2048 ! GFS does not have 2m pres, estimate it, also convert t to theta
2049  Do j=jsta,jend
2050  Do i=ista,iend
2051  pshltr(i,j)=pint(i,j,lm+1)*exp(-0.068283/tshltr(i,j))
2052  tshltr(i,j)= tshltr(i,j)*(p1000/pshltr(i,j))**capa ! convert to theta
2053 ! if (j == jm/2 .and. mod(i,50) == 0)
2054 ! + print*,'sample 2m T and P after scatter= '
2055 ! + ,i,j,tshltr(i,j),pshltr(i,j)
2056  end do
2057  end do
2058 
2059 ! 2m specific humidity using nemsio
2060  varname='spfh2m'
2061  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2062  spval,varname,qshltr)
2063  if(debugprint)print*,'sample ',varname,' = ',qshltr(isa,jsa)
2064 
2065 ! time averaged column cloud fractionusing nemsio
2066  varname='tcdc_aveclm'
2067  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2068  spval,varname,avgtcdc)
2069 ! where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction
2070 !$omp parallel do private(i,j)
2071  do j=jsta,jend
2072  do i=ista,iend
2073  if (avgtcdc(i,j) /= spval) avgtcdc(i,j) = avgtcdc(i,j) * 0.01
2074  enddo
2075  enddo
2076  if(debugprint)print*,'sample ',varname,' = ',avgtcdc(isa,jsa)
2077 
2078 ! GFS probably does not use zenith angle
2079 !$omp parallel do private(i,j)
2080  do j=jsta_2l,jend_2u
2081  do i=ista_2l,iend_2u
2082  czen(i,j) = spval
2083  czmean(i,j) = spval
2084  enddo
2085  enddo
2086 
2087 ! maximum snow albedo in fraction using nemsio
2088  varname='snoalb'
2089  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2090  spval,varname,mxsnal)
2091 
2092 ! land fraction
2093  varname='lfrac'
2094  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2095  spval,varname,landfrac)
2096 
2097 ! GFS probably does not use sigt4, set it to sig*t^4
2098 !$omp parallel do private(i,j,tlmh)
2099  Do j=jsta,jend
2100  Do i=ista,iend
2101  tlmh = t(i,j,lm) * t(i,j,lm)
2102  sigt4(i,j) = 5.67e-8 * tlmh * tlmh
2103  End do
2104  End do
2105 
2106 ! TG is not used, skip it for now
2107 
2108 ! GFS does not have inst cloud fraction for high, middle, and low cloud
2109 !$omp parallel do private(i,j)
2110  do j=jsta_2l,jend_2u
2111  do i=ista_2l,iend_2u
2112  cfrach(i,j) = spval
2113  cfracl(i,j) = spval
2114  cfracm(i,j) = spval
2115  enddo
2116  enddo
2117 
2118 ! ave high cloud fraction using nemsio
2119  varname='tcdc_avehcl'
2120  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2121  spval,varname,avgcfrach)
2122 ! where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction
2123 !$omp parallel do private(i,j)
2124  do j=jsta,jend
2125  do i=ista,iend
2126  if (avgcfrach(i,j) /= spval) avgcfrach(i,j) = avgcfrach(i,j) * 0.01
2127  enddo
2128  enddo
2129  if(debugprint)print*,'sample ',varname,' = ',avgcfrach(isa,jsa)
2130 
2131 ! ave low cloud fraction using nemsio
2132  varname='tcdc_avelcl'
2133  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2134  spval,varname,avgcfracl)
2135 ! where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction
2136 !$omp parallel do private(i,j)
2137  do j=jsta,jend
2138  do i=ista,iend
2139  if (avgcfracl(i,j) /= spval) avgcfracl(i,j) = avgcfracl(i,j) * 0.01
2140  enddo
2141  enddo
2142  if(debugprint)print*,'sample ',varname,' = ',avgcfracl(isa,jsa)
2143 
2144 ! ave middle cloud fraction using nemsio
2145  varname='tcdc_avemcl'
2146  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2147  spval,varname,avgcfracm)
2148 ! where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction
2149 !$omp parallel do private(i,j)
2150  do j=jsta,jend
2151  do i=ista,iend
2152  if (avgcfracm(i,j) /= spval) avgcfracm(i,j) = avgcfracm(i,j) * 0.01
2153  enddo
2154  enddo
2155  if(debugprint)print*,'sample ',varname,' = ',avgcfracm(isa,jsa)
2156 
2157 ! inst convective cloud fraction using nemsio
2158  varname='tcdccnvcl'
2159  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2160  spval,varname,cnvcfr)
2161 ! where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction
2162 !$omp parallel do private(i,j)
2163  do j=jsta,jend
2164  do i=ista,iend
2165  if (cnvcfr(i,j) /= spval) cnvcfr(i,j)= cnvcfr(i,j) * 0.01
2166  enddo
2167  enddo
2168 ! if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(isa,jsa)
2169 
2170 ! slope type using nemsio
2171  varname='sltyp'
2172  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2173  spval,varname,buf)
2174 !$omp parallel do private(i,j)
2175  do j = jsta_2l, jend_2u
2176  do i=ista,iend
2177  if (buf(i,j) < spval) then
2178  islope(i,j) = nint(buf(i,j))
2179  else
2180  islope(i,j) = 0
2181  endif
2182  enddo
2183  enddo
2184 ! if(debugprint)print*,'sample ',VarName,' = ',islope(isa,jsa)
2185 
2186 ! plant canopy sfc wtr in m
2187  varname='cnwat'
2188  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2189  spval,varname,cmc)
2190 !$omp parallel do private(i,j)
2191  do j=jsta,jend
2192  do i=ista,iend
2193  if (cmc(i,j) /= spval) cmc(i,j) = cmc(i,j) * 0.001
2194  if (sm(i,j) /= 0.0) cmc(i,j) = spval
2195  enddo
2196  enddo
2197 ! if(debugprint)print*,'sample ',VarName,' = ',cmc(isa,jsa)
2198 
2199 !$omp parallel do private(i,j)
2200  do j=jsta_2l,jend_2u
2201  do i=ista_2l,iend_2u
2202  grnflx(i,j) = spval ! GFS does not have inst ground heat flux
2203  enddo
2204  enddo
2205 
2206 ! frozen precip fraction using nemsio
2207  varname='cpofp'
2208  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2209  spval,varname,sr)
2210 !$omp parallel do private(i,j)
2211  do j=jsta,jend
2212  do i=ista,iend
2213  if(sr(i,j) /= spval) then
2214 !set range within (0,1)
2215  sr(i,j)=min(1.,max(0.,sr(i,j)))
2216  endif
2217  enddo
2218  enddo
2219 
2220 ! sea ice skin temperature
2221  varname='tisfc'
2222  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2223  spval,varname,ti)
2224 !$omp parallel do private(i,j)
2225  do j=jsta,jend
2226  do i=ista,iend
2227  if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval
2228  enddo
2229  enddo
2230 
2231 ! vegetation fraction in fraction. using nemsio
2232  varname='veg'
2233  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2234  spval,varname,vegfrc)
2235 !$omp parallel do private(i,j)
2236  do j=jsta,jend
2237  do i=ista,iend
2238  if (vegfrc(i,j) /= spval) then
2239  vegfrc(i,j) = vegfrc(i,j) * 0.01
2240  else
2241  vegfrc(i,j) = 0.0
2242  endif
2243  enddo
2244  enddo
2245 ! mask water areas
2246 !$omp parallel do private(i,j)
2247  do j=jsta,jend
2248  do i=ista,iend
2249  if (sm(i,j) /= 0.0) vegfrc(i,j) = spval
2250  enddo
2251  enddo
2252 ! if(debugprint)print*,'sample ',VarName,' = ',vegfrc(isa,jsa)
2253 
2254 ! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam
2255 
2256  sldpth(1) = 0.10
2257  sldpth(2) = 0.3
2258  sldpth(3) = 0.6
2259  sldpth(4) = 1.0
2260 
2261 ! Eric James, 1 Oct 2021: Because FV3 does not have 1d var "zs", used to
2262 ! assign soil depths for RUC LSM, hard wire 9 soil depths here
2263 ! so they aren't missing.
2264 
2265  IF (nsoil==9) THEN
2266  sllevel(1) = 0.0
2267  sllevel(2) = 0.01
2268  sllevel(3) = 0.04
2269  sllevel(4) = 0.1
2270  sllevel(5) = 0.3
2271  sllevel(6) = 0.6
2272  sllevel(7) = 1.0
2273  sllevel(8) = 1.6
2274  sllevel(9) = 3.0
2275  END IF
2276 
2277 ! liquid volumetric soil mpisture in fraction using nemsio
2278  varname='soill1'
2279  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2280  spval,varname,sh2o(ista_2l,jsta_2l,1))
2281 ! mask water areas
2282 !$omp parallel do private(i,j)
2283  do j=jsta,jend
2284  do i=ista,iend
2285  if (sm(i,j) /= 0.0) sh2o(i,j,1) = spval
2286  enddo
2287  enddo
2288  if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,1)
2289 
2290  varname='soill2'
2291  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2292  spval,varname,sh2o(ista_2l,jsta_2l,2))
2293 ! mask water areas
2294 !$omp parallel do private(i,j)
2295  do j=jsta,jend
2296  do i=ista,iend
2297  if (sm(i,j) /= 0.0) sh2o(i,j,2) = spval
2298  enddo
2299  enddo
2300  if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,2)
2301 
2302  varname='soill3'
2303  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2304  spval,varname,sh2o(ista_2l,jsta_2l,3))
2305 ! mask water areas
2306 !$omp parallel do private(i,j)
2307  do j=jsta,jend
2308  do i=ista,iend
2309  if (sm(i,j) /= 0.0) sh2o(i,j,3) = spval
2310  enddo
2311  enddo
2312  if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,3)
2313 
2314  varname='soill4'
2315  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2316  spval,varname,sh2o(ista_2l,jsta_2l,4))
2317 ! mask water areas
2318 !$omp parallel do private(i,j)
2319  do j=jsta,jend
2320  do i=ista,iend
2321  if (sm(i,j) /= 0.0) sh2o(i,j,4) = spval
2322  enddo
2323  enddo
2324  if(debugprint)print*,'sample l',varname,' = ',1,sh2o(isa,jsa,4)
2325 
2326 ! volumetric soil moisture using nemsio
2327  varname='soilw1'
2328  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2329  spval,varname,smc(ista_2l,jsta_2l,1))
2330 ! mask water areas
2331 !$omp parallel do private(i,j)
2332  do j=jsta,jend
2333  do i=ista,iend
2334  if (sm(i,j) /= 0.0) smc(i,j,1) = spval
2335  enddo
2336  enddo
2337  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,1)
2338 
2339  varname='soilw2'
2340  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2341  spval,varname,smc(ista_2l,jsta_2l,2))
2342 ! mask water areas
2343 !$omp parallel do private(i,j)
2344  do j=jsta,jend
2345  do i=ista,iend
2346  if (sm(i,j) /= 0.0) smc(i,j,2) = spval
2347  enddo
2348  enddo
2349  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,2)
2350 
2351  varname='soilw3'
2352  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2353  spval,varname,smc(ista_2l,jsta_2l,3))
2354 ! mask water areas
2355 !$omp parallel do private(i,j)
2356  do j=jsta,jend
2357  do i=ista,iend
2358  if (sm(i,j) /= 0.0) smc(i,j,3) = spval
2359  enddo
2360  enddo
2361  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,3)
2362 
2363  varname='soilw4'
2364  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2365  spval,varname,smc(ista_2l,jsta_2l,4))
2366 ! mask water areas
2367 !$omp parallel do private(i,j)
2368  do j=jsta,jend
2369  do i=ista,iend
2370  if (sm(i,j) /= 0.0) smc(i,j,4) = spval
2371  enddo
2372  enddo
2373  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,4)
2374 
2375  IF (nsoil==9) THEN
2376 
2377  varname='soilw5'
2378  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2379  spval,varname,smc(ista_2l,jsta_2l,5))
2380 ! mask water areas
2381 !$omp parallel do private(i,j)
2382  do j=jsta,jend
2383  do i=ista,iend
2384  if (sm(i,j) /= 0.0) smc(i,j,5) = spval
2385  enddo
2386  enddo
2387  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,5)
2388 
2389  varname='soilw6'
2390  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2391  spval,varname,smc(ista_2l,jsta_2l,6))
2392 ! mask water areas
2393 !$omp parallel do private(i,j)
2394  do j=jsta,jend
2395  do i=ista,iend
2396  if (sm(i,j) /= 0.0) smc(i,j,6) = spval
2397  enddo
2398  enddo
2399  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,6)
2400 
2401  varname='soilw7'
2402  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2403  spval,varname,smc(ista_2l,jsta_2l,7))
2404 ! mask water areas
2405 !$omp parallel do private(i,j)
2406  do j=jsta,jend
2407  do i=ista,iend
2408  if (sm(i,j) /= 0.0) smc(i,j,7) = spval
2409  enddo
2410  enddo
2411  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,7)
2412 
2413  varname='soilw8'
2414  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2415  spval,varname,smc(ista_2l,jsta_2l,8))
2416 ! mask water areas
2417 !$omp parallel do private(i,j)
2418  do j=jsta,jend
2419  do i=ista,iend
2420  if (sm(i,j) /= 0.0) smc(i,j,8) = spval
2421  enddo
2422  enddo
2423  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,8)
2424 
2425  varname='soilw9'
2426  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2427  spval,varname,smc(ista_2l,jsta_2l,9))
2428 ! mask water areas
2429 !$omp parallel do private(i,j)
2430  do j=jsta,jend
2431  do i=ista,iend
2432  if (sm(i,j) /= 0.0) smc(i,j,9) = spval
2433  enddo
2434  enddo
2435  if(debugprint)print*,'sample l',varname,' = ',1,smc(isa,jsa,9)
2436 
2437  END IF
2438 
2439 ! soil temperature using nemsio
2440  varname='soilt1'
2441  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2442  spval,varname,stc(ista_2l,jsta_2l,1))
2443 ! mask open water areas, combine with sea ice tmp
2444 !$omp parallel do private(i,j)
2445  do j=jsta,jend
2446  do i=ista,iend
2447  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval
2448  !if (sm(i,j) /= 0.0) stc(i,j,1) = spval
2449  enddo
2450  enddo
2451  if(debugprint)print*,'sample l','stc',' = ',1,stc(isa,jsa,1)
2452 
2453  varname='soilt2'
2454  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2455  spval,varname,stc(ista_2l,jsta_2l,2))
2456 ! mask open water areas, combine with sea ice tmp
2457 !$omp parallel do private(i,j)
2458  do j=jsta,jend
2459  do i=ista,iend
2460  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval
2461  !if (sm(i,j) /= 0.0) stc(i,j,2) = spval
2462  enddo
2463  enddo
2464  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,2)
2465 
2466  varname='soilt3'
2467  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2468  spval,varname,stc(ista_2l,jsta_2l,3))
2469 ! mask open water areas, combine with sea ice tmp
2470 !$omp parallel do private(i,j)
2471  do j=jsta,jend
2472  do i=ista,iend
2473  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval
2474  !if (sm(i,j) /= 0.0) stc(i,j,3) = spval
2475  enddo
2476  enddo
2477  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,3)
2478 
2479  varname='soilt4'
2480  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2481  spval,varname,stc(ista_2l,jsta_2l,4))
2482 ! mask open water areas, combine with sea ice tmp
2483 !$omp parallel do private(i,j)
2484  do j=jsta,jend
2485  do i=ista,iend
2486  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval
2487  !if (sm(i,j) /= 0.0) stc(i,j,4) = spval
2488  enddo
2489  enddo
2490  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,4)
2491 
2492  IF (nsoil==9) THEN
2493 
2494  varname='soilt5'
2495  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2496  spval,varname,stc(ista_2l,jsta_2l,5))
2497 ! mask open water areas, combine with sea ice tmp
2498 !$omp parallel do private(i,j)
2499  do j=jsta,jend
2500  do i=ista,iend
2501  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,5) = spval
2502  !if (sm(i,j) /= 0.0) stc(i,j,5) = spval
2503  enddo
2504  enddo
2505  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,5)
2506 
2507  varname='soilt6'
2508  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2509  spval,varname,stc(ista_2l,jsta_2l,6))
2510 ! mask open water areas, combine with sea ice tmp
2511 !$omp parallel do private(i,j)
2512  do j=jsta,jend
2513  do i=ista,iend
2514  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,6) = spval
2515  !if (sm(i,j) /= 0.0) stc(i,j,6) = spval
2516  enddo
2517  enddo
2518  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,6)
2519 
2520  varname='soilt7'
2521  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2522  spval,varname,stc(ista_2l,jsta_2l,7))
2523 ! mask open water areas, combine with sea ice tmp
2524 !$omp parallel do private(i,j)
2525  do j=jsta,jend
2526  do i=ista,iend
2527  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,7) = spval
2528  !if (sm(i,j) /= 0.0) stc(i,j,7) = spval enddo
2529  enddo
2530  enddo
2531  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,7)
2532 
2533  varname='soilt8'
2534  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2535  spval,varname,stc(ista_2l,jsta_2l,8))
2536 ! mask open water areas, combine with sea ice tmp
2537 !$omp parallel do private(i,j)
2538  do j=jsta,jend
2539  do i=ista,iend
2540  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,8) = spval
2541  !if (sm(i,j) /= 0.0) stc(i,j,8) = spval
2542  enddo
2543  enddo
2544  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,8)
2545 
2546  varname='soilt9'
2547  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2548  spval,varname,stc(ista_2l,jsta_2l,9))
2549 ! mask open water areas, combine with sea ice tmp
2550 !$omp parallel do private(i,j)
2551  do j=jsta,jend
2552  do i=ista,iend
2553  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,9) = spval
2554  !if (sm(i,j) /= 0.0) stc(i,j,9) = spval
2555  enddo
2556  enddo
2557  if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,9)
2558 
2559  END IF
2560 
2561 !$omp parallel do private(i,j)
2562  do j=jsta,jend
2563  do i=ista,iend
2564  acfrcv(i,j) = spval ! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1
2565  ncfrcv(i,j) = 1.0
2566  acfrst(i,j) = spval ! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1
2567  ncfrst(i,j) = 1.0
2568  bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF
2569  rlwtoa(i,j) = spval ! GFS does not have inst model top outgoing longwave
2570  enddo
2571  enddo
2572 ! trdlw(i,j) = 6.0
2573  ardlw = 1.0 ! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1
2574 
2575 ! time averaged incoming sfc longwave
2576  varname='dlwrf_ave'
2577  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2578  spval,varname,alwin)
2579 
2580 ! inst incoming sfc longwave
2581  varname='dlwrf'
2582  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2583  spval,varname,rlwin)
2584 
2585 ! time averaged outgoing sfc longwave
2586  varname='ulwrf_ave'
2587  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2588  spval,varname,alwout)
2589 
2590 ! inst outgoing sfc longwave
2591  varname='ulwrf'
2592  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2593  spval,varname,radot)
2594 
2595 ! where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing
2596 !$omp parallel do private(i,j)
2597  do j=jsta,jend
2598  do i=ista,iend
2599  if (alwout(i,j) /= spval) alwout(i,j) = -alwout(i,j)
2600  enddo
2601  enddo
2602 ! if(debugprint)print*,'sample l',VarName,' = ',1,alwout(isa,jsa)
2603 
2604 ! time averaged outgoing model top longwave using gfsio
2605  varname='ulwrf_avetoa'
2606  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2607  spval,varname,alwtoa)
2608 ! if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa)
2609 
2610 ! GFS incoming sfc longwave has been averaged, set ARDLW to 1
2611  ardsw=1.0
2612 ! trdsw=6.0
2613 
2614 ! time averaged incoming sfc shortwave
2615  varname='dswrf_ave'
2616  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2617  spval,varname,aswin)
2618 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa)
2619 
2620 ! inst incoming sfc shortwave
2621  varname='dswrf'
2622  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2623  spval,varname,rswin)
2624 
2625 ! inst incoming clear sky sfc shortwave
2626 ! FV3 do not output instant incoming clear sky sfc shortwave
2627  !$omp parallel do private(i,j)
2628  do j=jsta_2l,jend_2u
2629  do i=ista_2l,iend_2u
2630  rswinc(i,j) = spval
2631  enddo
2632  enddo
2633 
2634 ! time averaged incoming sfc uv-b using getgb
2635  varname='duvb_ave'
2636  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2637  spval,varname,auvbin)
2638 ! if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa)
2639 
2640 ! time averaged incoming sfc clear sky uv-b using getgb
2641  varname='cduvb_ave'
2642  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2643  spval,varname,auvbinc)
2644 ! if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa)
2645 
2646 ! time averaged outgoing sfc shortwave using gfsio
2647  varname='uswrf_ave'
2648  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2649  spval,varname,aswout)
2650 ! where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing
2651 !$omp parallel do private(i,j)
2652  do j=jsta,jend
2653  do i=ista,iend
2654  if (aswout(i,j) /= spval) aswout(i,j) = -aswout(i,j)
2655  enddo
2656  enddo
2657 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswout(isa,jsa)
2658 
2659 ! inst outgoing sfc shortwave using gfsio
2660  varname='uswrf'
2661  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2662  spval,varname,rswout)
2663 
2664 ! time averaged model top incoming shortwave
2665  varname='dswrf_avetoa'
2666  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2667  spval,varname,aswintoa)
2668 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa)
2669 
2670 ! time averaged model top outgoing shortwave
2671  varname='uswrf_avetoa'
2672  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2673  spval,varname,aswtoa)
2674 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(isa,jsa)
2675 
2676 ! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux
2677 ! has reversed sign convention using gfsio
2678  varname='shtfl_ave'
2679  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2680  spval,varname,sfcshx)
2681 ! where (sfcshx /= spval)sfcshx=-sfcshx
2682 !$omp parallel do private(i,j)
2683  do j=jsta,jend
2684  do i=ista,iend
2685  if (sfcshx(i,j) /= spval) sfcshx(i,j) = -sfcshx(i,j)
2686  enddo
2687  enddo
2688 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(isa,jsa)
2689 
2690 ! inst surface sensible heat flux
2691  varname='shtfl'
2692  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2693  spval,varname,twbs)
2694 !$omp parallel do private(i,j)
2695  do j=jsta,jend
2696  do i=ista,iend
2697  if (twbs(i,j) /= spval) twbs(i,j) = -twbs(i,j)
2698  enddo
2699  enddo
2700 
2701 ! GFS surface flux has been averaged, set ASRFC to 1
2702  asrfc=1.0
2703 ! tsrfc=6.0
2704 
2705 ! time averaged surface latent heat flux, multiplied by -1 because wrf model flux
2706 ! has reversed sign vonvention using gfsio
2707  varname='lhtfl_ave'
2708  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2709  spval,varname,sfclhx)
2710 ! where (sfclhx /= spval)sfclhx=-sfclhx
2711 !$omp parallel do private(i,j)
2712  do j=jsta,jend
2713  do i=ista,iend
2714  if (sfclhx(i,j) /= spval) sfclhx(i,j) = -sfclhx(i,j)
2715  enddo
2716  enddo
2717 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(isa,jsa)
2718 
2719 ! inst surface latent heat flux
2720  varname='lhtfl'
2721  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2722  spval,varname,qwbs)
2723 !$omp parallel do private(i,j)
2724  do j=jsta,jend
2725  do i=ista,iend
2726  if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j)
2727  enddo
2728  enddo
2729 
2730  if(me==0)print*,'rdaod= ',rdaod
2731 ! inst aod550 optical depth
2732  if(rdaod) then
2733  varname='aod550'
2734  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2735  spval,varname,aod550)
2736 
2737  varname='du_aod550'
2738  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2739  spval,varname,du_aod550)
2740 
2741  varname='ss_aod550'
2742  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2743  spval,varname,ss_aod550)
2744 
2745  varname='su_aod550'
2746  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2747  spval,varname,su_aod550)
2748 
2749  varname='oc_aod550'
2750  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2751  spval,varname,oc_aod550)
2752 
2753  varname='bc_aod550'
2754  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2755  spval,varname,bc_aod550)
2756  end if
2757 
2758 ! time averaged ground heat flux using nemsio
2759  varname='gflux_ave'
2760  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2761  spval,varname,subshx)
2762 ! mask water areas
2763 !$omp parallel do private(i,j)
2764  do j=jsta,jend
2765  do i=ista,iend
2766  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) subshx(i,j) = spval
2767  enddo
2768  enddo
2769 ! if(debugprint)print*,'sample l',VarName,' = ',1,subshx(isa,jsa)
2770 
2771 ! inst ground heat flux using nemsio
2772  varname='gflux'
2773  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2774  spval,varname,grnflx)
2775 ! mask water areas
2776 !$omp parallel do private(i,j)
2777  do j=jsta,jend
2778  do i=ista,iend
2779  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) grnflx(i,j) = spval
2780  enddo
2781  enddo
2782 
2783 ! time averaged zonal momentum flux using gfsio
2784  varname='uflx_ave'
2785  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2786  spval,varname,sfcux)
2787 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa)
2788 
2789 ! time averaged meridional momentum flux using nemsio
2790  varname='vflx_ave'
2791  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2792  spval,varname,sfcvx)
2793 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa)
2794 
2795 ! dong read in inst surface flux
2796 ! inst zonal momentum flux using gfsio
2797  varname='uflx'
2798  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2799  spval,varname,sfcuxi)
2800  if(debugprint)print*,'sample l',varname,' = ',1,sfcuxi(isa,jsa)
2801 
2802 ! inst meridional momentum flux using nemsio
2803  varname='vflx'
2804  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2805  spval,varname,sfcvxi)
2806  if(debugprint)print*,'sample l',varname,' = ',1,sfcvxi(isa,jsa)
2807 
2808 
2809 !$omp parallel do private(i,j)
2810  do j=jsta_2l,jend_2u
2811  do i=ista_2l,iend_2u
2812  sfcuvx(i,j) = spval ! GFS does not use total momentum flux
2813  enddo
2814  enddo
2815 
2816 ! time averaged zonal gravity wave stress using nemsio
2817  varname='u-gwd_ave'
2818  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2819  spval,varname,gtaux)
2820 ! if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa)
2821 
2822 ! time averaged meridional gravity wave stress using getgb
2823  varname='v-gwd_ave'
2824  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2825  spval,varname,gtauy)
2826 ! if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa)
2827 
2828 ! time averaged accumulated potential evaporation
2829  varname='pevpr_ave'
2830  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2831  spval,varname,avgpotevp)
2832 ! mask water areas
2833 !$omp parallel do private(i,j)
2834  do j=jsta,jend
2835  do i=ista,iend
2836  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) avgpotevp(i,j) = spval
2837  enddo
2838  enddo
2839 ! if(debugprint)print*,'sample l',VarName,' = ',1,potevp(isa,jsa)
2840 
2841 ! inst potential evaporation
2842  varname='pevpr'
2843  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2844  spval,varname,potevp)
2845 ! mask water areas
2846 !$omp parallel do private(i,j)
2847  do j=jsta,jend
2848  do i=ista,iend
2849  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) potevp(i,j) = spval
2850  enddo
2851  enddo
2852 
2853  do l=1,lm
2854 !$omp parallel do private(i,j)
2855  do j=jsta_2l,jend_2u
2856  do i=ista_2l,iend_2u
2857 ! GFS does not have temperature tendency due to long wave radiation
2858  rlwtt(i,j,l) = spval
2859 ! GFS does not have temperature tendency due to short wave radiation
2860  rswtt(i,j,l) = spval
2861 ! GFS does not have temperature tendency due to latent heating from convection
2862  tcucn(i,j,l) = spval
2863  tcucns(i,j,l) = spval
2864 ! GFS does not have temperature tendency due to latent heating from grid scale
2865  train(i,j,l) = spval
2866  enddo
2867  enddo
2868  enddo
2869 
2870 ! set avrain to 1
2871  avrain=1.0
2872  avcnvc=1.0
2873  theat=6.0 ! just in case GFS decides to output T tendency
2874 
2875 ! 10 m u using nemsio
2876  varname='ugrd10m'
2877  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2878  spval,varname,u10)
2879 
2880  do j=jsta,jend
2881  do i=ista,iend
2882  u10h(i,j)=u10(i,j)
2883  end do
2884  end do
2885 ! if(debugprint)print*,'sample l',VarName,' = ',1,u10(isa,jsa)
2886 
2887 ! 10 m v using gfsio
2888  varname='vgrd10m'
2889  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2890  spval,varname,v10)
2891 
2892  do j=jsta,jend
2893  do i=ista,iend
2894  v10h(i,j)=v10(i,j)
2895  end do
2896  end do
2897 ! if(debugprint)print*,'sample l',VarName,' = ',1,v10(isa,jsa)
2898 
2899 ! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon
2900  varname='vtype'
2901  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2902  spval,varname,buf)
2903 ! where (buf /= spval)
2904 ! ivgtyp=nint(buf)
2905 ! elsewhere
2906 ! ivgtyp=0 !need to feed reasonable value to crtm
2907 ! end where
2908 !$omp parallel do private(i,j)
2909  do j = jsta_2l, jend_2u
2910  do i=ista,iend
2911  if (buf(i,j) < spval) then
2912  ivgtyp(i,j) = nint(buf(i,j))
2913  else
2914  ivgtyp(i,j) = 0
2915  endif
2916  enddo
2917  enddo
2918 ! if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(isa,jsa)
2919 
2920 ! soil type, it's in GFS surface file, hopefully will merge into gfsio soon
2921  varname='sotyp'
2922  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2923  spval,varname,buf)
2924  l=1
2925 !$omp parallel do private(i,j)
2926  do j = jsta_2l, jend_2u
2927  do i=ista,iend
2928  if (buf(i,j) < spval) then
2929  isltyp(i,j) = nint(buf(i,j))
2930  else
2931  isltyp(i,j) = 0 !need to feed reasonable value to crtm
2932  endif
2933  enddo
2934  enddo
2935 ! if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(isa,jsa)
2936 
2937 !$omp parallel do private(i,j)
2938  do j=jsta_2l,jend_2u
2939  do i=ista_2l,iend_2u
2940  smstav(i,j) = spval ! GFS does not have soil moisture availability
2941 ! smstot(i,j) = spval ! GFS does not have total soil moisture
2942  sfcevp(i,j) = spval ! GFS does not have accumulated surface evaporation
2943  acsnow(i,j) = spval ! GFS does not have averaged accumulated snow
2944  acsnom(i,j) = spval ! GFS does not have snow melt
2945 ! sst(i,j) = spval ! GFS does not have sst????
2946  thz0(i,j) = ths(i,j) ! GFS does not have THZ0, use THS to substitute
2947  qz0(i,j) = spval ! GFS does not output humidity at roughness length
2948  uz0(i,j) = spval ! GFS does not output u at roughness length
2949  vz0(i,j) = spval ! GFS does not output humidity at roughness length
2950  enddo
2951  enddo
2952  do l=1,lm
2953 !$omp parallel do private(i,j)
2954  do j=jsta_2l,jend_2u
2955  do i=ista_2l,iend_2u
2956  el_pbl(i,j,l) = spval ! GFS does not have mixing length
2957  exch_h(i,j,l) = spval ! GFS does not output exchange coefficient
2958  enddo
2959  enddo
2960  enddo
2961 ! if(debugprint)print*,'sample l',VarName,' = ',1,thz0(isa,jsa)
2962 
2963 ! retrieve inst convective cloud top, GFS has cloud top pressure instead of index,
2964 ! will need to modify CLDRAD.f to use pressure directly instead of index
2965 ! VarName='pres'
2966 ! VcoordName='convect-cld top'
2967 ! l=1
2968 ! if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa)
2969  varname='prescnvclt'
2970  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
2971  spval,varname,ptop)
2972 
2973 
2974 !$omp parallel do private(i,j)
2975  do j=jsta,jend
2976  do i=ista,iend
2977  htop(i,j) = spval
2978  if(ptop(i,j) <= 0.0) ptop(i,j) = spval
2979  enddo
2980  enddo
2981  do j=jsta,jend
2982  do i=ista,iend
2983  if(ptop(i,j) < spval)then
2984  do l=1,lm
2985  if(ptop(i,j) <= pmid(i,j,l))then
2986  htop(i,j) = l
2987 ! if(i==ii .and. j==jj)print*,'sample ptop,pmid pmid-1,pint= ', &
2988 ! ptop(i,j),pmid(i,j,l),pmid(i,j,l-1),pint(i,j,l),htop(i,j)
2989  exit
2990  end if
2991  end do
2992  end if
2993  end do
2994  end do
2995 
2996 ! retrieve inst convective cloud bottom, GFS has cloud top pressure instead of index,
2997 ! will need to modify CLDRAD.f to use pressure directly instead of index
2998  varname='prescnvclb'
2999  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3000  spval,varname,pbot)
3001 ! if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa)
3002 !$omp parallel do private(i,j)
3003  do j=jsta,jend
3004  do i=ista,iend
3005  hbot(i,j) = spval
3006  if(pbot(i,j) <= 0.0) pbot(i,j) = spval
3007  enddo
3008  enddo
3009  do j=jsta,jend
3010  do i=ista,iend
3011 ! if(.not.lb(i,j))print*,'false bitmask for pbot at '
3012 ! + ,i,j,pbot(i,j)
3013  if(pbot(i,j) < spval)then
3014  do l=lm,1,-1
3015  if(pbot(i,j) >= pmid(i,j,l)) then
3016  hbot(i,j) = l
3017 ! if(i==ii .and. j==jj)print*,'sample pbot,pmid= ', &
3018 ! pbot(i,j),pmid(i,j,l),hbot(i,j)
3019  exit
3020  end if
3021  end do
3022  end if
3023  end do
3024  end do
3025  if(debugprint)print*,'sample hbot = ',hbot(isa,jsa)
3026 ! retrieve time averaged low cloud top pressure using nemsio
3027  varname='pres_avelct'
3028  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3029  spval,varname,ptopl)
3030 ! if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa)
3031 
3032 ! retrieve time averaged low cloud bottom pressure using nemsio
3033  varname='pres_avelcb'
3034  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3035  spval,varname,pbotl)
3036 ! if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa)
3037 
3038 ! retrieve time averaged low cloud top temperature using nemsio
3039  varname='tmp_avelct'
3040  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3041  spval,varname,ttopl)
3042 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa)
3043 
3044 ! retrieve time averaged middle cloud top pressure using nemsio
3045  varname='pres_avemct'
3046  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3047  spval,varname,ptopm)
3048 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa)
3049 
3050 ! retrieve time averaged middle cloud bottom pressure using nemsio
3051  varname='pres_avemcb'
3052  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3053  spval,varname,pbotm)
3054 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa)
3055 
3056 ! retrieve time averaged middle cloud top temperature using nemsio
3057  varname='tmp_avemct'
3058  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3059  spval,varname,ttopm)
3060 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa)
3061 
3062 ! retrieve time averaged high cloud top pressure using nemsio *********
3063  varname='pres_avehct'
3064  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3065  spval,varname,ptoph)
3066 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa)
3067 
3068 ! retrieve time averaged high cloud bottom pressure using nemsio
3069  varname='pres_avehcb'
3070  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3071  spval,varname,pboth)
3072 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa)
3073 
3074 ! retrieve time averaged high cloud top temperature using nemsio
3075  varname='tmp_avehct'
3076  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3077  spval,varname,ttoph)
3078 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa)
3079 
3080 ! retrieve boundary layer cloud cover using nemsio
3081  varname='tcdc_avebndcl'
3082  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3083  spval,varname,pblcfr)
3084 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa)
3085 ! where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction
3086 !$omp parallel do private(i,j)
3087  do j = jsta_2l, jend_2u
3088  do i=ista,iend
3089  if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01
3090  enddo
3091  enddo
3092 
3093 ! retrieve cloud work function
3094  varname='cwork_aveclm'
3095  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3096  spval,varname,cldwork)
3097 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa)
3098 
3099 ! accumulated total (base+surface) runoff
3100  varname='watr_acc'
3101  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3102  spval,varname,runoff)
3103 ! mask water areas
3104 !$omp parallel do private(i,j)
3105  do j=jsta,jend
3106  do i=ista,iend
3107  if (sm(i,j) /= 0.0) runoff(i,j) = spval
3108  enddo
3109  enddo
3110 
3111 ! total water storage in aquifer
3112  varname='wa_acc'
3113  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3114  spval,varname,twa)
3115 ! mask water areas
3116 !$omp parallel do private(i,j)
3117  do j=jsta,jend
3118  do i=ista,iend
3119  if (sm(i,j) /= 0.0) twa(i,j) = spval
3120  enddo
3121  enddo
3122 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa)
3123 
3124 ! accumulated evaporation of intercepted water
3125  varname='ecan_acc'
3126  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3127  spval,varname,tecan)
3128 ! mask water areas
3129 !$omp parallel do private(i,j)
3130  do j=jsta,jend
3131  do i=ista,iend
3132  if (sm(i,j) /= 0.0) tecan(i,j) = spval
3133  enddo
3134  enddo
3135 
3136 ! accumulated plant transpiration
3137  varname='etran_acc'
3138  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3139  spval,varname,tetran)
3140 ! mask water areas
3141 !$omp parallel do private(i,j)
3142  do j=jsta,jend
3143  do i=ista,iend
3144  if (sm(i,j) /= 0.0) tetran(i,j) = spval
3145  enddo
3146  enddo
3147 
3148 ! accumulated soil surface evaporation
3149  varname='edir_acc'
3150  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3151  spval,varname,tedir)
3152 ! mask water areas
3153 !$omp parallel do private(i,j)
3154  do j=jsta,jend
3155  do i=ista,iend
3156  if (sm(i,j) /= 0.0) tedir(i,j) = spval
3157  enddo
3158  enddo
3159 
3160 ! retrieve shelter max temperature using nemsio
3161  varname='t02max'
3162  if(modelname=='GFS') varname='tmax_max2m'
3163  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3164  spval,varname,maxtshltr)
3165 
3166 ! retrieve shelter min temperature using nemsio
3167  varname='t02min'
3168  if(modelname=='GFS') varname='tmin_min2m'
3169  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3170  spval,varname,mintshltr)
3171 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
3172 ! 1,mintshltr((ista+iend)/2,(jsta+jend)/2)
3173 
3174 ! retrieve shelter max RH
3175  varname='rh02max'
3176  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3177  spval,varname,maxrhshltr)
3178 
3179 ! retrieve shelter min temperature using nemsio
3180  varname='rh02min'
3181  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3182  spval,varname,minrhshltr)
3183 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
3184 ! 1,mintshltr((ista+iend)/2,(jsta+jend)/2)
3185 
3186 ! retrieve shelter max specific humidity using nemsio
3187  varname='spfhmax_max2m'
3188  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3189  spval,varname,maxqshltr)
3190 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',
3191 ! 1,maxqshltr(isa,jsa)
3192 
3193 ! retrieve shelter min temperature using nemsio
3194  varname='spfhmin_min2m'
3195  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3196  spval,varname,minqshltr)
3197 
3198 ! retrieve ice thickness using nemsio
3199  varname='icetk'
3200  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3201  spval,varname,dzice)
3202 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa)
3203 
3204 ! retrieve wilting point using nemsio
3205  varname='wilt'
3206  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3207  spval,varname,smcwlt)
3208 ! mask water areas
3209 !$omp parallel do private(i,j)
3210  do j=jsta,jend
3211  do i=ista,iend
3212  if (sm(i,j) /= 0.0) smcwlt(i,j) = spval
3213  enddo
3214  enddo
3215 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,smcwlt(isa,jsa)
3216 
3217 ! retrieve sunshine duration using nemsio
3218  varname='sunsd_acc'
3219  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3220  spval,varname,suntime)
3221 
3222 ! retrieve field capacity using nemsio
3223  varname='fldcp'
3224  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3225  spval,varname,fieldcapa)
3226 ! mask water areas
3227 !$omp parallel do private(i,j)
3228  do j=jsta,jend
3229  do i=ista,iend
3230  if (sm(i,j) /= 0.0) fieldcapa(i,j) = spval
3231  enddo
3232  enddo
3233 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,fieldcapa(isa,jsa)
3234 
3235 ! retrieve time averaged surface visible beam downward solar flux
3236  varname='vbdsf_ave'
3237  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3238  spval,varname,avisbeamswin)
3239  vcoordname='sfc'
3240  l=1
3241 
3242 ! retrieve time averaged surface visible diffuse downward solar flux
3243  varname='vddsf_ave'
3244  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3245  spval,varname,avisdiffswin)
3246 
3247 ! retrieve time averaged surface near IR beam downward solar flux
3248  varname='nbdsf_ave'
3249  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3250  spval,varname,airbeamswin)
3251 
3252 ! retrieve time averaged surface near IR diffuse downward solar flux
3253  varname='nddsf_ave'
3254  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3255  spval,varname,airdiffswin)
3256 
3257 ! retrieve time averaged surface clear sky outgoing LW
3258  varname='csulf'
3259  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3260  spval,varname,alwoutc)
3261 
3262 ! retrieve time averaged TOA clear sky outgoing LW
3263  varname='csulftoa'
3264  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3265  spval,varname,alwtoac)
3266 
3267 ! retrieve time averaged surface clear sky outgoing SW
3268  varname='csusf'
3269  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3270  spval,varname,aswoutc)
3271 
3272 ! retrieve time averaged TOA clear sky outgoing LW
3273  varname='csusftoa'
3274  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3275  spval,varname,aswtoac)
3276 
3277 ! retrieve time averaged surface clear sky incoming LW
3278  varname='csdlf'
3279  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3280  spval,varname,alwinc)
3281 
3282 ! retrieve time averaged surface clear sky incoming SW
3283  varname='csdsf'
3284  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3285  spval,varname,aswinc)
3286 
3287 ! retrieve storm runoff using nemsio
3288  varname='ssrun_acc'
3289  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3290  spval,varname,ssroff)
3291 ! mask water areas
3292 !$omp parallel do private(i,j)
3293  do j=jsta,jend
3294  do i=ista,iend
3295  if (sm(i,j) /= 0.0) ssroff(i,j) = spval
3296  enddo
3297  enddo
3298 
3299 ! retrieve direct soil evaporation
3300  varname='evbs_ave'
3301  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3302  spval,varname,avgedir)
3303 ! mask water areas
3304 !$omp parallel do private(i,j)
3305  do j=jsta,jend
3306  do i=ista,iend
3307  if (sm(i,j) /= 0.0) avgedir(i,j) = spval
3308  enddo
3309  enddo
3310 
3311 ! retrieve CANOPY WATER EVAP
3312  varname='evcw_ave'
3313  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3314  spval,varname,avgecan)
3315 ! mask water areas
3316 !$omp parallel do private(i,j)
3317  do j=jsta,jend
3318  do i=ista,iend
3319  if (sm(i,j) /= 0.0) avgecan(i,j) = spval
3320  enddo
3321  enddo
3322 
3323 ! retrieve AVERAGED PRECIP ADVECTED HEAT FLUX
3324  varname='pah_ave'
3325  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3326  spval,varname,paha)
3327 ! mask water areas
3328 !$omp parallel do private(i,j)
3329  do j=jsta,jend
3330  do i=ista,iend
3331  if (sm(i,j) /= 0.0) paha(i,j) = spval
3332  enddo
3333  enddo
3334 
3335 ! retrieve instantaneous PRECIP ADVECTED HEAT FLUX
3336  varname='pahi'
3337  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3338  spval,varname,pahi)
3339 ! mask water areas
3340 !$omp parallel do private(i,j)
3341  do j=jsta,jend
3342  do i=ista,iend
3343  if (sm(i,j) /= 0.0) pahi(i,j) = spval
3344  enddo
3345  enddo
3346 
3347 ! retrieve PLANT TRANSPIRATION
3348  varname='trans_ave'
3349  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3350  spval,varname,avgetrans)
3351 ! mask water areas
3352 !$omp parallel do private(i,j)
3353  do j=jsta,jend
3354  do i=ista,iend
3355  if (sm(i,j) /= 0.0) avgetrans(i,j) = spval
3356  enddo
3357  enddo
3358 
3359 ! retrieve snow sublimation
3360  varname='sbsno_ave'
3361  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3362  spval,varname,avgesnow)
3363 ! mask water areas
3364 !$omp parallel do private(i,j)
3365  do j=jsta,jend
3366  do i=ista,iend
3367  if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval
3368  enddo
3369  enddo
3370 
3371 ! retrive total soil moisture
3372  varname='soilm'
3373  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3374  spval,varname,smstot)
3375 ! mask water areas
3376 !$omp parallel do private(i,j)
3377  do j=jsta,jend
3378  do i=ista,iend
3379  if (sm(i,j) /= 0.0) smstot(i,j) = spval
3380  enddo
3381  enddo
3382 
3383 ! retrieve snow phase change heat flux
3384  varname='snohf'
3385  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3386  spval,varname,snopcx)
3387 ! mask water areas
3388 !$omp parallel do private(i,j)
3389  do j=jsta,jend
3390  do i=ista,iend
3391  if (sm(i,j) /= 0.0) snopcx(i,j) = spval
3392  enddo
3393  enddo
3394 
3395 ! retrieve pwater
3396  varname='pwat'
3397  call read_netcdf_2d_para(ncid2d,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3398  spval,varname,pwat)
3399 
3400 ! GFS does not have deep convective cloud top and bottom fields
3401 
3402 !$omp parallel do private(i,j)
3403  do j=jsta,jend
3404  do i=ista,iend
3405  htopd(i,j) = spval
3406  hbotd(i,j) = spval
3407  htops(i,j) = spval
3408  hbots(i,j) = spval
3409  cuppt(i,j) = spval
3410  enddo
3411  enddo
3412 
3413 ! done with flux file, close it for now
3414  status=nf90_close(ncid2d)
3415 ! deallocate(tmp,recname,reclevtyp,reclev)
3416 
3417 ! pos east
3418 ! call collect_loc(gdlat,dummy)
3419 ! if(me == 0)then
3420 ! latstart = nint(dummy(1,1)*gdsdegr)
3421 ! latlast = nint(dummy(im,jm)*gdsdegr)
3422 ! print*,'laststart,latlast B bcast= ',latstart,latlast,'gdsdegr=',gdsdegr,&
3423 ! 'dummy(1,1)=',dummy(1,1),dummy(im,jm),'gdlat=',gdlat(1,1)
3424 ! end if
3425 ! call mpi_bcast(latstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
3426 ! call mpi_bcast(latlast,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
3427 ! write(6,*) 'laststart,latlast,me A calling bcast=',latstart,latlast,me
3428 ! call collect_loc(gdlon,dummy)
3429 ! if(me == 0)then
3430 ! lonstart = nint(dummy(1,1)*gdsdegr)
3431 ! lonlast = nint(dummy(im,jm)*gdsdegr)
3432 ! end if
3433 ! call mpi_bcast(lonstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
3434 ! call mpi_bcast(lonlast, 1,MPI_INTEGER,0,mpi_comm_comp,irtn)
3435 
3436 ! write(6,*)'lonstart,lonlast A calling bcast=',lonstart,lonlast
3437 !
3438 
3439 ! generate look up table for lifted parcel calculations
3440 
3441  thl = 210.
3442  plq = 70000.
3443  pt_tbl = 10000. ! this is for 100 hPa added by Moorthi
3444 
3445  CALL table(ptbl,ttbl,pt_tbl, &
3446  rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
3447 
3448  CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
3449 
3450 !
3451 !
3452  IF(me == 0)THEN
3453  WRITE(6,*)' SPL (POSTED PRESSURE LEVELS) BELOW: '
3454  WRITE(6,51) (spl(l),l=1,lsm)
3455  50 FORMAT(14(f4.1,1x))
3456  51 FORMAT(8(f8.1,1x))
3457  ENDIF
3458 !
3459 !$omp parallel do private(l)
3460  DO l = 1,lsm
3461  alsl(l) = log(spl(l))
3462  END DO
3463 !
3464 !HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
3465  if(me == 0)then
3466  print*,'writing out igds'
3467  igdout = 110
3468 ! open(igdout,file='griddef.out',form='unformatted'
3469 ! + ,status='unknown')
3470  if(maptype == 1)THEN ! Lambert conformal
3471  WRITE(igdout)3
3472  WRITE(6,*)'igd(1)=',3
3473  WRITE(igdout)im
3474  WRITE(igdout)jm
3475  WRITE(igdout)latstart
3476  WRITE(igdout)lonstart
3477  WRITE(igdout)8
3478  WRITE(igdout)cenlon
3479  WRITE(igdout)dxval
3480  WRITE(igdout)dyval
3481  WRITE(igdout)0
3482  WRITE(igdout)64
3483  WRITE(igdout)truelat2
3484  WRITE(igdout)truelat1
3485  WRITE(igdout)255
3486  ELSE IF(maptype == 2)THEN !Polar stereographic
3487  WRITE(igdout)5
3488  WRITE(igdout)im
3489  WRITE(igdout)jm
3490  WRITE(igdout)latstart
3491  WRITE(igdout)lonstart
3492  WRITE(igdout)8
3493  WRITE(igdout)cenlon
3494  WRITE(igdout)dxval
3495  WRITE(igdout)dyval
3496  WRITE(igdout)0
3497  WRITE(igdout)64
3498  WRITE(igdout)truelat2 !Assume projection at +-90
3499  WRITE(igdout)truelat1
3500  WRITE(igdout)255
3501  ! Note: The calculation of the map scale factor at the standard
3502  ! lat/lon and the PSMAPF
3503  ! Get map factor at 60 degrees (N or S) for PS projection, which will
3504  ! be needed to correctly define the DX and DY values in the GRIB GDS
3505  if (truelat1 < 0.) THEN
3506  lat = -60.
3507  else
3508  lat = 60.
3509  end if
3510 
3511  CALL msfps(lat,truelat1*0.001,psmapf)
3512 
3513  ELSE IF(maptype == 3) THEN !Mercator
3514  WRITE(igdout)1
3515  WRITE(igdout)im
3516  WRITE(igdout)jm
3517  WRITE(igdout)latstart
3518  WRITE(igdout)lonstart
3519  WRITE(igdout)8
3520  WRITE(igdout)latlast
3521  WRITE(igdout)lonlast
3522  WRITE(igdout)truelat1
3523  WRITE(igdout)0
3524  WRITE(igdout)64
3525  WRITE(igdout)dxval
3526  WRITE(igdout)dyval
3527  WRITE(igdout)255
3528  ELSE IF(maptype == 0 .OR. maptype == 203)THEN !A STAGGERED E-GRID
3529  WRITE(igdout)203
3530  WRITE(igdout)im
3531  WRITE(igdout)jm
3532  WRITE(igdout)latstart
3533  WRITE(igdout)lonstart
3534  WRITE(igdout)136
3535  WRITE(igdout)cenlat
3536  WRITE(igdout)cenlon
3537  WRITE(igdout)dxval
3538  WRITE(igdout)dyval
3539  WRITE(igdout)64
3540  WRITE(igdout)0
3541  WRITE(igdout)0
3542  WRITE(igdout)0
3543  END IF
3544  end if
3545 !
3546 !
3547 
3548  RETURN
3549  END
3550 
3551 
3552  subroutine read_netcdf_3d_para(ncid,im,jm,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3553  spval,varname,buf,lm)
3554 
3555  use netcdf
3556  use ctlblk_mod, only : me
3557  use params_mod, only : small
3558  implicit none
3559  include "mpif.h"
3560 
3561  character(len=20),intent(in) :: varname
3562  real,intent(in) :: spval
3563  integer,intent(in) :: ncid,im,jm,lm,jsta_2l,jend_2u,jsta,jend
3564  integer,intent(in) :: ista_2l,iend_2u,ista,iend
3565  real,intent(out) :: buf(ista_2l:iend_2u,jsta_2l:jend_2u,lm)
3566  integer :: varid,iret,ii,jj,i,j,l,kk
3567  integer :: start(3), count(3), stride(3)
3568  real,parameter :: spval_netcdf=9.99e+20
3569  real :: fill_value
3570 
3571  iret = nf90_inq_varid(ncid,trim(varname),varid)
3572  if (iret /= 0) then
3573  if (me == 0) print*,varname," not found -Assigned missing values"
3574 !$omp parallel do private(i,j,l)
3575  do l=1,lm
3576  do j=jsta,jend
3577  do i=ista,iend
3578  buf(i,j,l)=spval
3579  enddo
3580  enddo
3581  enddo
3582  else
3583  iret = nf90_get_att(ncid,varid,"_FillValue",fill_value)
3584  if (iret /= 0) fill_value = spval_netcdf
3585  start = (/ista,jsta,1/)
3586  ii=iend-ista+1
3587  jj=jend-jsta+1
3588  count = (/ii,jj,lm/)
3589  iret = nf90_get_var(ncid,varid,buf(ista:iend,jsta:jend,1:lm),start=start,count=count)
3590  if (iret /= 0) then
3591  print*," iret /=0, Error in reading varid "
3592  endif
3593  do l=1,lm
3594  do j=jsta,jend
3595  do i=ista,iend
3596  if(abs(buf(i,j,l)-fill_value)<small)buf(i,j,l)=spval
3597  end do
3598  end do
3599  end do
3600  endif
3601 
3602  end subroutine read_netcdf_3d_para
3603 
3604  subroutine read_netcdf_2d_para(ncid,ista,ista_2l,iend,iend_2u,jsta,jsta_2l,jend,jend_2u, &
3605  spval,varname,buf)
3606 
3607  use netcdf
3608  use ctlblk_mod, only : me
3609  use params_mod, only : small
3610  implicit none
3611  include "mpif.h"
3612 
3613  character(len=20),intent(in) :: varname
3614  real,intent(in) :: spval
3615  integer,intent(in) :: ncid,jsta_2l,jend_2u,jsta,jend,ista_2l,iend_2u,ista,iend
3616  real,intent(out) :: buf(ista_2l:iend_2u,jsta_2l:jend_2u)
3617  integer :: varid,iret,ii,jj,i,j,l,kk
3618  integer :: start(2), count(2)
3619  real,parameter :: spval_netcdf=9.99e+20
3620  real :: fill_value
3621 
3622  iret = nf90_inq_varid(ncid,trim(varname),varid)
3623  if (iret /= 0) then
3624  if (me==0) print*,varname," not found -Assigned missing values"
3625 !$omp parallel do private(i,j)
3626  do j=jsta,jend
3627  do i=ista,iend
3628  buf(i,j)=spval
3629  enddo
3630  enddo
3631  else
3632  iret = nf90_get_att(ncid,varid,"_FillValue",fill_value)
3633  if (iret /= 0) fill_value = spval_netcdf
3634  start = (/ista,jsta/)
3635  ii=iend-ista+1
3636  jj=jend-jsta+1
3637  count = (/ii,jj/)
3638  iret = nf90_get_var(ncid,varid,buf(ista:iend,jsta:jend),start=start,count=count)
3639  if (iret /= 0) then
3640  print*," iret /=0, Error in reading varid "
3641  endif
3642  do j=jsta,jend
3643  do i=ista,iend
3644  if(abs(buf(i,j)-fill_value)<small)buf(i,j)=spval
3645  end do
3646  end do
3647  endif
3648 
3649  end subroutine read_netcdf_2d_para
masks
Definition: MASKS_mod.f:1
physcons_post
Definition: physcons.f:1
params_mod
Definition: params.F:1
soil
Definition: SOIL_mod.f:1
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:345
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27
lookup_mod
Definition: LOOKUP.f:1