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