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