UPP  001
 All Data Structures Files Functions Pages
INITPOST_GFS_NEMS_MPIIO.f
Go to the documentation of this file.
1 
5 
29  SUBROUTINE initpost_gfs_nems_mpiio(iostatusAER)
30 
31 
32  use vrbls4d, only: dust, salt, suso, soot, waso, pp25, pp10
33  use vrbls3d, only: t, q, uh, vh,wh,pmid,pint,alpint, dpres,zint,zmid,o3, &
34  qqr, qqs, cwm, qqi, qqw, omga, rhomid, q2, cfr, rlwtt, rswtt, tcucn, &
35  tcucns, train, el_pbl, exch_h, vdifftt, vdiffmois, dconvmois, nradtt, &
36  o3vdiff, o3prod, o3tndy, mwpv, qqg, vdiffzacce, zgdrag,cnvctummixing, &
37  vdiffmacce, mgdrag, cnvctvmmixing, ncnvctcfrac, cnvctumflx, cnvctdmflx, &
38  cnvctzgdrag, sconvmois, cnvctmgdrag, cnvctdetmflx, duwt, duem, dusd, dudp, &
39  dusv,ssem,sssd,ssdp,sswt,sssv,bcem,bcsd,bcdp,bcwt,bcsv,ocem,ocsd,ocdp, &
40  ocwt,ocsv, ref_10cm
41  use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate, &
42  cprate, avgprec, prec, lspa, sno, si, cldefi, th10, q10, tshltr, pshltr, &
43  tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, radot, sigt4, &
44  cfrach, cfracl, cfracm, avgcfrach, qshltr, avgcfracl, avgcfracm, cnvcfr, &
45  islope, cmc, grnflx, vegfrc, acfrcv, ncfrcv, acfrst, ncfrst, ssroff, &
46  bgroff, rlwin, rlwtoa, cldwork, alwin, alwout, alwtoa, rswin, rswinc, &
47  rswout, aswin, auvbin, auvbinc, aswout, aswtoa, sfcshx, sfclhx, subshx, &
48  snopcx, sfcux, sfcvx, sfcuvx, gtaux, gtauy, potevp, u10, v10, smstav, &
49  smstot, ivgtyp, isltyp, sfcevp, sfcexc, acsnow, acsnom, sst, thz0, qz0, &
50  uz0, vz0, ptop, htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm, pbotm, ttopm, &
51  ptoph, pboth, pblcfr, ttoph, runoff, maxtshltr, mintshltr, maxrhshltr, &
52  minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots, &
53  cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa, &
54  maxqshltr, minqshltr, acond, sr, u10h, v10h, &
55  avgedir,avgecan,avgetrans,avgesnow,avgprec_cont,avgcprate_cont, &
56  avisbeamswin,avisdiffswin,airbeamswin,airdiffswin, &
57  alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, &
58  dustcb,bccb,occb,sulfcb,sscb,dustallcb,ssallcb,dustpm,dustpm10,sspm,pp25cb, &
59  pp10cb,maod,ti
60  use soil, only: sldpth, sh2o, smc, stc
61  use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
62 ! use kinds, only: i_llong
63 ! use nemsio_module, only: nemsio_gfile, nemsio_getfilehead, nemsio_getheadvar, nemsio_close
64  use physcons_post, only: grav => con_g, fv => con_fvirt, rgas => con_rd, &
65  eps => con_eps, epsm1 => con_epsm1
66  use params_mod, only: erad, dtr, tfrz, h1, d608, rd, p1000, capa
67  use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe, &
68  ttblq, rdpq, rdtheq, stheq, the0q, the0
69  use ctlblk_mod, only: me, mpi_comm_comp, icnt, idsp, jsta, jend, ihrst, idat, sdat, ifhr, &
70  ifmin, filename, tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, td3d, restrt, sdat, &
71  jend_m, imin, imp_physics, dt, spval, pdtop, pt, qmin, nbin_du, nphs, dtq2, ardlw,&
72  ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm, &
73  jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc, &
74  nbin_oc, nbin_su, gocart_on, pt_tbl, hyb_sigp, filenameflux, filenameaer, &
75  isf_surface_physics
76  use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, &
77  dxval, dyval, truelat2, truelat1, psmapf, cenlat
78  use nemsio_module_mpi
79  use upp_physics, only: fpvsnew
80 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81  implicit none
82 !
83  type(nemsio_gfile) :: nfile,ffile,rfile
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  integer,intent(in) :: iostatusaer
104  character(len=20) :: varname, vcoordname
105  integer :: status, fldsize, fldst, recn
106  integer :: recn_vvel,recn_delz,recn_dpres
107  character startdate*19,sysdepinfo*80,cgar*1
108  character startdate2(19)*4,lprecip_accu*3
109 !
110 ! NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
111 ! AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
112 !
113 ! ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
114 ! INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
115  LOGICAL runb,singlrst,subpost,nest,hydro,ioomg,ioall
116  logical, parameter :: debugprint = .false., zerout = .false.
117 ! logical, parameter :: debugprint = .true., zerout = .false.
118  logical :: reduce_grid = .true. ! set default to true for ops GSM
119  CHARACTER*32 label
120  CHARACTER*40 contrl,filall,filmst,filtmp,filtke,filunv,filcld,filrad,filsfc
121  CHARACTER*4 resthr
122  CHARACTER fname*255,envar*50
123  INTEGER idate(8),jdate(8),jpds(200),jgds(200),kpds(200),kgds(200)
124 ! LOGICAL*1 LB(IM,JM)
125 !
126 ! INCLUDE COMMON BLOCKS.
127 !
128 ! DECLARE VARIABLES.
129 !
130 ! REAL fhour
131  integer nfhour ! forecast hour from nems io file
132  integer fhzero ! bucket
133  real dtp !physics time step
134  REAL rinc(5)
135 
136  REAL dummy(im,jm)
137  real, allocatable :: fi(:,:,:)
138 !jw
139  integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus, &
140  i,j,l,ll,k,kf,irtn,igdout,n,index,nframe, &
141  impf,jmpf,nframed2,iunitd3d,ierr,idum,iret,nrec,idrt
142  real tstart,tlmh,tsph,es,fact,soilayert,soilayerb,zhour,dum, &
143  tvll,pmll,tv, tx1, tx2
144 
145  character*16,allocatable :: recname(:)
146  character*16,allocatable :: reclevtyp(:)
147  character*6 :: modelname_nemsio
148  integer, allocatable :: reclev(:), kmsk(:,:)
149  real, allocatable :: glat1d(:), glon1d(:), qstl(:)
150  real, allocatable :: wrk1(:,:), wrk2(:,:)
151  real, allocatable :: p2d(:,:), t2d(:,:), q2d(:,:), &
152  qs2d(:,:), cw2d(:,:), cfr2d(:,:)
153  real(kind=4),allocatable :: vcoord4(:,:,:)
154  real, dimension(lm+1) :: ak5, bk5
155  real*8, allocatable :: pm2d(:,:), pi2d(:,:)
156  real, allocatable :: tmp(:)
157  real :: buf(im,jsta_2l:jend_2u)
158  integer :: lonsperlat(jm/2), numi(jm)
159 
160 ! real buf(im,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u) &
161 ! ,buf3d(im,jsta_2l:jend_2u,lm),buf3d2(im,lp1,jsta_2l:jend_2u)
162 
163  real lat
164  integer isa, jsa, latghf, jtem, idvc, idsl, nvcoord, ip1, nn, npass
165 ! REAL, PARAMETER :: QMIN = 1.E-15
166 
167 ! DATA BLANK/' '/
168 !
169 ! integer, parameter :: npass2=2, npass3=3
170 ! integer, parameter :: npass2=20, npass3=30
171  integer, parameter :: npass2=5, npass3=30
172  real, parameter :: third=1.0/3.0
173  INTEGER, DIMENSION(2) :: ij4min, ij4max
174  REAL :: omgmin, omgmax
175  real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:)
176  REAL, ALLOCATABLE :: ps2d(:,:),psx2d(:,:),psy2d(:,:)
177  real, allocatable :: div3d(:,:,:)
178  real(kind=4),allocatable :: vcrd(:,:)
179  real :: omg1(im), omg2(im+2)
180 
181 !***********************************************************************
182 ! START INIT HERE.
183 !
184  WRITE(6,*)'INITPOST: ENTER INITPOST_GFS_NEMS_MPIIO'
185  WRITE(6,*)'me=',me,'LMV=',size(lmv,1),size(lmv,2),'LMH=', &
186  size(lmh,1),size(lmh,2),'jsta_2l=',jsta_2l,'jend_2u=', &
187  jend_2u,'im=',im
188 !
189  isa = im / 2
190  jsa = (jsta+jend) / 2
191 
192 !$omp parallel do private(i,j)
193  do j = jsta_2l, jend_2u
194  do i=1,im
195  buf(i,j) = spval
196  enddo
197  enddo
198 
199 ! initialize nemsio using mpi io module
200  call nemsio_init()
201  call nemsio_open(nfile,trim(filename),'read',mpi_comm_comp,iret=status)
202  if ( status /= 0 ) then
203  print*,'error opening ',filename, ' Status = ', status ; stop
204  endif
205  call nemsio_getfilehead(nfile,iret=status,nrec=nrec,idrt=idrt)
206 
207 ! open flux file early yo read imp_physics
208  call nemsio_open(ffile,trim(filenameflux),'read',mpi_comm_comp &
209  ,iret=status)
210  if ( status /= 0 ) then
211  print*,'error opening ',filenameflux, ' Status = ', status
212  endif
213 
214 !
215 ! STEP 1. READ MODEL OUTPUT FILE
216 !
217 ! LMH and LMV always = LM for sigma-type vert coord
218 
219 !$omp parallel do private(i,j)
220  do j = jsta_2l, jend_2u
221  do i = 1, im
222  lmv(i,j) = lm
223  lmh(i,j) = lm
224  end do
225  end do
226 
227 ! HTM VTM all 1 for sigma-type vert coord
228 
229 !$omp parallel do private(i,j,l)
230  do l = 1, lm
231  do j = jsta_2l, jend_2u
232  do i = 1, im
233  htm(i,j,l) = 1.0
234  vtm(i,j,l) = 1.0
235  end do
236  end do
237  end do
238 
239 ! write(0,*)'nrec=',nrec
240  allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
241  allocate(glat1d(im*jm),glon1d(im*jm))
242  allocate(vcoord4(lm+1,3,2))
243 ! get start date
244  call nemsio_getfilehead(nfile,iret=iret &
245  ,idate=idate(1:7),nfhour=nfhour,recname=recname &
246  ,reclevtyp=reclevtyp,reclev=reclev,lat=glat1d &
247  ,lon=glon1d,nframe=nframe,vcoord=vcoord4,idrt=maptype &
248  ,modelname=modelname_nemsio)
249  if(iret/=0)print*,'error getting idate,nfhour'
250  print *,'latstar1=',glat1d(1),glat1d(im*jm)
251 !
252 ! modelname_nemsio='FV3GFS'
253  print*,'modelname = ',modelname_nemsio
254  if(trim(modelname_nemsio)=='FV3GFS')reduce_grid=.false.
255  if(reduce_grid)then
256 !------------------------------
257  if (idrt == 4) then
258 !------------------------------
259 ! read lonsperlat
260  open (201,file='lonsperlat.dat',status='old',form='formatted', &
261  action='read',iostat=iret)
262  rewind(201)
263  read (201,*,iostat=iret) latghf,(lonsperlat(i),i=1,latghf)
264  close (201)
265  print*,'finished reading lonsperlat'
266 
267  if (jm /= latghf+latghf) then
268  write(0,*)' wrong reduced grid - execution skipped'
269  stop 777
270  endif
271  do j=1,jm/2
272  numi(j) = lonsperlat(j)
273  enddo
274  do j=jm/2+1,jm
275  numi(j) = lonsperlat(jm+1-j)
276  enddo
277 !------------------------------
278  else
279 !------------------------------
280  do j=1,jm
281  numi(j) = im
282  enddo
283 !------------------------------
284  endif
285 !------------------------------
286  end if
287 
288 ! Specigy grid staggering type
289  gridtype = 'A'
290  if (me == 0) print *, 'maptype and gridtype is ', &
291  maptype,gridtype
292 
293  if(debugprint)then
294  if (me == 0)then
295  do i=1,nrec
296  print *,'recname,reclevtyp,reclev=',trim(recname(i)),' ', &
297  trim(reclevtyp(i)),reclev(i)
298  end do
299  end if
300  end if
301 
302 !$omp parallel do private(i,j,js)
303  do j=jsta,jend
304  js = (j-1)*im
305  do i=1,im
306  gdlat(i,j) = glat1d(js+i)
307  gdlon(i,j) = glon1d(js+i)
308  end do
309  end do
310 !
311 ! if (hyb_sigp) then
312  do l=1,lm+1
313  ak5(l) = vcoord4(l,1,1)
314  bk5(l) = vcoord4(l,2,1)
315  enddo
316 ! endif
317 
318 !--Fanglin Yang: nemsio file created from FV3 does not have vcoord.
319  if ( minval(ak5) <0 .or. minval(bk5) <0 ) then
320  open (202,file='global_hyblev.txt',status='old',form='formatted', &
321  action='read',iostat=iret)
322  rewind(202)
323  read(202,*)
324  do l=1,lm+1
325  read (202,*,iostat=iret) ak5(l),bk5(l)
326  enddo
327  close (202)
328 
329  if (iret == 0 ) then
330  do l=1,lm+1
331  vcoord4(l,1,1)=ak5(l)
332  vcoord4(l,2,1)=bk5(l)
333  enddo
334  else
335  print *, 'ak5 and bk5 not found, stop !'
336  stop
337  endif
338  endif
339 
340  if (me == 0)then
341  print *,"ak5",ak5
342  print *,"bk5",bk5
343  endif
344 
345 ! deallocate(glat1d,glon1d,vcoord4)
346  deallocate(glat1d,glon1d)
347 
348  print*,'idate = ',(idate(i),i=1,7)
349  print*,'idate after broadcast = ',(idate(i),i=1,4)
350  print*,'nfhour = ',nfhour
351 
352 ! sample print point
353  ii = im/2
354  jj = jm/2
355 
356  print *,me,'max(gdlat)=', maxval(gdlat), &
357  'max(gdlon)=', maxval(gdlon)
358  CALL exch(gdlat(1,jsta_2l))
359  print *,'after call EXCH,me=',me
360 
361 !$omp parallel do private(i,j,ip1)
362  do j = jsta, jend_m
363  do i = 1, im
364  ip1 = i + 1
365  if (ip1 > im) ip1 = ip1 - im
366  dx(i,j) = erad*cos(gdlat(i,j)*dtr) *(gdlon(ip1,j)-gdlon(i,j))*dtr
367  dy(i,j) = erad*(gdlat(i,j)-gdlat(i,j+1))*dtr ! like A*DPH
368 ! F(I,J)=1.454441e-4*sin(gdlat(i,j)*DTR) ! 2*omeg*sin(phi)
369 ! if (i == ii .and. j == jj) print*,'sample LATLON, DY, DY=' &
370 ! ,i,j,GDLAT(I,J),GDLON(I,J),DX(I,J),DY(I,J)
371  end do
372  end do
373 
374 !$omp parallel do private(i,j)
375  do j=jsta,jend
376  do i=1,im
377  f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr) ! 2*omeg*sin(phi)
378  end do
379  end do
380 
381  impf = im
382  jmpf = jm
383  print*,'impf,jmpf,nframe= ',impf,jmpf,nframe
384 
385  iyear = idate(1)
386  imn = idate(2) ! ask Jun
387  iday = idate(3) ! ask Jun
388  ihrst = idate(4)
389  imin = idate(5)
390  jdate = 0
391  idate = 0
392 !
393  print*,'start yr mo day hr min =',iyear,imn,iday,ihrst,imin
394  print*,'processing yr mo day hr min=' &
395  ,idat(3),idat(1),idat(2),idat(4),idat(5)
396 !
397  idate(1) = iyear
398  idate(2) = imn
399  idate(3) = iday
400  idate(5) = ihrst
401  idate(6) = imin
402  sdat(1) = imn
403  sdat(2) = iday
404  sdat(3) = iyear
405  jdate(1) = idat(3)
406  jdate(2) = idat(1)
407  jdate(3) = idat(2)
408  jdate(5) = idat(4)
409  jdate(6) = idat(5)
410 !
411  print *,' idate=',idate
412  print *,' jdate=',jdate
413 !
414  CALL w3difdat(jdate,idate,0,rinc)
415 !
416  print *,' rinc=',rinc
417  ifhr = nint(rinc(2)+rinc(1)*24.)
418  print *,' ifhr=',ifhr
419  ifmin = nint(rinc(3))
420 ! if(ifhr /= nint(fhour))print*,'find wrong Grib file';stop
421  print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,filename
422 
423 ! Getting tstart
424  tstart = 0.
425  print*,'tstart= ',tstart
426 
427 ! Getiing restart
428 
429  restrt = .true. ! set RESTRT as default
430 
431  IF(tstart > 1.0e-2)THEN
432  ifhr = ifhr+nint(tstart)
433  rinc = 0
434  idate = 0
435  rinc(2) = -1.0*ifhr
436  call w3movdat(rinc,jdate,idate)
437  sdat(1) = idate(2)
438  sdat(2) = idate(3)
439  sdat(3) = idate(1)
440  ihrst = idate(5)
441  print*,'new forecast hours for restrt run= ',ifhr
442  print*,'new start yr mo day hr min =',sdat(3),sdat(1) &
443  ,sdat(2),ihrst,imin
444  END IF
445 
446  varname='imp_physics'
447  call nemsio_getheadvar(ffile,trim(varname),imp_physics,iret)
448  if (iret /= 0) then
449  if(me==0)print*,varname, &
450  " not found in file-Assigned 99 for Zhao"
451  imp_physics=99 !set GFS mp physics to 99 for Zhao scheme
452  end if
453 
454  if(me==0)print*,'MP_PHYSICS= ',imp_physics
455 
456  varname='sf_surface_physi'
457  call nemsio_getheadvar(ffile,trim(varname),imp_physics,iret)
458  if (iret /= 0) then
459  if(me==0)print*,varname, &
460  " not found in file-Assigned 2 for NOAH"
461  isf_surface_physics=2 !set GFS LSM physics to 2 for NOAH
462  end if
463 
464  if(me==0)print*,'SF_SURFACE_PHYSICS= ',isf_surface_physics
465 
466 ! read bucket
467  varname='fhzero'
468  call nemsio_getheadvar(ffile,trim(varname),fhzero,iret)
469  if (iret /= 0) then
470  if(me==0)print*,varname, &
471  " not found in file-Assign 6 or 12 hours precip bucket"
472  tprec = 6.
473  if(ifhr>240)tprec=12.
474  tclod = tprec
475  trdlw = tprec
476  trdsw = tprec
477  tsrfc = tprec
478  tmaxmin = tprec
479  td3d = tprec
480  else
481  tprec=float(fhzero)
482  tclod = tprec
483  trdlw = tprec
484  trdsw = tprec
485  tsrfc = tprec
486  tmaxmin = tprec
487  td3d = tprec
488  end if
489 
490 
491 ! read meta data to see if precip has zero bucket
492 ! VarName='lprecip_accu'
493 ! call nemsio_getheadvar(ffile,trim(VarName),lprecip_accu,iret)
494 ! if (iret /= 0) then
495 ! if(me==0)print*,VarName, &
496 ! " not found in file-Assign non-zero precip bucket"
497 ! lprecip_accu='no'
498 ! end if
499 ! if(lprecip_accu=='yes')tprec=float(ifhr)
500  print*,'tprec, tclod, trdlw = ',tprec,tclod,trdlw
501 
502 
503 ! Initializes constants for Ferrier microphysics
504  if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then
505  CALL microinit(imp_physics)
506  end if
507 
508 ! GFS does not need DT to compute accumulated fields, set it to one
509 ! VarName='DT'
510  dt = 1
511 
512  hbm2 = 1.0
513 
514 
515 ! start reading nemsio sigma files using parallel read
516  fldsize = (jend-jsta+1)*im
517  allocate(tmp(fldsize*nrec))
518  print*,'allocate tmp successfully'
519  tmp = 0.
520  call nemsio_denseread(nfile,1,im,jsta,jend,tmp,iret=iret)
521  if(iret /=0 ) then
522  print*,"fail to read sigma file using mpi io read, stopping"
523  stop
524  end if
525 !
526 ! covert from reduced grid to full grid
527 !
528  if(reduce_grid)then
529  print*,'performing reduced grid'
530  jtem = jend-jsta+1
531  allocate (kmsk(im,jtem))
532  kmsk = 0
533  do recn=1,nrec
534  fldst = (recn-1)*fldsize
535  do j=jsta,jend
536  js = fldst + (j-jsta)*im
537  do i=1,im
538  buf(i,j) = tmp(i+js)
539  enddo
540  enddo
541  call gg2rg(im,jtem,numi(jsta),buf(1,jsta))
542  call uninterpred(2,kmsk,numi(jsta),im,jtem,buf(1,jsta),tmp(fldst+1))
543  enddo
544  deallocate(kmsk)
545  end if
546 
547 ! Terrain height * G using nemsio
548  varname='hgt'
549  vcoordname = 'sfc'
550  l = 1
551  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
552  if(recn /=0 ) then
553  fldst = (recn-1)*fldsize
554 !$omp parallel do private(i,j,js)
555  do j=jsta,jend
556  js = fldst + (j-jsta)*im
557  do i=1,im
558  fis(i,j) = tmp(i+js)
559  enddo
560  enddo
561  else
562  if(me == 0) print*,'fail to read ', varname,vcoordname,l
563  endif
564 
565 ! call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
566 ! ,l,nrec,fldsize,spval,tmp &
567 ! ,recname,reclevtyp,reclev,VarName,VcoordName &
568 ! ,fis)
569 
570 !$omp parallel do private(i,j)
571  do j=jsta,jend
572  do i=1,im
573  if (fis(i,j) /= spval) then
574  zint(i,j,lp1) = fis(i,j)
575  fis(i,j) = fis(i,j) * grav
576  endif
577  enddo
578  enddo
579  if(debugprint) print*,'sample ',varname,' = ',fis(isa,jsa)
580 
581  vcoordname = 'sfc' ! surface fileds
582  l = 1
583 
584 ! Surface pressure using nemsio
585  varname='pres'
586  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
587  ,l,nrec,fldsize,spval,tmp &
588  ,recname,reclevtyp,reclev,varname,vcoordname &
589  ,pint(1,jsta_2l,lp1))
590 
591  if(debugprint)print*,'sample surface pressure = ',pint(isa,jsa,lp1)
592 
593 !
594 ! vertical loop for Layer 3d fields
595 ! --------------------------------
596  vcoordname = 'mid layer'
597 
598  do l=1,lm
599  ll = lm-l+1
600 ! model level T
601  !if (me == 0) print*,'start retrieving GFS T using nemsio'
602  varname='tmp'
603  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
604  if(recn /= 0) then
605  fldst = (recn-1)*fldsize
606 !$omp parallel do private(i,j,js)
607  do j=jsta,jend
608  js = fldst + (j-jsta)*im
609  do i=1,im
610  t(i,j,ll) = tmp(i+js)
611  enddo
612  enddo
613  else
614  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
615  stop
616  endif
617 
618  if(debugprint)print*,'sample ',ll,varname,' = ',ll,t(isa,jsa,ll)
619 
620 ! model level q
621  varname='spfh'
622  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
623  if(recn /= 0) then
624  fldst = (recn-1)*fldsize
625 !$omp parallel do private(i,j,js)
626  do j=jsta,jend
627  js = fldst + (j-jsta)*im
628  do i=1,im
629  q(i,j,ll) = tmp(i+js)
630  enddo
631  enddo
632  else
633  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
634  stop
635  endif
636 
637  if(debugprint)print*,'sample ',ll,varname,' = ',ll,q(isa,jsa,ll)
638 
639 ! model level u
640  varname='ugrd'
641  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
642  if(recn /= 0) then
643  fldst = (recn-1)*fldsize
644 !$omp parallel do private(i,j,js)
645  do j=jsta,jend
646  js = fldst + (j-jsta)*im
647  do i=1,im
648  uh(i,j,ll) = tmp(i+js)
649  enddo
650  enddo
651  else
652  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
653  stop
654  endif
655 
656  if(debugprint)print*,'sample ',ll,varname,' = ',ll,uh(isa,jsa,ll)
657 
658 ! model level v
659  varname='vgrd'
660  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
661  if(recn /= 0) then
662  fldst = (recn-1)*fldsize
663 !$omp parallel do private(i,j,js)
664  do j=jsta,jend
665  js = fldst + (j-jsta)*im
666  do i=1,im
667  vh(i,j,ll) = tmp(i+js)
668  enddo
669  enddo
670  else
671  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
672  stop
673  endif
674 
675  if(debugprint)print*,'sample ',ll,varname,' = ',ll,vh(isa,jsa,ll)
676 
677 ! model level pressure
678 ! if (.not. hyb_sigp) then
679 ! VarName='pres'
680 ! call getrecn(recname,reclevtyp,reclev,nrec,varname,VcoordName,l,recn)
681 ! if(recn /= 0) then
682 ! fldst = (recn-1)*fldsize
683 !!$omp parallel do private(i,j,js)
684 ! do j=jsta,jend
685 ! js = fldst + (j-jsta)*im
686 ! do i=1,im
687 ! pmid(i,j,ll) = tmp(i+js)
688 ! enddo
689 ! enddo
690 ! else
691 ! recn_pres=-9999
692 ! if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
693 ! 'will derive pressure using ak bk later'
694 ! stop
695 ! endif
696 
697 ! if(debugprint)print*,'sample ',ll,VarName,' = ',ll,pmid(isa,jsa,ll)
698 ! end if
699 ! GFS is on A grid and does not need PMIDV
700 
701 ! dp
702  varname='dpres'
703  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
704  if(recn /= 0) then
705  fldst = (recn-1)*fldsize
706 !$omp parallel do private(i,j,js)
707  do j=jsta,jend
708  js = fldst + (j-jsta)*im
709  do i=1,im
710  dpres(i,j,ll) = tmp(i+js)
711  enddo
712  enddo
713  else
714  recn_dpres = -9999
715  if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
716  'will derive pressure using ak bk later'
717  endif
718 ! ozone mixing ratio
719  varname='o3mr'
720  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
721  if(recn /= 0) then
722  fldst = (recn-1)*fldsize
723 !$omp parallel do private(i,j,js)
724  do j=jsta,jend
725  js = fldst + (j-jsta)*im
726  do i=1,im
727  o3(i,j,ll) = tmp(i+js)
728  enddo
729  enddo
730  else
731  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
732  stop
733  endif
734 
735 
736  if(debugprint)print*,'sample ',ll,varname,' = ',ll,o3(isa,jsa,ll)
737 
738 ! cloud water and ice mixing ratio for zhao scheme
739 ! need to look up old eta post to derive cloud water/ice from cwm
740 ! Zhao scheme does not produce suspended rain and snow
741 
742 !!$omp parallel do private(i,j)
743  do j = jsta, jend
744  do i=1,im
745  qqw(i,j,ll) = 0.
746  qqr(i,j,ll) = 0.
747  qqs(i,j,ll) = 0.
748  qqi(i,j,ll) = 0.
749  enddo
750  enddo
751 
752  if(imp_physics==99 .or. imp_physics==98)then
753  varname='clwmr'
754  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
755  if(recn /= 0) then
756  fldst = (recn-1)*fldsize
757 !$omp parallel do private(i,j,js)
758  do j=jsta,jend
759  js = fldst + (j-jsta)*im
760  do i=1,im
761  cwm(i,j,ll) = tmp(i+js)
762  enddo
763  enddo
764  else
765  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
766  stop
767  endif
768 
769  if(debugprint)print*,'sample ',ll,varname,' = ',ll,cwm(isa,jsa,ll)
770 
771 !$omp parallel do private(i,j)
772  do j=jsta,jend
773  do i=1,im
774  if(t(i,j,ll) < (tfrz-15.) )then ! dividing cloud water from ice
775  qqi(i,j,ll) = cwm(i,j,ll)
776  else
777  qqw(i,j,ll) = cwm(i,j,ll)
778  end if
779  end do
780  end do
781  else if(imp_physics==11 .or. imp_physics==8)then ! GFDL or Thompson MP scheme
782  varname='clwmr'
783  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
784  if(recn /= 0) then
785  fldst = (recn-1)*fldsize
786 !$omp parallel do private(i,j,js)
787  do j=jsta,jend
788  js = fldst + (j-jsta)*im
789  do i=1,im
790  qqw(i,j,ll) = tmp(i+js)
791  enddo
792  enddo
793  else
794  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
795  stop
796  endif
797  if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqw(isa,jsa,ll)
798 
799  varname='icmr'
800  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
801  if(recn /= 0) then
802  fldst = (recn-1)*fldsize
803 !$omp parallel do private(i,j,js)
804  do j=jsta,jend
805  js = fldst + (j-jsta)*im
806  do i=1,im
807  qqi(i,j,ll) = tmp(i+js)
808  enddo
809  enddo
810  else
811  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
812  stop
813  endif
814  if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqi(isa,jsa,ll)
815 
816  varname='rwmr'
817  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
818  if(recn /= 0) then
819  fldst = (recn-1)*fldsize
820 !$omp parallel do private(i,j,js)
821  do j=jsta,jend
822  js = fldst + (j-jsta)*im
823  do i=1,im
824  qqr(i,j,ll) = tmp(i+js)
825  enddo
826  enddo
827  else
828  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
829  stop
830  endif
831  if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqr(isa,jsa,ll)
832 
833  varname ='snmr'
834  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
835  if(recn /= 0) then
836  fldst = (recn-1)*fldsize
837 !$omp parallel do private(i,j,js)
838  do j=jsta,jend
839  js = fldst + (j-jsta)*im
840  do i=1,im
841  qqs(i,j,ll) = tmp(i+js)
842  enddo
843  enddo
844  else
845  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
846  stop
847  endif
848  if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqs(isa,jsa,ll)
849 
850  varname ='grle'
851  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
852  if(recn /= 0) then
853  fldst = (recn-1)*fldsize
854 !$omp parallel do private(i,j,js)
855  do j=jsta,jend
856  js = fldst + (j-jsta)*im
857  do i=1,im
858  qqg(i,j,ll) = tmp(i+js)
859  enddo
860  enddo
861  else
862  print*,'fail to read ', varname,' at lev ',ll, 'stopping'
863  stop
864  endif
865  if(debugprint)print*,'sample ',ll,varname,' = ',ll,qqg(isa,jsa,ll)
866 !define cwm
867  do j=jsta,jend
868  do i=1,im
869  cwm(i,j,ll)=qqg(i,j,ll)+qqs(i,j,ll)+qqr(i,j,ll)+qqi(i,j,ll)+qqw(i,j,ll)
870  enddo
871  enddo
872 
873  end if ! end of reading MP species for diff MP options
874 
875 ! if (iret /= 0)print*,'Error scattering array';stop
876 
877 ! pressure vertical velocity
878  if(trim(modelname_nemsio)=='FV3GFS')then
879  recn_vvel = 0 ! do not derive omega
880  varname='dzdt'
881  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
882  if(recn /= 0) then
883  fldst = (recn-1)*fldsize
884 !$omp parallel do private(i,j,js)
885  do j=jsta,jend
886  js = fldst + (j-jsta)*im
887  do i=1,im
888  wh(i,j,ll) = tmp(i+js)
889  enddo
890  enddo
891  if(debugprint)print*,'sample l ',varname,' = ',ll,wh(isa,jsa,ll)
892  else
893  if(me==0)print*,'fail to read ', varname,' at lev ',ll
894  end if
895  else
896  varname='vvel'
897  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
898  if(recn /= 0) then
899  fldst = (recn-1)*fldsize
900 !$omp parallel do private(i,j,js)
901  do j=jsta,jend
902  js = fldst + (j-jsta)*im
903  do i=1,im
904  omga(i,j,ll) = tmp(i+js)
905  enddo
906  enddo
907  if(debugprint)print*,'sample l ',varname,' = ',ll,omga(isa,jsa,ll)
908  else
909  recn_vvel = -9999
910  if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
911  'will derive omega later'
912  endif
913  end if
914 
915 ! pressure vertical velocity
916  varname='delz'
917  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
918  if(recn /= 0) then
919  fldst = (recn-1)*fldsize
920 ! make sure delz is positive.
921 !$omp parallel do private(i,j,js)
922  do j=jsta,jend
923  js = fldst + (j-jsta)*im
924  do i=1,im
925  zint(i,j,ll)=zint(i,j,ll+1)+abs(tmp(i+js))
926  if(recn_dpres /= -9999)pmid(i,j,ll)=rgas*dpres(i,j,ll)* &
927  t(i,j,ll)*(q(i,j,ll)*fv+1.0)/grav/abs(tmp(i+js))
928  enddo
929  enddo
930  if(debugprint)print*,'sample l ',varname,' = ',ll, &
931  zint(isa,jsa,ll)
932  if(trim(modelname_nemsio)=='FV3GFS' .and. &
933  recn_dpres /= -9999)then
934  do j=jsta,jend
935  js = fldst + (j-jsta)*im
936  do i=1,im
937  omga(i,j,ll)=(-1.)*wh(i,j,ll)*dpres(i,j,ll)/abs(tmp(i+js))
938  end do
939  end do
940  if(debugprint)print*,'sample l omga for FV3',ll, &
941  omga(isa,jsa,ll)
942  end if
943  else
944  recn_delz = -9999
945  if(me==0)print*,'fail to read ', varname,' at lev ',ll, &
946  'will derive height later'
947  endif
948 
949 ! cloud fraction
950  varname='cld_amt'
951  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
952  if(recn /= 0) then
953  fldst = (recn-1)*fldsize
954 !$omp parallel do private(i,j,js)
955  do j=jsta,jend
956  js = fldst + (j-jsta)*im
957  do i=1,im
958  cfr(i,j,ll)=tmp(i+js)
959  enddo
960  enddo
961 ! if(debugprint)print*,'sample l ',VarName,' = ',ll, &
962 ! cfr(isa,jsa,ll)
963  endif
964 
965  if(imp_physics == 99)then
966  allocate(p2d(im,lm),t2d(im,lm),q2d(im,lm),cw2d(im,lm), &
967  qs2d(im,lm),cfr2d(im,lm))
968  do j=jsta,jend
969  do k=1,lm
970  do i=1,im
971  p2d(i,k) = pmid(i,j,ll)*0.01
972  t2d(i,k) = t(i,j,ll)
973  q2d(i,k) = q(i,j,ll)
974  cw2d(i,k) = cwm(i,j,ll)
975  es = min(fpvsnew(t(i,j,ll)),pmid(i,j,ll))
976  qs2d(i,k) = eps*es/(pmid(i,j,ll)+epsm1*es)!saturation q for GFS
977  enddo
978  enddo
979  call progcld1 &
980 !...................................
981 ! --- inputs:
982  ( p2d,t2d,q2d,qs2d,cw2d,im,lm,0, &
983 ! --- outputs:
984  cfr2d &
985  )
986 !$omp parallel do private(i,k)
987  do k=1,lm
988  do i=1,im
989  cfr(i,j,k) = cfr2d(i,k)
990  enddo
991  end do
992  end do
993  deallocate(p2d,t2d,q2d,qs2d,cw2d,cfr2d)
994  end if
995 
996 ! With SHOC NEMS/GSM does output TKE now
997  varname='tke'
998  recn = 0
999  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
1000  if(recn /=0 ) then
1001  fldst = (recn-1)*fldsize
1002 !$omp parallel do private(i,j,js)
1003  do j=jsta,jend
1004  js = fldst + (j-jsta)*im
1005  do i=1,im
1006  q2(i,j,ll) = tmp(i+js)
1007  enddo
1008  enddo
1009  else
1010  if(me==0)print*,'fail to read ', varname,' at lev ',ll
1011 !$omp parallel do private(i,j)
1012  do j=jsta,jend
1013  do i=1,im
1014  q2(i,j,ll) = spval
1015  end do
1016  end do
1017  endif
1018  if(debugprint)print*,'sample l ',varname,' = ',ll,q2(isa,jsa,ll)
1019 
1020 ! Read model derived radar ref.
1021  varname='ref3D'
1022  recn = 0
1023  call getrecn(recname,reclevtyp,reclev,nrec,varname,vcoordname,l,recn)
1024  if(recn /=0 ) then
1025 !$omp parallel do private(i,j,js)
1026  do j=jsta,jend
1027  js = fldst + (j-jsta)*im
1028  do i=1,im
1029  ref_10cm(i,j,ll) = tmp(i+js)
1030  enddo
1031  enddo
1032  else
1033 !$omp parallel do private(i,j)
1034  do j=jsta,jend
1035  do i=1,im
1036  ref_10cm(i,j,ll) = spval
1037  end do
1038  end do
1039  if(me==0)print*,'fail to read ', varname,' at lev ',ll
1040  endif
1041  if(debugprint)print*,'sample l ',varname,' = ',ll,ref_10cm(isa,jsa,ll)
1042 
1043 
1044  end do ! do loop for l
1045 
1046 ! construct interface pressure from model top (which is zero) and dp from top down PDTOP
1047 ! pdtop = spval
1048 ! pt = 0.
1049 ! pd = spval ! GFS does not output PD
1050  pt=ak5(lp1)
1051 
1052  ii = im/2
1053  jj = (jsta+jend)/2
1054 
1055 !!!!! COMPUTE Z, GFS integrates Z on mid-layer instead
1056 !!! use GFS contants to see if height becomes more aggreable to GFS pressure grib file
1057 
1058  if (recn_dpres == -9999) then
1059  do l=lm,1,-1
1060 !$omp parallel do private(i,j)
1061  do j=jsta,jend
1062  do i=1,im
1063  pint(i,j,l) = ak5(lm+2-l) + bk5(lm+2-l)*pint(i,j,lp1)
1064  if(recn_delz == -9999)pmid(i,j,l) = 0.5*(pint(i,j,l)+ &
1065  pint(i,j,l+1)) ! for now - Moorthi
1066  enddo
1067  enddo
1068  if (me == 0) print*,'sample pint,pmid' ,ii,jj,l,pint(ii,jj,l),pmid(ii,jj,l)
1069  enddo
1070  else
1071 ! compute pint using dpres from bot up
1072 ! do l=lm,1,-1
1073 ! do j=jsta,jend
1074 ! do i=1,im
1075 ! pint(i,j,l) = pint(i,j,l+1) - dpres(i,j,l)
1076 ! enddo
1077 ! enddo
1078 ! if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
1079 ! ,pint(ii,jj,l),pmid(ii,jj,l)
1080 ! end do
1081 
1082 !Feb 6 2018: per Jun, change to compute pint from top down
1083  do j=jsta,jend
1084  do i=1,im
1085  pint(i,j,1)=ak5(lp1)
1086  end do
1087  end do
1088 
1089  do l=2,lp1
1090  do j=jsta,jend
1091  do i=1,im
1092  pint(i,j,l) = pint(i,j,l-1) + dpres(i,j,l-1)
1093  enddo
1094  enddo
1095  if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l &
1096  ,pint(ii,jj,l)
1097  end do
1098  endif
1099 
1100 !
1101 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1102 ! Compute omega
1103 !sk05132016
1104 
1105  if (recn_vvel == -9999) then !only compute omega if it's not in model output
1106  allocate(ps2d(im,jsta_2l:jend_2u), psx2d(im,jsta_2l:jend_2u), &
1107  psy2d(im,jsta_2l:jend_2u))
1108  allocate(div3d(im,jsta:jend,lm))
1109 
1110 !$omp parallel do private(i,j)
1111  do j=jsta,jend
1112  do i=1,im
1113  ps2d(i,j) = log(pint(i,j,lm+1))
1114  enddo
1115  enddo
1116  call calgradps(ps2d,psx2d,psy2d)
1117 
1118  call caldiv(uh, vh, div3d)
1119 
1120 !----------------------------------------------------------------------
1121  allocate (vcrd(lm+1,2), d2d(im,lm), u2d(im,lm), v2d(im,lm), &
1122  pi2d(im,lm+1), pm2d(im,lm), omga2d(im,lm))
1123  omga2d=spval
1124  idvc = 2
1125  idsl = 2
1126  nvcoord = 2
1127  do l=1,lm+1
1128  vcrd(l,1) = vcoord4(l,1,1)
1129  vcrd(l,2) = vcoord4(l,2,1)
1130  enddo
1131 
1132 ! jtem = jm / 18 + 1
1133  jtem = jm / 20 + 1
1134  do j=jsta,jend
1135  npass = npass2
1136 ! if (j > jm-jtem+1 .or. j < jtem) npass = npass3
1137  if (j > jm-jtem+1) then
1138  npass = npass + nint(0.5*(j-jm+jtem-1))
1139  elseif (j < jtem) then
1140  npass = npass + nint(0.5*(jtem-j))
1141  endif
1142 ! npass = 0
1143 !$omp parallel do private(i,l,ll)
1144  do l=1,lm
1145  ll = lm-l+1
1146  do i=1,im
1147  u2d(i,l) = uh(i,j,ll) !flip u & v for calling modstuff
1148  v2d(i,l) = vh(i,j,ll)
1149  d2d(i,l) = div3d(i,j,ll)
1150  end do
1151  end do
1152 
1153  call modstuff2(im, im, lm, idvc, idsl, nvcoord, &
1154  vcrd, pint(1,j,lp1), psx2d(1,j), psy2d(1,j), &
1155  d2d, u2d, v2d, pi2d, pm2d, omga2d, me)
1156 ! if (j ==1 .or. j == jm) &
1157 ! write (0,*)' omga2d=',omga2d(1,:),' j=',j
1158 
1159  if (npass <= 0 ) then
1160 !$omp parallel do private(i,l,ll)
1161  do l=1,lm
1162  ll = lm-l+1
1163  do i=1,im
1164  omga(i,j,l) = omga2d(i,ll)
1165 ! pmid(i,j,l) = pm2d(i,ll)
1166 ! pint(i,j,l) = pi2d(i,ll+1)
1167  enddo
1168  enddo
1169  else
1170 !$omp parallel do private(i,l,ll,nn,omg1,omg2)
1171  do l=1,lm
1172  ll = lm-l+1
1173  do i=1,im
1174  omg1(i) = omga2d(i,ll)
1175  enddo
1176  do nn=1,npass
1177  do i=1,im
1178  omg2(i+1) = omg1(i)
1179  enddo
1180  omg2(1) = omg2(im+1)
1181  omg2(im+2) = omg2(2)
1182  do i=2,im+1
1183  omg1(i-1) = third * (omg2(i-1) + omg2(i) + omg2(i+1))
1184  enddo
1185  enddo
1186 
1187  do i=1,im
1188  omga(i,j,l) = omg1(i)
1189  enddo
1190 ! if (j ==1 .or. j == jm) &
1191 ! write (1000+me,*)' omga=',omga(:,j,l),' j=',j,' l=',l
1192 
1193  enddo
1194  endif
1195 
1196 ! Average omega for the last point near the poles - moorthi
1197  if (j ==1 .or. j == jm) then
1198  tx1 = 1.0 / im
1199 ! write(0,*)' j=',j,' omgamax=',maxval(omga(:,j,:)),' omgamin=',minval(omga(:,j,:))
1200 !$omp parallel do private(i,l,tx2)
1201  do l=1,lm
1202  tx2 = 0.0
1203  do i=1,im
1204  tx2 = tx2 + omga(i,j,l)
1205  enddo
1206  tx2 = tx2 * tx1
1207  do i=1,im
1208  omga(i,j,l) = tx2
1209  enddo
1210  enddo
1211  endif
1212 
1213  enddo ! end of j loop
1214 
1215  deallocate (vcrd,d2d,u2d,v2d,pi2d,pm2d,omga2d)
1216  deallocate (ps2d,psx2d,psy2d,div3d)
1217  end if
1218  deallocate (vcoord4)
1219 
1220 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1221 
1222 !
1223  allocate(wrk1(im,jsta:jend),wrk2(im,jsta:jend))
1224  allocate(fi(im,jsta:jend,2))
1225 
1226  do j=jsta,jend
1227  do i=1,im
1228  pd(i,j) = spval ! GFS does not output PD
1229  pint(i,j,1) = pt
1230  end do
1231  end do
1232 
1233  do l=lp1,1,-1
1234  do j=jsta,jend
1235  do i=1,im
1236  alpint(i,j,l)=log(pint(i,j,l))
1237  end do
1238  end do
1239  end do
1240 
1241  if (recn_delz == -9999) then !only compute H if it's not in model output
1242  do j=jsta,jend
1243  do i=1,im
1244  wrk1(i,j) = log(pmid(i,j,lm))
1245  wrk2(i,j) = t(i,j,lm)*(q(i,j,lm)*fv+1.0)
1246  fi(i,j,1) = fis(i,j) &
1247  + wrk2(i,j)*rgas*(alpint(i,j,lp1)-wrk1(i,j))
1248  zmid(i,j,lm) = fi(i,j,1) * gravi
1249  end do
1250  end do
1251 
1252  DO l=lm,2,-1 ! omit computing model top height because it's infinity
1253  ll = l - 1
1254  do j = jsta, jend
1255  do i = 1, im
1256  tvll = t(i,j,ll)*(q(i,j,ll)*fv+1.0)
1257  pmll = log(pmid(i,j,ll))
1258 
1259  fi(i,j,2) = fi(i,j,1) + (0.5*rgas)*(wrk2(i,j)+tvll) &
1260  * (wrk1(i,j)-pmll)
1261  zmid(i,j,ll) = fi(i,j,2) * gravi
1262 !
1263  fact = (alpint(i,j,l)-wrk1(i,j)) / (pmll-wrk1(i,j))
1264  zint(i,j,l) = zmid(i,j,l) +(zmid(i,j,ll)-zmid(i,j,l))*fact
1265  fi(i,j,1) = fi(i,j,2)
1266  wrk1(i,j) = pmll
1267  wrk2(i,j) = tvll
1268  ENDDO
1269  ENDDO
1270 
1271  if (me == 0) print*,'L ZINT= ',l,zint(ii,jj,l), &
1272  'alpint=',alpint(ii,jj,l),'pmid=',log(pmid(ii,jj,l)), &
1273  'pmid(l-1)=',log(pmid(ii,jj,l-1)),'zmd=',zmid(ii,jj,l), &
1274  'zmid(l-1)=',zmid(ii,jj,l-1)
1275  ENDDO
1276  deallocate(wrk1,wrk2,fi)
1277  else
1278  do l=lm,1,-1
1279  do j=jsta,jend
1280  do i=1,im
1281  zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* &
1282  (log(pmid(i,j,l))-alpint(i,j,l+1))/ &
1283  (alpint(i,j,l)-alpint(i,j,l+1))
1284  end do
1285  end do
1286  end do
1287  end if
1288 
1289 ! do j=jsta,jend
1290 ! do i=1,im
1291 ! pd(i,j) = spval ! GFS does not output PD
1292 ! pint(i,j,1) = PT
1293 ! alpint(i,j,lp1) = log(pint(i,j,lp1))
1294 ! wrk1(i,j) = log(PMID(I,J,LM))
1295 ! wrk2(i,j) = T(I,J,LM)*(Q(I,J,LM)*fv+1.0)
1296 ! FI(I,J,1) = FIS(I,J) &
1297 ! + wrk2(i,j)*rgas*(ALPINT(I,J,Lp1)-wrk1(i,j))
1298 ! ZMID(I,J,LM) = FI(I,J,1) * gravi
1299 ! end do
1300 ! end do
1301 
1302 ! SECOND, INTEGRATE HEIGHT HYDROSTATICLY, GFS integrate height on mid-layer
1303 
1304 ! DO L=LM,2,-1 ! omit computing model top height because it's infinity
1305 ! ll = l - 1
1306 ! do j = jsta, jend
1307 ! do i = 1, im
1308 ! alpint(i,j,l) = log(pint(i,j,l))
1309 ! tvll = T(I,J,LL)*(Q(I,J,LL)*fv+1.0)
1310 ! pmll = log(PMID(I,J,LL))
1311 
1312 ! FI(I,J,2) = FI(I,J,1) + (0.5*rgas)*(wrk2(i,j)+tvll) &
1313 ! * (wrk1(i,j)-pmll)
1314 ! ZMID(I,J,LL) = FI(I,J,2) * gravi
1315 !
1316 ! FACT = (ALPINT(I,J,L)-wrk1(i,j)) / (pmll-wrk1(i,j))
1317 ! ZINT(I,J,L) = ZMID(I,J,L) +(ZMID(I,J,LL)-ZMID(I,J,L))*FACT
1318 ! FI(I,J,1) = FI(I,J,2)
1319 ! wrk1(i,J) = pmll
1320 ! wrk2(i,j) = tvll
1321 ! ENDDO
1322 ! ENDDO
1323 
1324 ! if (me == 0) print*,'L ZINT= ',l,zint(ii,jj,l), &
1325 ! 'alpint=',ALPINT(ii,jj,l),'pmid=',LOG(PMID(Ii,Jj,L)), &
1326 ! 'pmid(l-1)=',LOG(PMID(Ii,Jj,L-1)),'zmd=',ZMID(Ii,Jj,L), &
1327 ! 'zmid(l-1)=',ZMID(Ii,Jj,L-1)
1328 ! ENDDO
1329 ! deallocate(wrk1,wrk2)
1330 
1331 
1332  print *, 'gocart_on2=',gocart_on
1333  if (gocart_on) then
1334 
1335 ! GFS output dust in nemsio (GOCART)
1336  dustcb=0.0
1337  dustallcb=0.0
1338  do n=1,nbin_du
1339  do l=1,lm
1340 !$omp parallel do private(i,j)
1341  do j=jsta_2l,jend_2u
1342  do i=1,im
1343  dust(i,j,l,n) = spval
1344  enddo
1345  enddo
1346  enddo
1347  enddo
1348 ! DUST = SPVAL
1349  !VarName='du001'
1350  varname='dust1'
1351  vcoordname='mid layer'
1352  do l=1,lm
1353  ll=lm-l+1
1354  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1355  ,l,nrec,fldsize,spval,tmp &
1356  ,recname,reclevtyp,reclev,varname,vcoordname &
1357  ,dust(1:im,jsta_2l:jend_2u,ll,1))
1358 
1359 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,1)
1360  end do ! do loop for l
1361 
1362  !VarName='du002'
1363  varname='dust2'
1364  vcoordname='mid layer'
1365  do l=1,lm
1366  ll=lm-l+1
1367  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1368  ,l,nrec,fldsize,spval,tmp &
1369  ,recname,reclevtyp,reclev,varname,vcoordname &
1370  ,dust(1:im,jsta_2l:jend_2u,ll,2))
1371 
1372  dustcb(1:im,jsta_2l:jend_2u)=dustcb(1:im,jsta_2l:jend_2u)+ &
1373  (dust(1:im,jsta_2l:jend_2u,ll,1)+0.38*dust(1:im,jsta_2l:jend_2u,ll,2))* &
1374  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1375 
1376 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,2)
1377  end do ! do loop for l
1378 
1379  !VarName='du003'
1380  varname='dust3'
1381  vcoordname='mid layer'
1382  do l=1,lm
1383  ll=lm-l+1
1384  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1385  ,l,nrec,fldsize,spval,tmp &
1386  ,recname,reclevtyp,reclev,varname,vcoordname &
1387  ,dust(1:im,jsta_2l:jend_2u,ll,3))
1388 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,3)
1389  end do ! do loop for l
1390 
1391  !VarName='du004'
1392  varname='dust4'
1393  vcoordname='mid layer'
1394  do l=1,lm
1395  ll=lm-l+1
1396  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1397  ,l,nrec,fldsize,spval,tmp &
1398  ,recname,reclevtyp,reclev,varname,vcoordname &
1399  ,dust(1:im,jsta_2l:jend_2u,ll,4))
1400 
1401 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,4)
1402  end do ! do loop for l
1403 
1404  !VarName='du005'
1405  varname='dust5'
1406  vcoordname='mid layer'
1407  do l=1,lm
1408  ll=lm-l+1
1409  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1410  ,l,nrec,fldsize,spval,tmp &
1411  ,recname,reclevtyp,reclev,varname,vcoordname &
1412  ,dust(1:im,jsta_2l:jend_2u,ll,5))
1413 
1414  dustallcb(1:im,jsta_2l:jend_2u)=dustallcb(1:im,jsta_2l:jend_2u)+ &
1415  (dust(1:im,jsta_2l:jend_2u,ll,1)+dust(1:im,jsta_2l:jend_2u,ll,2)+ &
1416  dust(1:im,jsta_2l:jend_2u,ll,3)+0.74*dust(1:im,jsta_2l:jend_2u,ll,4))* &
1417  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1418 
1419 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,dust(isa,jsa,ll,5)
1420  end do ! do loop for l
1421 !
1422 ! GFS output sea salt in nemsio (GOCART)
1423  sscb=0.0
1424  ssallcb=0.0
1425  do n=1,nbin_ss
1426  do l=1,lm
1427 !$omp parallel do private(i,j)
1428  do j=jsta_2l,jend_2u
1429  do i=1,im
1430  salt(i,j,l,n) = spval
1431  enddo
1432  enddo
1433  enddo
1434  enddo
1435 ! SALT = SPVAL
1436  !VarName='ss001'
1437  varname='seas1'
1438  vcoordname='mid layer'
1439  do l=1,lm
1440  ll=lm-l+1
1441  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1442  ,l,nrec,fldsize,spval,tmp &
1443  ,recname,reclevtyp,reclev,varname,vcoordname &
1444  ,salt(1:im,jsta_2l:jend_2u,ll,1))
1445 
1446 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,1)
1447  end do ! do loop for l
1448 
1449  !VarName='ss002'
1450  varname='seas2'
1451  vcoordname='mid layer'
1452  do l=1,lm
1453  ll=lm-l+1
1454  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1455  ,l,nrec,fldsize,spval,tmp &
1456  ,recname,reclevtyp,reclev,varname,vcoordname &
1457  ,salt(1:im,jsta_2l:jend_2u,ll,2))
1458 
1459 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,2)
1460  end do ! do loop for l
1461 
1462  !VarName='ss003'
1463  varname='seas3'
1464  vcoordname='mid layer'
1465  do l=1,lm
1466  ll=lm-l+1
1467  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1468  ,l,nrec,fldsize,spval,tmp &
1469  ,recname,reclevtyp,reclev,varname,vcoordname &
1470  ,salt(1:im,jsta_2l:jend_2u,ll,3))
1471 
1472  sscb(1:im,jsta_2l:jend_2u)=sscb(1:im,jsta_2l:jend_2u)+ &
1473  (salt(1:im,jsta_2l:jend_2u,ll,1)+ &
1474  salt(1:im,jsta_2l:jend_2u,ll,2)+0.83*salt(1:im,jsta_2l:jend_2u,ll,3))* &
1475  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1476 
1477 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,3)
1478  end do ! do loop for l
1479 
1480  !VarName='ss004'
1481  varname='seas4'
1482  vcoordname='mid layer'
1483  do l=1,lm
1484  ll=lm-l+1
1485  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1486  ,l,nrec,fldsize,spval,tmp &
1487  ,recname,reclevtyp,reclev,varname,vcoordname &
1488  ,salt(1:im,jsta_2l:jend_2u,ll,4))
1489 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,4)
1490  end do ! do loop for l
1491 
1492  !VarName='ss005'
1493  varname='seas5'
1494  vcoordname='mid layer'
1495  do l=1,lm
1496  ll=lm-l+1
1497  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1498  ,l,nrec,fldsize,spval,tmp &
1499  ,recname,reclevtyp,reclev,varname,vcoordname &
1500  ,salt(1:im,jsta_2l:jend_2u,ll,5))
1501 
1502  ssallcb(1:im,jsta_2l:jend_2u)=ssallcb(1:im,jsta_2l:jend_2u)+ &
1503  (salt(1:im,jsta_2l:jend_2u,ll,1)+salt(1:im,jsta_2l:jend_2u,ll,2)+ &
1504  salt(1:im,jsta_2l:jend_2u,ll,3)+ &
1505  salt(1:im,jsta_2l:jend_2u,ll,4))* &
1506  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1507 
1508 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,salt(isa,jsa,ll,5)
1509  end do ! do loop for l
1510 
1511 ! GFS output black carbon in nemsio (GOCART)
1512  bccb=0.0
1513  do n=1,nbin_bc
1514  do l=1,lm
1515 !$omp parallel do private(i,j)
1516  do j=jsta_2l,jend_2u
1517  do i=1,im
1518  soot(i,j,l,n) = spval
1519  enddo
1520  enddo
1521  enddo
1522  enddo
1523 ! SOOT = SPVAL
1524  !VarName='bcphobic'
1525  varname='bc1'
1526  vcoordname='mid layer'
1527  do l=1,lm
1528  ll=lm-l+1
1529  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1530  ,l,nrec,fldsize,spval,tmp &
1531  ,recname,reclevtyp,reclev,varname,vcoordname &
1532  ,soot(1:im,jsta_2l:jend_2u,ll,1))
1533 
1534 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,soot(isa,jsa,ll,1)
1535  end do ! do loop for l
1536 
1537  !VarName='bcphilic'
1538  varname='bc2'
1539  vcoordname='mid layer'
1540  do l=1,lm
1541  ll=lm-l+1
1542  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1543  ,l,nrec,fldsize,spval,tmp &
1544  ,recname,reclevtyp,reclev,varname,vcoordname &
1545  ,soot(1:im,jsta_2l:jend_2u,ll,2))
1546 
1547  bccb(1:im,jsta_2l:jend_2u)=bccb(1:im,jsta_2l:jend_2u)+ &
1548  (soot(1:im,jsta_2l:jend_2u,ll,1)+soot(1:im,jsta_2l:jend_2u,ll,2))* &
1549  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1550 
1551 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,soot(isa,jsa,ll,2)
1552  end do ! do loop for l
1553 
1554  occb=0.0
1555 ! GFS output organic carbon in nemsio (GOCART)
1556  do n=1,nbin_oc
1557  do l=1,lm
1558 !$omp parallel do private(i,j)
1559  do j=jsta_2l,jend_2u
1560  do i=1,im
1561  waso(i,j,l,n) = spval
1562  enddo
1563  enddo
1564  enddo
1565  enddo
1566 ! WASO = SPVAL
1567  !VarName='ocphobic'
1568  varname='oc1'
1569  vcoordname='mid layer'
1570  do l=1,lm
1571  ll=lm-l+1
1572  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1573  ,l,nrec,fldsize,spval,tmp &
1574  ,recname,reclevtyp,reclev,varname,vcoordname &
1575  ,waso(1:im,jsta_2l:jend_2u,ll,1))
1576 
1577 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,waso(isa,jsa,ll,1)
1578  end do ! do loop for l
1579 
1580  !VarName='ocphilic'
1581  varname='oc2'
1582  vcoordname='mid layer'
1583  do l=1,lm
1584  ll=lm-l+1
1585  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1586  ,l,nrec,fldsize,spval,tmp &
1587  ,recname,reclevtyp,reclev,varname,vcoordname &
1588  ,waso(1:im,jsta_2l:jend_2u,ll,2))
1589 
1590  occb(1:im,jsta_2l:jend_2u)=occb(1:im,jsta_2l:jend_2u)+ &
1591  (waso(1:im,jsta_2l:jend_2u,ll,1)+waso(1:im,jsta_2l:jend_2u,ll,2)) * &
1592  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1593 
1594 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,waso(isa,jsa,ll,2)
1595  end do ! do loop for l
1596 
1597 ! GFS output sulfate in nemsio (GOCART)
1598  sulfcb=0.0
1599  do n=1,nbin_su
1600  do l=1,lm
1601 !$omp parallel do private(i,j)
1602  do j=jsta_2l,jend_2u
1603  do i=1,im
1604  suso(i,j,l,n) = spval
1605  enddo
1606  enddo
1607  enddo
1608  enddo
1609 ! SUSO = SPVAL
1610  !VarName='so4'
1611  varname='sulf'
1612  vcoordname='mid layer'
1613  do l=1,lm
1614  ll=lm-l+1
1615  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1616  ,l,nrec,fldsize,spval,tmp &
1617  ,recname,reclevtyp,reclev,varname,vcoordname &
1618  ,suso(1:im,jsta_2l:jend_2u,ll,1))
1619 
1620  sulfcb(1:im,jsta_2l:jend_2u)=sulfcb(1:im,jsta_2l:jend_2u)+ &
1621  suso(1:im,jsta_2l:jend_2u,ll,1)* &
1622  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1623 
1624 ! if(debugprint)print*,'sample l ',VarName,' = ',ll,suso(isa,jsa,ll,1)
1625  end do ! do loop for l
1626 
1627 ! GFS output pp25 in nemsio (GOCART)
1628  pp25cb=0.0
1629  do n=1,nbin_su
1630  do l=1,lm
1631 !$omp parallel do private(i,j)
1632  do j=jsta_2l,jend_2u
1633  do i=1,im
1634  pp25(i,j,l,n) = spval
1635  enddo
1636  enddo
1637  enddo
1638  enddo
1639 ! PP25 = SPVAL
1640  !VarName='so4'
1641  varname='pp25'
1642  vcoordname='mid layer'
1643  do l=1,lm
1644  ll=lm-l+1
1645  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1646  ,l,nrec,fldsize,spval,tmp &
1647  ,recname,reclevtyp,reclev,varname,vcoordname &
1648  ,pp25(1:im,jsta_2l:jend_2u,ll,1))
1649  pp25cb(1:im,jsta_2l:jend_2u)=pp25cb(1:im,jsta_2l:jend_2u)+ &
1650  pp25(1:im,jsta_2l:jend_2u,ll,1)* &
1651  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1652 ! if(debugprint)print*,'sample l ',VarName,' =
1653 ! ',ll,suso(isa,jsa,ll,1)
1654  end do ! do loop for l
1655 ! GFS output pp10 in nemsio (GOCART)
1656  pp10cb=0.0
1657  do n=1,nbin_su
1658  do l=1,lm
1659 !$omp parallel do private(i,j)
1660  do j=jsta_2l,jend_2u
1661  do i=1,im
1662  pp10(i,j,l,n) = spval
1663  enddo
1664  enddo
1665  enddo
1666  enddo
1667 ! PP10 = SPVAL
1668  !VarName='so4'
1669  varname='pp10'
1670  vcoordname='mid layer'
1671  do l=1,lm
1672  ll=lm-l+1
1673  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1674  ,l,nrec,fldsize,spval,tmp &
1675  ,recname,reclevtyp,reclev,varname,vcoordname &
1676  ,pp10(1:im,jsta_2l:jend_2u,ll,1))
1677  pp10cb(1:im,jsta_2l:jend_2u)=pp10cb(1:im,jsta_2l:jend_2u)+ &
1678  pp10(1:im,jsta_2l:jend_2u,ll,1)* &
1679  dpres(1:im,jsta_2l:jend_2u,ll)/grav
1680 ! if(debugprint)print*,'sample l ',VarName,' =
1681 ! ',ll,suso(isa,jsa,ll,1)
1682  end do ! do loop for l
1683 
1684 ! -- compute air density RHOMID and remove negative tracer values
1685  do l=1,lm
1686 !$omp parallel do private(i,j,n,tv)
1687  do j=jsta,jend
1688  do i=1,im
1689 
1690  tv = t(i,j,l) * (h1+d608*max(q(i,j,l),qmin))
1691  rhomid(i,j,l) = pmid(i,j,l) / (rd*tv)
1692  do n = 1, nbin_du
1693  IF ( dust(i,j,l,n) < spval) THEN
1694  dust(i,j,l,n) = max(dust(i,j,l,n), 0.0)
1695  ENDIF
1696  enddo
1697  do n = 1, nbin_ss
1698  IF ( salt(i,j,l,n) < spval) THEN
1699  salt(i,j,l,n) = max(salt(i,j,l,n), 0.0)
1700  ENDIF
1701  enddo
1702  do n = 1, nbin_oc
1703  IF ( waso(i,j,l,n) < spval) THEN
1704  waso(i,j,l,n) = max(waso(i,j,l,n), 0.0)
1705  ENDIF
1706  enddo
1707  do n = 1, nbin_bc
1708  IF ( soot(i,j,l,n) < spval) THEN
1709  soot(i,j,l,n) = max(soot(i,j,l,n), 0.0)
1710  ENDIF
1711  enddo
1712  do n = 1, nbin_su
1713  IF ( suso(i,j,l,n) < spval) THEN
1714  suso(i,j,l,n) = max(suso(i,j,l,n), 0.0)
1715  ENDIF
1716  enddo
1717 
1718  end do
1719  end do
1720  end do
1721  l=lm
1722 !$omp parallel do private(i,j)
1723  do j=jsta,jend
1724  do i=1,im
1725  dustcb(i,j) = max(dustcb(i,j), 0.0)
1726  dustallcb(i,j) = max(dustallcb(i,j), 0.0)
1727  sscb(i,j) = max(sscb(i,j), 0.0)
1728  ssallcb(i,j) = max(ssallcb(i,j), 0.0)
1729  bccb(i,j) = max(bccb(i,j), 0.0)
1730  occb(i,j) = max(occb(i,j), 0.0)
1731  sulfcb(i,j) = max(sulfcb(i,j), 0.0)
1732  pp25cb(i,j) = max(pp25cb(i,j), 0.0)
1733  pp10cb(i,j) = max(pp10cb(i,j), 0.0)
1734 ! PM10 concentration
1735  dusmass(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
1736  0.74*dust(i,j,l,4)+salt(i,j,l,1)+salt(i,j,l,2)+salt(i,j,l,3)+ &
1737  salt(i,j,l,4) + soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
1738  waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1)+pp10(i,j,l,1)) &
1739  *rhomid(i,j,l) !ug/m3
1740 ! PM25 dust and seasalt
1741  dustpm(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2))*rhomid(i,j,l) !ug/m3
1742  dustpm10(i,j)=(dust(i,j,l,1)+dust(i,j,l,2)+dust(i,j,l,3)+ &
1743  0.74*dust(i,j,l,4))*rhomid(i,j,l) !ug/m3
1744  sspm(i,j)=(salt(i,j,l,1)+salt(i,j,l,2)+ &
1745  0.83*salt(i,j,l,3))*rhomid(i,j,l) !ug/m3
1746 ! PM25 concentration
1747  dusmass25(i,j)=(dust(i,j,l,1)+0.38*dust(i,j,l,2)+ &
1748  salt(i,j,l,1)+salt(i,j,l,2)+0.83*salt(i,j,l,3) + &
1749  soot(i,j,l,1)+soot(i,j,l,2)+waso(i,j,l,1)+ &
1750  waso(i,j,l,2) +suso(i,j,l,1)+pp25(i,j,l,1))*rhomid(i,j,l) !ug/m3
1751 ! PM10 column
1752  ducmass(i,j)=dustallcb(i,j)+ssallcb(i,j)+bccb(i,j)+ &
1753  occb(i,j)+sulfcb(i,j)+pp25cb(i,j)+pp10cb(i,j)
1754 ! PM25 column
1755  ducmass25(i,j)=dustcb(i,j)+sscb(i,j)+bccb(i,j)+occb(i,j) &
1756  +sulfcb(i,j)+pp25cb(i,j)
1757  end do
1758  end do
1759  endif ! endif for gocart_on
1760 !
1761 ! done with sigma file, close it for now
1762  call nemsio_close(nfile,iret=status)
1763  deallocate(tmp,recname,reclevtyp,reclev)
1764 
1765 ! open flux file
1766 ! has to move it to beginning to read imp_physics
1767 ! call nemsio_open(ffile,trim(fileNameFlux),'read',mpi_comm_comp &
1768 ! ,iret=status)
1769 ! if ( Status /= 0 ) then
1770 ! print*,'error opening ',fileNameFlux, ' Status = ', Status
1771 ! print*,'skip reading of flux file'
1772 ! endif
1773  call nemsio_getfilehead(ffile,iret=status,nrec=nrec)
1774  print*,'nrec for flux file=',nrec
1775  allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
1776  call nemsio_getfilehead(ffile,iret=iret &
1777  ,recname=recname ,reclevtyp=reclevtyp,reclev=reclev)
1778  if(debugprint)then
1779  if (me == 0)then
1780  do i=1,nrec
1781  print *,'recname,reclevtyp,reclev=',trim(recname(i)),' ', &
1782  trim(reclevtyp(i)),reclev(i)
1783  end do
1784  end if
1785  end if
1786 
1787 ! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD
1788  varname='IVEGSRC'
1789  call nemsio_getheadvar(ffile,trim(varname),ivegsrc,iret)
1790  if (iret /= 0) then
1791  print*,varname,' not found in file-use 1 for IGBP as default'
1792  ivegsrc=1
1793  end if
1794  if (me == 0) print*,'IVEGSRC= ',ivegsrc
1795 
1796 ! set novegtype based on vegetation classification
1797  if(ivegsrc==2)then
1798  novegtype=13
1799  else if(ivegsrc==1)then
1800  novegtype=20
1801  else if(ivegsrc==0)then
1802  novegtype=24
1803  end if
1804  if (me == 0) print*,'novegtype= ',novegtype
1805 
1806  varname='CU_PHYSICS'
1807  call nemsio_getheadvar(ffile,trim(varname),icu_physics,iret)
1808  if (iret /= 0) then
1809  print*,varname," not found in file-Assigned 4 for SAS as default"
1810  icu_physics=4
1811  end if
1812  if (me == 0) print*,'CU_PHYSICS= ',icu_physics
1813 
1814  varname='dtp'
1815  call nemsio_getheadvar(ffile,trim(varname),dtp,iret)
1816  if (iret /= 0) then
1817  print*,varname," not found in file-Assigned 225. for dtp as default"
1818  dtp=225.
1819  end if
1820  if (me == 0) print*,'dtp= ',dtp
1821 
1822 ! Chuang: zhour is when GFS empties bucket last so using this
1823 ! to compute buket will result in changing bucket with forecast time.
1824 ! set default bucket for now
1825 
1826 ! call nemsio_getheadvar(ffile,'zhour',zhour,iret=iret)
1827 ! if(iret == 0) then
1828 ! tprec = 1.0*ifhr-zhour
1829 ! tclod = tprec
1830 ! trdlw = tprec
1831 ! trdsw = tprec
1832 ! tsrfc = tprec
1833 ! tmaxmin = tprec
1834 ! td3d = tprec
1835 ! print*,'tprec from flux file header= ',tprec
1836 ! else
1837 ! print*,'Error reading accumulation bucket from flux file', &
1838 ! 'header - will try to read from env variable FHZER'
1839 ! CALL GETENV('FHZER',ENVAR)
1840 ! read(ENVAR, '(I2)')idum
1841 ! tprec = idum*1.0
1842 ! tclod = tprec
1843 ! trdlw = tprec
1844 ! trdsw = tprec
1845 ! tsrfc = tprec
1846 ! tmaxmin = tprec
1847 ! td3d = tprec
1848 ! print*,'TPREC from FHZER= ',tprec
1849 ! end if
1850 
1851 
1852 ! start reading nemsio flux files using parallel read
1853  fldsize = (jend-jsta+1)*im
1854  allocate(tmp(fldsize*nrec))
1855  print*,'allocate tmp successfully'
1856  tmp=0.
1857  call nemsio_denseread(ffile,1,im,jsta,jend,tmp,iret=iret)
1858 ! can't stop because anl does not have flux file
1859 ! if(iret/=0)then
1860 ! print*,"fail to read flux file using mpi io read, stopping"
1861 ! stop
1862 ! end if
1863 
1864  vcoordname='sfc' ! surface fileds
1865  varname='land'
1866  l=1
1867  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1868  ,l,nrec,fldsize,spval,tmp &
1869  ,recname,reclevtyp,reclev,varname,vcoordname,sm)
1870  if(debugprint)print*,'sample ',varname,' =',sm(im/2,(jsta+jend)/2)
1871 
1872 !$omp parallel do private(i,j)
1873  do j=jsta,jend
1874  do i=1,im
1875  if (sm(i,j) /= spval) sm(i,j) = 1.0 - sm(i,j)
1876  enddo
1877  enddo
1878 
1879 ! sea ice mask
1880 
1881  varname='icec'
1882  vcoordname = 'sfc'
1883  l=1
1884  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1885  ,l,nrec,fldsize,spval,tmp &
1886  ,recname,reclevtyp,reclev,varname,vcoordname,sice)
1887 
1888  if(debugprint)print*,'sample ',varname,' = ',sice(isa,jsa)
1889 
1890 ! where(sice /=spval .and. sice >=1.0)sm=0.0 !sea ice has sea
1891 ! mask=0
1892 ! GFS flux files have land points with non-zero sea ice, per Iredell,
1893 ! these
1894 ! points have sea ice changed to zero, i.e., trust land mask more than
1895 ! sea ice
1896 ! where(sm/=spval .and. sm==0.0)sice=0.0 !specify sea ice=0 at land
1897 
1898 !$omp parallel do private(i,j)
1899  do j=jsta,jend
1900  do i=1,im
1901  if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0
1902  enddo
1903  enddo
1904 
1905 
1906 ! PBL height using nemsio
1907  varname='hpbl'
1908  vcoordname = 'sfc'
1909  l=1
1910  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1911  ,l,nrec,fldsize,spval,tmp &
1912  ,recname,reclevtyp,reclev,varname,vcoordname &
1913  ,pblh)
1914  if(debugprint)print*,'sample ',varname,' = ',pblh(isa,jsa)
1915 
1916 ! frictional velocity using nemsio
1917  varname='fricv'
1918 ! VcoordName='sfc'
1919 ! l=1
1920  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1921  ,l,nrec,fldsize,spval,tmp &
1922  ,recname,reclevtyp,reclev,varname,vcoordname &
1923  ,ustar)
1924 ! if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa)
1925 
1926 ! roughness length using getgb
1927  varname='sfcr'
1928 ! VcoordName='sfc'
1929 ! l=1
1930  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1931  ,l,nrec,fldsize,spval,tmp &
1932  ,recname,reclevtyp,reclev,varname,vcoordname &
1933  ,z0)
1934 ! if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa)
1935 
1936 ! sfc exchange coeff
1937  varname='sfexc'
1938 ! VcoordName='sfc'
1939 ! l=1
1940  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1941  ,l,nrec,fldsize,spval,tmp &
1942  ,recname,reclevtyp,reclev,varname,vcoordname &
1943  ,sfcexc)
1944 
1945 ! aerodynamic conductance
1946  varname='acond'
1947 ! VcoordName='sfc'
1948 ! l=1
1949  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1950  ,l,nrec,fldsize,spval,tmp &
1951  ,recname,reclevtyp,reclev,varname,vcoordname &
1952  ,acond)
1953 
1954 ! surface potential T using getgb
1955  varname='tmp'
1956 ! VcoordName='sfc'
1957 ! l=1
1958  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
1959  ,l,nrec,fldsize,spval,tmp &
1960  ,recname,reclevtyp,reclev,varname,vcoordname &
1961  ,ths)
1962 
1963 ! where(ths/=spval)ths=ths*(p1000/pint(:,:,lp1))**CAPA ! convert to THS
1964 
1965 !$omp parallel do private(i,j)
1966  do j=jsta,jend
1967  do i=1,im
1968  if (ths(i,j) /= spval) then
1969 ! write(0,*)' i=',i,' j=',j,' ths=',ths(i,j),' pint=',pint(i,j,lp1)
1970  ths(i,j) = ths(i,j) * (p1000/pint(i,j,lp1))**capa
1971  endif
1972  qs(i,j) = spval ! GFS does not have surface specific humidity
1973  twbs(i,j) = spval ! GFS does not have inst sensible heat flux
1974  qwbs(i,j) = spval ! GFS does not have inst latent heat flux
1975 !assign sst
1976  if (sm(i,j) /= 0.0) then
1977  if (sice(i,j) >= 0.15) then
1978  sst(i,j)=271.4
1979  else
1980  sst(i,j) = ths(i,j) * (pint(i,j,lp1)/p1000)**capa
1981  endif
1982  else
1983  sst(i,j) = spval
1984  endif
1985  enddo
1986  enddo
1987 ! if(debugprint)print*,'sample ',VarName,' = ',ths(isa,jsa)
1988 
1989 
1990 ! GFS does not have time step and physics time step, make up ones since they
1991 ! are not really used anyway
1992 ! NPHS=2.
1993 ! DT=80.
1994  dtq2 = dtp !MEB need to get physics DT
1995  nphs=2.
1996  dt=dtq2/nphs
1997  tsph = 3600./dt !MEB need to get DT
1998 
1999 ! convective precip in m per physics time step using getgb
2000 ! read 6 hour bucket
2001  varname='cpratb_ave'
2002 ! VcoordName='sfc'
2003 ! l=1
2004  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2005  ,l,nrec,fldsize,spval,tmp &
2006  ,recname,reclevtyp,reclev,varname,vcoordname &
2007  ,avgcprate)
2008 ! where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m
2009 !$omp parallel do private(i,j)
2010  do j=jsta,jend
2011  do i=1,im
2012  if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001)
2013 !wm cprate(i,j) = avgcprate(i,j)
2014  enddo
2015  enddo
2016 ! read continuous bucket
2017  varname='cprat_ave'
2018 ! VcoordName='sfc'
2019 ! l=1
2020  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2021  ,l,nrec,fldsize,spval,tmp &
2022  ,recname,reclevtyp,reclev,varname,vcoordname &
2023  ,avgcprate_cont)
2024 !$omp parallel do private(i,j)
2025  do j=jsta,jend
2026  do i=1,im
2027  if (avgcprate_cont(i,j) /= spval) avgcprate_cont(i,j) = &
2028  avgcprate_cont(i,j) * (dtq2*0.001)
2029  enddo
2030  enddo
2031 
2032 ! if(debugprint)print*,'sample ',VarName,' = ',avgcprate(isa,jsa)
2033 
2034 ! print*,'maxval CPRATE: ', maxval(CPRATE)
2035 
2036 ! time averaged bucketed precip rate
2037  varname='prateb_ave'
2038 ! VcoordName='sfc'
2039 ! l=1
2040  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2041  ,l,nrec,fldsize,spval,tmp &
2042  ,recname,reclevtyp,reclev,varname,vcoordname &
2043  ,avgprec)
2044 ! where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
2045 !$omp parallel do private(i,j)
2046  do j=jsta,jend
2047  do i=1,im
2048  if (avgprec(i,j) /= spval) avgprec(i,j) = avgprec(i,j) * (dtq2*0.001)
2049  enddo
2050  enddo
2051 
2052 ! if(debugprint)print*,'sample ',VarName,' = ',avgprec(isa,jsa)
2053 
2054 ! time averaged continuous precip rate in m per physics time step using getgb
2055  varname='prate_ave'
2056 ! VcoordName='sfc'
2057 ! l=1
2058  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2059  ,l,nrec,fldsize,spval,tmp &
2060  ,recname,reclevtyp,reclev,varname,vcoordname &
2061  ,avgprec_cont)
2062 ! where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m
2063 !$omp parallel do private(i,j)
2064  do j=jsta,jend
2065  do i=1,im
2066  if (avgprec_cont(i,j) /= spval) avgprec_cont(i,j) = avgprec_cont(i,j) &
2067  * (dtq2*0.001)
2068  enddo
2069  enddo
2070 
2071 ! if(debugprint)print*,'sample ',VarName,' = ',avgprec(isa,jsa)
2072 
2073 ! precip rate in m per physics time step
2074  varname='tprcp'
2075 ! VcoordName='sfc'
2076 ! l=1
2077  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2078  ,l,nrec,fldsize,spval,tmp &
2079  ,recname,reclevtyp,reclev,varname,vcoordname &
2080  ,prec)
2081 !$omp parallel do private(i,j)
2082 ! unit of prec and cprate in post is supposed to be m per physics time step
2083 ! it will be converted back to kg/m^2/s by multiplying by density in SURFCE
2084  do j=jsta,jend
2085  do i=1,im
2086  if (prec(i,j) /= spval) prec(i,j) = prec(i,j) * (dtq2*0.001) &
2087  * 1000. / dtp
2088  enddo
2089  enddo
2090 
2091 ! convective precip rate in m per physics time step
2092  varname='cnvprcp'
2093 ! VcoordName='sfc'
2094 ! l=1
2095  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2096  ,l,nrec,fldsize,spval,tmp &
2097  ,recname,reclevtyp,reclev,varname,vcoordname &
2098  ,cprate)
2099 !$omp parallel do private(i,j)
2100  do j=jsta,jend
2101  do i=1,im
2102  if (cprate(i,j) /= spval) then
2103  cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) * 1000. / dtp
2104  else
2105  cprate(i,j) = 0.
2106  endif
2107  enddo
2108  enddo
2109  if(debugprint)print*,'sample ',varname,' = ',cprate(isa,jsa)
2110 
2111 ! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f
2112 
2113 
2114 ! inst snow water eqivalent using nemsio
2115  varname='weasd'
2116 ! VcoordName='sfc'
2117 ! l=1
2118  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2119  ,l,nrec,fldsize,spval,tmp &
2120  ,recname,reclevtyp,reclev,varname,vcoordname &
2121  ,sno)
2122 ! mask water areas
2123 !$omp parallel do private(i,j)
2124  do j=jsta,jend
2125  do i=1,im
2126  if (sm(i,j) == 1.0 .and. sice(i,j)==0.) sno(i,j) = spval
2127  enddo
2128  enddo
2129 ! if(debugprint)print*,'sample ',VarName,' = ',sno(isa,jsa)
2130 
2131 ! ave snow cover
2132  varname='snowc_ave'
2133 ! VcoordName='sfc'
2134 ! l=1
2135  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2136  ,l,nrec,fldsize,spval,tmp &
2137  ,recname,reclevtyp,reclev,varname,vcoordname &
2138  ,snoavg)
2139 ! snow cover is multipled by 100 in SURFCE before writing it out
2140  do j=jsta,jend
2141  do i=1,im
2142  if (sm(i,j)==1.0 .and. sice(i,j)==0.) snoavg(i,j)=spval
2143  if(snoavg(i,j)/=spval)snoavg(i,j)=snoavg(i,j)/100.
2144  end do
2145  end do
2146 
2147 ! snow depth in mm using nemsio
2148  varname='snod'
2149 ! VcoordName='sfc'
2150 ! l=1
2151  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2152  ,l,nrec,fldsize,spval,tmp &
2153  ,recname,reclevtyp,reclev,varname,vcoordname &
2154  ,si)
2155 ! where(si /= spval)si=si*1000. ! convert to mm
2156 !$omp parallel do private(i,j)
2157  do j=jsta,jend
2158  do i=1,im
2159  if (sm(i,j)==1.0 .and. sice(i,j)==0.) si(i,j)=spval
2160  if (si(i,j) /= spval) si(i,j) = si(i,j) * 1000.0
2161  cldefi(i,j) = spval ! GFS does not have convective cloud efficiency
2162  lspa(i,j) = spval ! GFS does not have similated precip
2163  th10(i,j) = spval ! GFS does not have 10 m theta
2164  th10(i,j) = spval ! GFS does not have 10 m theta
2165  q10(i,j) = spval ! GFS does not have 10 m humidity
2166  albase(i,j) = spval ! GFS does not have snow free albedo
2167  enddo
2168  enddo
2169 ! if(debugprint)print*,'sample ',VarName,' = ',si(isa,jsa)
2170 
2171 
2172 ! 2m T using nemsio
2173  varname='tmp'
2174  vcoordname='2 m above gnd'
2175 ! l=1
2176  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2177  ,l,nrec,fldsize,spval,tmp &
2178  ,recname,reclevtyp,reclev,varname,vcoordname &
2179  ,tshltr)
2180 ! if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa)
2181 
2182 ! GFS does not have 2m pres, estimate it, also convert t to theta
2183  Do j=jsta,jend
2184  Do i=1,im
2185  pshltr(i,j)=pint(i,j,lm+1)*exp(-0.068283/tshltr(i,j))
2186  tshltr(i,j)= tshltr(i,j)*(p1000/pshltr(i,j))**capa ! convert to theta
2187 ! if (j == jm/2 .and. mod(i,50) == 0)
2188 ! + print*,'sample 2m T and P after scatter= '
2189 ! + ,i,j,tshltr(i,j),pshltr(i,j)
2190  end do
2191  end do
2192 
2193 ! 2m specific humidity using nemsio
2194  varname='spfh'
2195  vcoordname='2 m above gnd'
2196 ! l=1
2197  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2198  ,l,nrec,fldsize,spval,tmp &
2199  ,recname,reclevtyp,reclev,varname,vcoordname &
2200  ,qshltr)
2201 ! if(debugprint)print*,'sample ',VarName,' = ',qshltr(isa,jsa)
2202 
2203 ! mid day avg albedo in fraction using nemsio
2204  varname='albdo_ave'
2205  vcoordname='sfc'
2206 ! l=1
2207  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2208  ,l,nrec,fldsize,spval,tmp &
2209  ,recname,reclevtyp,reclev,varname,vcoordname &
2210  ,avgalbedo)
2211 ! where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction
2212 !$omp parallel do private(i,j)
2213  do j=jsta,jend
2214  do i=1,im
2215  if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01
2216  enddo
2217  enddo
2218 ! if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa)
2219 
2220 ! time averaged column cloud fractionusing nemsio
2221  varname='tcdc_ave'
2222  vcoordname='atmos col'
2223 ! l=1
2224  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2225  ,l,nrec,fldsize,spval,tmp &
2226  ,recname,reclevtyp,reclev,varname,vcoordname &
2227  ,avgtcdc)
2228 ! where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction
2229 !$omp parallel do private(i,j)
2230  do j=jsta,jend
2231  do i=1,im
2232  if (avgtcdc(i,j) /= spval) avgtcdc(i,j) = avgtcdc(i,j) * 0.01
2233  enddo
2234  enddo
2235 ! if(debugprint)print*,'sample ',VarName,' = ',avgtcdc(isa,jsa)
2236 
2237 ! GFS probably does not use zenith angle
2238 !$omp parallel do private(i,j)
2239  do j=jsta_2l,jend_2u
2240  do i=1,im
2241  czen(i,j) = spval
2242  czmean(i,j) = spval
2243  enddo
2244  enddo
2245 
2246 ! maximum snow albedo in fraction using nemsio
2247  varname='mxsalb'
2248  vcoordname='sfc'
2249 ! l=1
2250  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2251  ,l,nrec,fldsize,spval,tmp &
2252  ,recname,reclevtyp,reclev,varname,vcoordname &
2253  ,mxsnal)
2254 ! where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction
2255 !$omp parallel do private(i,j)
2256  do j=jsta,jend
2257  do i=1,im
2258  if (mxsnal(i,j) /= spval) mxsnal(i,j) = mxsnal(i,j) * 0.01
2259  enddo
2260  enddo
2261 ! if(debugprint)print*,'sample ',VarName,' = ',mxsnal(isa,jsa)
2262 
2263 !$omp parallel do private(i,j)
2264  do j=jsta_2l,jend_2u
2265  do i=1,im
2266  radot(i,j) = spval ! GFS does not have inst surface outgoing longwave
2267  enddo
2268  enddo
2269 
2270 ! GFS probably does not use sigt4, set it to sig*t^4
2271 !$omp parallel do private(i,j,tlmh)
2272  Do j=jsta,jend
2273  Do i=1,im
2274  tlmh = t(i,j,lm) * t(i,j,lm)
2275  sigt4(i,j) = 5.67e-8 * tlmh * tlmh
2276  End do
2277  End do
2278 
2279 ! TG is not used, skip it for now
2280 
2281 ! GFS does not have inst cloud fraction for high, middle, and low cloud
2282 !$omp parallel do private(i,j)
2283  do j=jsta_2l,jend_2u
2284  do i=1,im
2285  cfrach(i,j) = spval
2286  cfracl(i,j) = spval
2287  cfracm(i,j) = spval
2288  enddo
2289  enddo
2290 
2291 ! ave high cloud fraction using nemsio
2292  varname='tcdc_ave'
2293  vcoordname='high cld lay'
2294  l=1
2295  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2296  ,l,nrec,fldsize,spval,tmp &
2297  ,recname,reclevtyp,reclev,varname,vcoordname &
2298  ,avgcfrach)
2299 ! where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction
2300 !$omp parallel do private(i,j)
2301  do j=jsta,jend
2302  do i=1,im
2303  if (avgcfrach(i,j) /= spval) avgcfrach(i,j) = avgcfrach(i,j) * 0.01
2304  enddo
2305  enddo
2306 ! if(debugprint)print*,'sample ',VarName,' = ',avgcfrach(isa,jsa)
2307 
2308 ! ave low cloud fraction using nemsio
2309  varname='tcdc_ave'
2310  vcoordname='low cld lay'
2311  l=1
2312  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2313  ,l,nrec,fldsize,spval,tmp &
2314  ,recname,reclevtyp,reclev,varname,vcoordname &
2315  ,avgcfracl)
2316 ! where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction
2317 !$omp parallel do private(i,j)
2318  do j=jsta,jend
2319  do i=1,im
2320  if (avgcfracl(i,j) /= spval) avgcfracl(i,j) = avgcfracl(i,j) * 0.01
2321  enddo
2322  enddo
2323 ! if(debugprint)print*,'sample ',VarName,' = ',avgcfracl(isa,jsa)
2324 
2325 ! ave middle cloud fraction using nemsio
2326  varname='tcdc_ave'
2327  vcoordname='mid cld lay'
2328  l=1
2329  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2330  ,l,nrec,fldsize,spval,tmp &
2331  ,recname,reclevtyp,reclev,varname,vcoordname &
2332  ,avgcfracm)
2333 ! where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction
2334 !$omp parallel do private(i,j)
2335  do j=jsta,jend
2336  do i=1,im
2337  if (avgcfracm(i,j) /= spval) avgcfracm(i,j) = avgcfracm(i,j) * 0.01
2338  enddo
2339  enddo
2340 ! if(debugprint)print*,'sample ',VarName,' = ',avgcfracm(isa,jsa)
2341 
2342 ! inst convective cloud fraction using nemsio
2343  varname='tcdc'
2344  vcoordname='convect-cld laye'
2345  l=1
2346  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2347  ,l,nrec,fldsize,spval,tmp &
2348  ,recname,reclevtyp,reclev,varname,vcoordname &
2349  ,cnvcfr)
2350 ! where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction
2351 !$omp parallel do private(i,j)
2352  do j=jsta,jend
2353  do i=1,im
2354  if (cnvcfr(i,j) /= spval) cnvcfr(i,j)= cnvcfr(i,j) * 0.01
2355  enddo
2356  enddo
2357 ! if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(isa,jsa)
2358 
2359 ! slope type using nemsio
2360  varname='sltyp'
2361  vcoordname='sfc'
2362  l=1
2363  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2364  ,l,nrec,fldsize,spval,tmp &
2365  ,recname,reclevtyp,reclev,varname,vcoordname &
2366  ,buf)
2367 ! where(buf /= spval)islope=nint(buf)
2368 !$omp parallel do private(i,j)
2369  do j = jsta_2l, jend_2u
2370  do i=1,im
2371  if (buf(i,j) < spval) then
2372  islope(i,j) = nint(buf(i,j))
2373  else
2374  islope(i,j) = 0
2375  endif
2376  enddo
2377  enddo
2378 ! if(debugprint)print*,'sample ',VarName,' = ',islope(isa,jsa)
2379 
2380 ! plant canopy sfc wtr in m using nemsio
2381  varname='cnwat'
2382  vcoordname='sfc'
2383  l=1
2384  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2385  ,l,nrec,fldsize,spval,tmp &
2386  ,recname,reclevtyp,reclev,varname,vcoordname &
2387  ,cmc)
2388 ! where(cmc /= spval)cmc=cmc/1000. ! convert from kg*m^2 to m
2389 !$omp parallel do private(i,j)
2390  do j=jsta,jend
2391  do i=1,im
2392  if (cmc(i,j) /= spval) cmc(i,j) = cmc(i,j) * 0.001
2393  if (sm(i,j) /= 0.0) cmc(i,j) = spval
2394  enddo
2395  enddo
2396 ! if(debugprint)print*,'sample ',VarName,' = ',cmc(isa,jsa)
2397 
2398 !$omp parallel do private(i,j)
2399  do j=jsta_2l,jend_2u
2400  do i=1,im
2401  grnflx(i,j) = spval ! GFS does not have inst ground heat flux
2402  enddo
2403  enddo
2404 
2405 ! frozen precip fraction using nemsio
2406  varname='cpofp'
2407  vcoordname='sfc'
2408  l=1
2409  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2410  ,l,nrec,fldsize,spval,tmp &
2411  ,recname,reclevtyp,reclev,varname,vcoordname &
2412  ,sr)
2413 !$omp parallel do private(i,j)
2414  do j=jsta,jend
2415  do i=1,im
2416  if(sr(i,j) /= spval) then
2417 !set range within (0,1)
2418  sr(i,j)=min(1.,max(0.,sr(i,j)))
2419  endif
2420  enddo
2421  enddo
2422 
2423 ! sea ice skin temperature
2424  varname='ti'
2425  vcoordname='sfc'
2426  l=1
2427  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2428  ,l,nrec,fldsize,spval,tmp &
2429  ,recname,reclevtyp,reclev,varname,vcoordname &
2430  ,ti)
2431 !$omp parallel do private(i,j)
2432  do j=jsta,jend
2433  do i=1,im
2434  if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval
2435  enddo
2436  enddo
2437 
2438 ! vegetation fraction in fraction. using nemsio
2439  varname='veg'
2440  vcoordname='sfc'
2441  l=1
2442  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2443  ,l,nrec,fldsize,spval,tmp &
2444  ,recname,reclevtyp,reclev,varname,vcoordname &
2445  ,vegfrc)
2446 !$omp parallel do private(i,j)
2447  do j=jsta,jend
2448  do i=1,im
2449  if (vegfrc(i,j) /= spval) then
2450  vegfrc(i,j) = vegfrc(i,j) * 0.01
2451  else
2452  vegfrc(i,j) = 0.0
2453  endif
2454  enddo
2455  enddo
2456 ! mask water areas
2457 !$omp parallel do private(i,j)
2458  do j=jsta,jend
2459  do i=1,im
2460  if (sm(i,j) /= 0.0) vegfrc(i,j) = spval
2461  enddo
2462  enddo
2463 ! if(debugprint)print*,'sample ',VarName,' = ',vegfrc(isa,jsa)
2464 
2465 ! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam
2466 
2467  sldpth(1) = 0.10
2468  sldpth(2) = 0.3
2469  sldpth(3) = 0.6
2470  sldpth(4) = 1.0
2471 
2472 ! liquid volumetric soil mpisture in fraction using nemsio
2473  varname='soill'
2474  vcoordname='0-10 cm down'
2475  l=1
2476  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2477  ,l,nrec,fldsize,spval,tmp &
2478  ,recname,reclevtyp,reclev,varname,vcoordname &
2479  ,sh2o(1,jsta_2l,1))
2480 ! mask water areas
2481 !$omp parallel do private(i,j)
2482  do j=jsta,jend
2483  do i=1,im
2484  if (sm(i,j) /= 0.0) sh2o(i,j,1) = spval
2485  enddo
2486  enddo
2487 ! if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,1)
2488 
2489  varname='soill'
2490  vcoordname='10-40 cm down'
2491  l=1
2492  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2493  ,l,nrec,fldsize,spval,tmp &
2494  ,recname,reclevtyp,reclev,varname,vcoordname &
2495  ,sh2o(1,jsta_2l,2))
2496 ! mask water areas
2497 !$omp parallel do private(i,j)
2498  do j=jsta,jend
2499  do i=1,im
2500  if (sm(i,j) /= 0.0) sh2o(i,j,2) = spval
2501  enddo
2502  enddo
2503 ! if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,2)
2504 
2505  varname='soill'
2506  vcoordname='40-100 cm down'
2507  l=1
2508  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2509  ,l,nrec,fldsize,spval,tmp &
2510  ,recname,reclevtyp,reclev,varname,vcoordname &
2511  ,sh2o(1,jsta_2l,3))
2512 ! mask water areas
2513 !$omp parallel do private(i,j)
2514  do j=jsta,jend
2515  do i=1,im
2516  if (sm(i,j) /= 0.0) sh2o(i,j,3) = spval
2517  enddo
2518  enddo
2519 ! if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,3)
2520 
2521  varname='soill'
2522  vcoordname='100-200 cm down'
2523  l=1
2524  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2525  ,l,nrec,fldsize,spval,tmp &
2526  ,recname,reclevtyp,reclev,varname,vcoordname &
2527  ,sh2o(1,jsta_2l,4))
2528 ! mask water areas
2529 !$omp parallel do private(i,j)
2530  do j=jsta,jend
2531  do i=1,im
2532  if (sm(i,j) /= 0.0) sh2o(i,j,4) = spval
2533  enddo
2534  enddo
2535 ! if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,4)
2536 
2537 ! volumetric soil moisture using nemsio
2538  varname='soilw'
2539  vcoordname='0-10 cm down'
2540  l=1
2541  l=1
2542  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2543  ,l,nrec,fldsize,spval,tmp &
2544  ,recname,reclevtyp,reclev,varname,vcoordname &
2545  ,smc(1,jsta_2l,1))
2546 ! mask water areas
2547 !$omp parallel do private(i,j)
2548  do j=jsta,jend
2549  do i=1,im
2550  if (sm(i,j) /= 0.0) smc(i,j,1) = spval
2551  enddo
2552  enddo
2553 ! if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,1)
2554 
2555  varname='soilw'
2556  vcoordname='10-40 cm down'
2557  l=1
2558  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2559  ,l,nrec,fldsize,spval,tmp &
2560  ,recname,reclevtyp,reclev,varname,vcoordname &
2561  ,smc(1,jsta_2l,2))
2562 ! mask water areas
2563 !$omp parallel do private(i,j)
2564  do j=jsta,jend
2565  do i=1,im
2566  if (sm(i,j) /= 0.0) smc(i,j,2) = spval
2567  enddo
2568  enddo
2569 ! if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,2)
2570 
2571  varname='soilw'
2572  vcoordname='40-100 cm down'
2573  l=1
2574  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2575  ,l,nrec,fldsize,spval,tmp &
2576  ,recname,reclevtyp,reclev,varname,vcoordname &
2577  ,smc(1,jsta_2l,3))
2578 ! mask water areas
2579 !$omp parallel do private(i,j)
2580  do j=jsta,jend
2581  do i=1,im
2582  if (sm(i,j) /= 0.0) smc(i,j,3) = spval
2583  enddo
2584  enddo
2585 ! if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,3)
2586 
2587  varname='soilw'
2588  vcoordname='100-200 cm down'
2589  l=1
2590  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2591  ,l,nrec,fldsize,spval,tmp &
2592  ,recname,reclevtyp,reclev,varname,vcoordname &
2593  ,smc(1,jsta_2l,4))
2594 ! mask water areas
2595 !$omp parallel do private(i,j)
2596  do j=jsta,jend
2597  do i=1,im
2598  if (sm(i,j) /= 0.0) smc(i,j,4) = spval
2599  enddo
2600  enddo
2601 ! if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,4)
2602 
2603 ! soil temperature using nemsio
2604  varname='tmp'
2605  vcoordname='0-10 cm down'
2606  l=1
2607  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2608  ,l,nrec,fldsize,spval,tmp &
2609  ,recname,reclevtyp,reclev,varname,vcoordname &
2610  ,stc(1,jsta_2l,1))
2611 ! mask open water areas, combine with sea ice tmp
2612 !$omp parallel do private(i,j)
2613  do j=jsta,jend
2614  do i=1,im
2615  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval
2616  !if (sm(i,j) /= 0.0) stc(i,j,1) = spval
2617  enddo
2618  enddo
2619 ! if(debugprint)print*,'sample l','stc',' = ',1,stc(isa,jsa,1)
2620 
2621  varname='tmp'
2622  vcoordname='10-40 cm down'
2623  l=1
2624  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2625  ,l,nrec,fldsize,spval,tmp &
2626  ,recname,reclevtyp,reclev,varname,vcoordname &
2627  ,stc(1,jsta_2l,2))
2628 ! mask open water areas, combine with sea ice tmp
2629 !$omp parallel do private(i,j)
2630  do j=jsta,jend
2631  do i=1,im
2632  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval
2633  !if (sm(i,j) /= 0.0) stc(i,j,2) = spval
2634  enddo
2635  enddo
2636 ! if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,2)
2637 
2638  varname='tmp'
2639  vcoordname='40-100 cm down'
2640  l=1
2641  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2642  ,l,nrec,fldsize,spval,tmp &
2643  ,recname,reclevtyp,reclev,varname,vcoordname &
2644  ,stc(1,jsta_2l,3))
2645 ! mask open water areas, combine with sea ice tmp
2646 !$omp parallel do private(i,j)
2647  do j=jsta,jend
2648  do i=1,im
2649  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval
2650  !if (sm(i,j) /= 0.0) stc(i,j,3) = spval
2651  enddo
2652  enddo
2653 ! if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,3)
2654 
2655  varname='tmp'
2656  vcoordname='100-200 cm down'
2657  l=1
2658  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2659  ,l,nrec,fldsize,spval,tmp &
2660  ,recname,reclevtyp,reclev,varname,vcoordname &
2661  ,stc(1,jsta_2l,4))
2662 ! mask open water areas, combine with sea ice tmp
2663 !$omp parallel do private(i,j)
2664  do j=jsta,jend
2665  do i=1,im
2666  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval
2667  !if (sm(i,j) /= 0.0) stc(i,j,4) = spval
2668  enddo
2669  enddo
2670 ! if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,4)
2671 
2672 !$omp parallel do private(i,j)
2673  do j=jsta,jend
2674  do i=1,im
2675  acfrcv(i,j) = spval ! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1
2676  ncfrcv(i,j) = 1.0
2677  acfrst(i,j) = spval ! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1
2678  ncfrst(i,j) = 1.0
2679  ssroff(i,j) = spval ! GFS does not have storm runoff
2680  bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF
2681  rlwin(i,j) = spval ! GFS does not have inst incoming sfc longwave
2682  rlwtoa(i,j) = spval ! GFS does not have inst model top outgoing longwave
2683  enddo
2684  enddo
2685 ! trdlw(i,j) = 6.0
2686  ardlw = 1.0 ! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1
2687 
2688 ! time averaged incoming sfc longwave using nemsio
2689  varname='dlwrf_ave'
2690  vcoordname='sfc'
2691  l=1
2692  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2693  ,l,nrec,fldsize,spval,tmp &
2694  ,recname,reclevtyp,reclev,varname,vcoordname &
2695  ,alwin)
2696 
2697 ! inst incoming sfc longwave using nemsio
2698  varname='dlwrf'
2699  vcoordname='sfc'
2700  l=1
2701  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2702  ,l,nrec,fldsize,spval,tmp &
2703  ,recname,reclevtyp,reclev,varname,vcoordname &
2704  ,rlwin)
2705 
2706 ! time averaged outgoing sfc longwave using gfsio
2707  varname='ulwrf_ave'
2708  vcoordname='sfc'
2709  l=1
2710  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2711  ,l,nrec,fldsize,spval,tmp &
2712  ,recname,reclevtyp,reclev,varname,vcoordname &
2713  ,alwout)
2714 ! inst outgoing sfc longwave using nemsio
2715  varname='ulwrf'
2716  vcoordname='sfc'
2717  l=1
2718  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2719  ,l,nrec,fldsize,spval,tmp &
2720  ,recname,reclevtyp,reclev,varname,vcoordname &
2721  ,radot)
2722 
2723 ! where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing
2724 !$omp parallel do private(i,j)
2725  do j=jsta,jend
2726  do i=1,im
2727  if (alwout(i,j) /= spval) alwout(i,j) = -alwout(i,j)
2728  enddo
2729  enddo
2730 ! if(debugprint)print*,'sample l',VarName,' = ',1,alwout(isa,jsa)
2731 
2732 ! time averaged outgoing model top longwave using gfsio
2733  varname='ulwrf_ave'
2734  vcoordname='nom. top'
2735  l=1
2736  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2737  ,l,nrec,fldsize,spval,tmp &
2738  ,recname,reclevtyp,reclev,varname,vcoordname &
2739  ,alwtoa)
2740 ! if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa)
2741 
2742 !$omp parallel do private(i,j)
2743  do j=jsta_2l,jend_2u
2744  do i=1,im
2745  rswin(i,j) = spval ! GFS does not have inst incoming sfc shortwave
2746  rswinc(i,j) = spval ! GFS does not have inst incoming clear sky sfc shortwave
2747  rswout(i,j) = spval ! GFS does not have inst outgoing sfc shortwave
2748  enddo
2749  enddo
2750 
2751 ! GFS incoming sfc longwave has been averaged, set ARDLW to 1
2752  ardsw=1.0
2753 ! trdsw=6.0
2754 
2755 ! time averaged incoming sfc shortwave using gfsio
2756  varname='dswrf_ave'
2757  vcoordname='sfc'
2758  l=1
2759  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2760  ,l,nrec,fldsize,spval,tmp &
2761  ,recname,reclevtyp,reclev,varname,vcoordname &
2762  ,aswin)
2763 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa)
2764 
2765 ! inst incoming sfc shortwave using nemsio
2766  varname='dswrf'
2767  vcoordname='sfc'
2768  l=1
2769  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2770  ,l,nrec,fldsize,spval,tmp &
2771  ,recname,reclevtyp,reclev,varname,vcoordname &
2772  ,rswin)
2773 
2774 ! time averaged incoming sfc uv-b using getgb
2775  varname='duvb_ave'
2776  vcoordname='sfc'
2777  l=1
2778  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2779  ,l,nrec,fldsize,spval,tmp &
2780  ,recname,reclevtyp,reclev,varname,vcoordname &
2781  ,auvbin)
2782 ! if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa)
2783 
2784 ! time averaged incoming sfc clear sky uv-b using getgb
2785  varname='cduvb_ave'
2786  vcoordname='sfc'
2787  l=1
2788  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2789  ,l,nrec,fldsize,spval,tmp &
2790  ,recname,reclevtyp,reclev,varname,vcoordname &
2791  ,auvbinc)
2792 ! if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa)
2793 
2794 ! time averaged outgoing sfc shortwave using gfsio
2795  varname='uswrf_ave'
2796  vcoordname='sfc'
2797  l=1
2798  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2799  ,l,nrec,fldsize,spval,tmp &
2800  ,recname,reclevtyp,reclev,varname,vcoordname &
2801  ,aswout)
2802 ! where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing
2803 !$omp parallel do private(i,j)
2804  do j=jsta,jend
2805  do i=1,im
2806  if (aswout(i,j) /= spval) aswout(i,j) = -aswout(i,j)
2807  enddo
2808  enddo
2809 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswout(isa,jsa)
2810 
2811 ! inst outgoing sfc shortwave using gfsio
2812  varname='uswrf'
2813  vcoordname='sfc'
2814  l=1
2815  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2816  ,l,nrec,fldsize,spval,tmp &
2817  ,recname,reclevtyp,reclev,varname,vcoordname &
2818  ,rswout)
2819 
2820 ! time averaged model top incoming shortwave
2821  varname='dswrf_ave'
2822  vcoordname='nom. top'
2823  l=1
2824  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2825  ,l,nrec,fldsize,spval,tmp &
2826  ,recname,reclevtyp,reclev,varname,vcoordname &
2827  ,aswintoa)
2828 
2829 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa)
2830 
2831 ! time averaged model top outgoing shortwave
2832  varname='uswrf_ave'
2833  vcoordname='nom. top'
2834  l=1
2835  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2836  ,l,nrec,fldsize,spval,tmp &
2837  ,recname,reclevtyp,reclev,varname,vcoordname &
2838  ,aswtoa)
2839 ! if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(isa,jsa)
2840 
2841 ! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux
2842 ! has reversed sign convention using gfsio
2843  varname='shtfl_ave'
2844  vcoordname='sfc'
2845  l=1
2846  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2847  ,l,nrec,fldsize,spval,tmp &
2848  ,recname,reclevtyp,reclev,varname,vcoordname &
2849  ,sfcshx)
2850 ! where (sfcshx /= spval)sfcshx=-sfcshx
2851 !$omp parallel do private(i,j)
2852  do j=jsta,jend
2853  do i=1,im
2854  if (sfcshx(i,j) /= spval) sfcshx(i,j) = -sfcshx(i,j)
2855  enddo
2856  enddo
2857 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(isa,jsa)
2858 
2859 ! inst surface sensible heat flux
2860  varname='shtfl'
2861  vcoordname='sfc'
2862  l=1
2863  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2864  ,l,nrec,fldsize,spval,tmp &
2865  ,recname,reclevtyp,reclev,varname,vcoordname &
2866  ,twbs)
2867 !$omp parallel do private(i,j)
2868  do j=jsta,jend
2869  do i=1,im
2870  if (twbs(i,j) /= spval) twbs(i,j) = -twbs(i,j)
2871  enddo
2872  enddo
2873 
2874 ! GFS surface flux has been averaged, set ASRFC to 1
2875  asrfc=1.0
2876 ! tsrfc=6.0
2877 
2878 ! time averaged surface latent heat flux, multiplied by -1 because wrf model flux
2879 ! has reversed sign vonvention using gfsio
2880  varname='lhtfl_ave'
2881  vcoordname='sfc'
2882  l=1
2883  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2884  ,l,nrec,fldsize,spval,tmp &
2885  ,recname,reclevtyp,reclev,varname,vcoordname &
2886  ,sfclhx)
2887 ! where (sfclhx /= spval)sfclhx=-sfclhx
2888 !$omp parallel do private(i,j)
2889  do j=jsta,jend
2890  do i=1,im
2891  if (sfclhx(i,j) /= spval) sfclhx(i,j) = -sfclhx(i,j)
2892  enddo
2893  enddo
2894 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(isa,jsa)
2895 
2896 ! inst surface latent heat flux
2897  varname='lhtfl'
2898  vcoordname='sfc'
2899  l=1
2900  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2901  ,l,nrec,fldsize,spval,tmp &
2902  ,recname,reclevtyp,reclev,varname,vcoordname &
2903  ,qwbs)
2904 ! where (sfclhx /= spval)sfclhx=-sfclhx
2905 !$omp parallel do private(i,j)
2906  do j=jsta,jend
2907  do i=1,im
2908  if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j)
2909  enddo
2910  enddo
2911 
2912 ! time averaged ground heat flux using nemsio
2913  varname='gflux_ave'
2914  vcoordname='sfc'
2915  l=1
2916  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2917  ,l,nrec,fldsize,spval,tmp &
2918  ,recname,reclevtyp,reclev,varname,vcoordname &
2919  ,subshx)
2920 ! mask water areas
2921 !$omp parallel do private(i,j)
2922  do j=jsta,jend
2923  do i=1,im
2924  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) subshx(i,j) = spval
2925  enddo
2926  enddo
2927 ! if(debugprint)print*,'sample l',VarName,' = ',1,subshx(isa,jsa)
2928 
2929 ! inst ground heat flux using nemsio
2930  varname='gflux'
2931  vcoordname='sfc'
2932  l=1
2933  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2934  ,l,nrec,fldsize,spval,tmp &
2935  ,recname,reclevtyp,reclev,varname,vcoordname &
2936  ,grnflx)
2937 ! mask water areas
2938 !$omp parallel do private(i,j)
2939  do j=jsta,jend
2940  do i=1,im
2941  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) grnflx(i,j) = spval
2942  enddo
2943  enddo
2944 ! time averaged zonal momentum flux using gfsio
2945  varname='uflx_ave'
2946  vcoordname='sfc'
2947  l=1
2948  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2949  ,l,nrec,fldsize,spval,tmp &
2950  ,recname,reclevtyp,reclev,varname,vcoordname &
2951  ,sfcux)
2952 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa)
2953 
2954 ! time averaged meridional momentum flux using nemsio
2955  varname='vflx_ave'
2956  vcoordname='sfc'
2957  l=1
2958  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2959  ,l,nrec,fldsize,spval,tmp &
2960  ,recname,reclevtyp,reclev,varname,vcoordname &
2961  ,sfcvx)
2962 ! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa)
2963 
2964 !$omp parallel do private(i,j)
2965  do j=jsta_2l,jend_2u
2966  do i=1,im
2967 ! snopcx(i,j) =spval ! GFS does not have snow phase change heat flux
2968  sfcuvx(i,j) = spval ! GFS does not use total momentum flux
2969  enddo
2970  enddo
2971 
2972 ! time averaged zonal gravity wave stress using nemsio
2973  varname='u-gwd_ave'
2974  vcoordname='sfc'
2975  l=1
2976  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2977  ,l,nrec,fldsize,spval,tmp &
2978  ,recname,reclevtyp,reclev,varname,vcoordname &
2979  ,gtaux)
2980 
2981 ! if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa)
2982 
2983 ! time averaged meridional gravity wave stress using getgb
2984  varname='v-gwd_ave'
2985  vcoordname='sfc'
2986  l=1
2987  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2988  ,l,nrec,fldsize,spval,tmp &
2989  ,recname,reclevtyp,reclev,varname,vcoordname &
2990  ,gtauy)
2991 ! if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa)
2992 
2993 ! time averaged accumulated potential evaporation
2994  varname='pevpr_ave'
2995  vcoordname='sfc'
2996  l=1
2997  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
2998  ,l,nrec,fldsize,spval,tmp &
2999  ,recname,reclevtyp,reclev,varname,vcoordname &
3000  ,avgpotevp)
3001 ! mask water areas
3002 !$omp parallel do private(i,j)
3003  do j=jsta,jend
3004  do i=1,im
3005  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) avgpotevp(i,j) = spval
3006  enddo
3007  enddo
3008 ! if(debugprint)print*,'sample l',VarName,' = ',1,potevp(isa,jsa)
3009 
3010 ! inst potential evaporation
3011  varname='pevpr'
3012  vcoordname='sfc'
3013  l=1
3014  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3015  ,l,nrec,fldsize,spval,tmp &
3016  ,recname,reclevtyp,reclev,varname,vcoordname &
3017  ,potevp)
3018 ! mask water areas
3019 !$omp parallel do private(i,j)
3020  do j=jsta,jend
3021  do i=1,im
3022  if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) potevp(i,j) = spval
3023  enddo
3024  enddo
3025 
3026  do l=1,lm
3027 !$omp parallel do private(i,j)
3028  do j=jsta_2l,jend_2u
3029  do i=1,im
3030 ! GFS does not have temperature tendency due to long wave radiation
3031  rlwtt(i,j,l) = spval
3032 ! GFS does not have temperature tendency due to short wave radiation
3033  rswtt(i,j,l) = spval
3034 ! GFS does not have temperature tendency due to latent heating from convection
3035  tcucn(i,j,l) = spval
3036  tcucns(i,j,l) = spval
3037 ! GFS does not have temperature tendency due to latent heating from grid scale
3038  train(i,j,l) = spval
3039  enddo
3040  enddo
3041  enddo
3042 
3043 ! set avrain to 1
3044  avrain=1.0
3045  avcnvc=1.0
3046  theat=6.0 ! just in case GFS decides to output T tendency
3047 
3048 ! 10 m u using nemsio
3049  varname='ugrd'
3050  vcoordname='10 m above gnd'
3051  l=1
3052  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3053  ,l,nrec,fldsize,spval,tmp &
3054  ,recname,reclevtyp,reclev,varname,vcoordname &
3055  ,u10)
3056 
3057  do j=jsta,jend
3058  do i=1,im
3059  u10h(i,j)=u10(i,j)
3060  end do
3061  end do
3062 ! if(debugprint)print*,'sample l',VarName,' = ',1,u10(isa,jsa)
3063 
3064 ! 10 m v using gfsio
3065  varname='vgrd'
3066  vcoordname='10 m above gnd'
3067  l=1
3068  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3069  ,l,nrec,fldsize,spval,tmp &
3070  ,recname,reclevtyp,reclev,varname,vcoordname &
3071  ,v10)
3072 
3073  do j=jsta,jend
3074  do i=1,im
3075  v10h(i,j)=v10(i,j)
3076  end do
3077  end do
3078 ! if(debugprint)print*,'sample l',VarName,' = ',1,v10(isa,jsa)
3079 
3080 ! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon
3081 ! VarName='vgtyp'
3082 !Use for fv3 model output
3083  varname='vtype'
3084  vcoordname='sfc'
3085  l=1
3086  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3087  ,l,nrec,fldsize,spval,tmp &
3088  ,recname,reclevtyp,reclev,varname,vcoordname &
3089  ,buf)
3090 ! where (buf /= spval)
3091 ! ivgtyp=nint(buf)
3092 ! elsewhere
3093 ! ivgtyp=0 !need to feed reasonable value to crtm
3094 ! end where
3095 !$omp parallel do private(i,j)
3096  do j = jsta_2l, jend_2u
3097  do i=1,im
3098  if (buf(i,j) < spval) then
3099  ivgtyp(i,j) = nint(buf(i,j))
3100  else
3101  ivgtyp(i,j) = 0
3102  endif
3103  enddo
3104  enddo
3105 ! if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(isa,jsa)
3106 
3107 ! soil type, it's in GFS surface file, hopefully will merge into gfsio soon
3108  varname='sotyp'
3109  vcoordname='sfc'
3110  l=1
3111  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3112  ,l,nrec,fldsize,spval,tmp &
3113  ,recname,reclevtyp,reclev,varname,vcoordname &
3114  ,buf)
3115 ! where (buf /= spval)
3116 ! isltyp=nint(buf)
3117 ! elsewhere
3118 ! isltyp=0 !need to feed reasonable value to crtm
3119 ! end where
3120 !$omp parallel do private(i,j)
3121  do j = jsta_2l, jend_2u
3122  do i=1,im
3123  if (buf(i,j) < spval) then
3124  isltyp(i,j) = nint(buf(i,j))
3125  else
3126  isltyp(i,j) = 0 !need to feed reasonable value to crtm
3127  endif
3128  enddo
3129  enddo
3130 ! if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(isa,jsa)
3131 
3132 !$omp parallel do private(i,j)
3133  do j=jsta_2l,jend_2u
3134  do i=1,im
3135  smstav(i,j) = spval ! GFS does not have soil moisture availability
3136 ! smstot(i,j) = spval ! GFS does not have total soil moisture
3137  sfcevp(i,j) = spval ! GFS does not have accumulated surface evaporation
3138  acsnow(i,j) = spval ! GFS does not have averaged accumulated snow
3139  acsnom(i,j) = spval ! GFS does not have snow melt
3140 ! sst(i,j) = spval ! GFS does not have sst????
3141  thz0(i,j) = ths(i,j) ! GFS does not have THZ0, use THS to substitute
3142  qz0(i,j) = spval ! GFS does not output humidity at roughness length
3143  uz0(i,j) = spval ! GFS does not output u at roughness length
3144  vz0(i,j) = spval ! GFS does not output humidity at roughness length
3145  enddo
3146  enddo
3147  do l=1,lm
3148 !$omp parallel do private(i,j)
3149  do j=jsta_2l,jend_2u
3150  do i=1,im
3151  el_pbl(i,j,l) = spval ! GFS does not have mixing length
3152  exch_h(i,j,l) = spval ! GFS does not output exchange coefficient
3153  enddo
3154  enddo
3155  enddo
3156 ! if(debugprint)print*,'sample l',VarName,' = ',1,thz0(isa,jsa)
3157 
3158 ! retrieve inst convective cloud top, GFS has cloud top pressure instead of index,
3159 ! will need to modify CLDRAD.f to use pressure directly instead of index
3160  varname='pres'
3161  vcoordname='convect-cld top'
3162  l=1
3163  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3164  ,l,nrec,fldsize,spval,tmp &
3165  ,recname,reclevtyp,reclev,varname,vcoordname &
3166  ,ptop)
3167 ! if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa)
3168 
3169 !$omp parallel do private(i,j)
3170  do j=jsta,jend
3171  do i=1,im
3172  htop(i,j) = spval
3173  if(ptop(i,j) <= 0.0) ptop(i,j) = spval
3174  enddo
3175  enddo
3176  do j=jsta,jend
3177  do i=1,im
3178  if(ptop(i,j) < spval)then
3179  do l=1,lm
3180  if(ptop(i,j) <= pmid(i,j,l))then
3181  htop(i,j) = l
3182 ! if(i==ii .and. j==jj)print*,'sample ptop,pmid pmid-1,pint= ', &
3183 ! ptop(i,j),pmid(i,j,l),pmid(i,j,l-1),pint(i,j,l),htop(i,j)
3184  exit
3185  end if
3186  end do
3187  end if
3188  end do
3189  end do
3190 
3191 ! retrieve inst convective cloud bottom, GFS has cloud top pressure instead of index,
3192 ! will need to modify CLDRAD.f to use pressure directly instead of index
3193  varname='pres'
3194  vcoordname='convect-cld bot'
3195  l=1
3196  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3197  ,l,nrec,fldsize,spval,tmp &
3198  ,recname,reclevtyp,reclev,varname,vcoordname &
3199  ,pbot)
3200 ! if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa)
3201 
3202 !$omp parallel do private(i,j)
3203  do j=jsta,jend
3204  do i=1,im
3205  hbot(i,j) = spval
3206  if(pbot(i,j) <= 0.0) pbot(i,j) = spval
3207  enddo
3208  enddo
3209  do j=jsta,jend
3210  do i=1,im
3211 ! if(.not.lb(i,j))print*,'false bitmask for pbot at '
3212 ! + ,i,j,pbot(i,j)
3213  if(pbot(i,j) < spval)then
3214  do l=lm,1,-1
3215  if(pbot(i,j) >= pmid(i,j,l)) then
3216  hbot(i,j) = l
3217 ! if(i==ii .and. j==jj)print*,'sample pbot,pmid= ', &
3218 ! pbot(i,j),pmid(i,j,l),hbot(i,j)
3219  exit
3220  end if
3221  end do
3222  end if
3223  end do
3224  end do
3225 
3226 ! retrieve time averaged low cloud top pressure using nemsio
3227  varname='pres_ave'
3228  vcoordname='low cld top'
3229  l=1
3230  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3231  ,l,nrec,fldsize,spval,tmp &
3232  ,recname,reclevtyp,reclev,varname,vcoordname &
3233  ,ptopl)
3234 ! if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa)
3235 
3236 ! retrieve time averaged low cloud bottom pressure using nemsio
3237  varname='pres_ave'
3238  vcoordname='low cld bot'
3239  l=1
3240  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3241  ,l,nrec,fldsize,spval,tmp &
3242  ,recname,reclevtyp,reclev,varname,vcoordname &
3243  ,pbotl)
3244 ! if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa)
3245 
3246 ! retrieve time averaged low cloud top temperature using nemsio
3247  varname='tmp_ave'
3248  vcoordname='low cld top'
3249  l=1
3250  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3251  ,l,nrec,fldsize,spval,tmp &
3252  ,recname,reclevtyp,reclev,varname,vcoordname &
3253  ,ttopl)
3254 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa)
3255 
3256 ! retrieve time averaged middle cloud top pressure using nemsio
3257  varname='pres_ave'
3258  vcoordname='mid cld top'
3259  l=1
3260  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3261  ,l,nrec,fldsize,spval,tmp &
3262  ,recname,reclevtyp,reclev,varname,vcoordname &
3263  ,ptopm)
3264 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa)
3265 
3266 ! retrieve time averaged middle cloud bottom pressure using nemsio
3267  varname='pres_ave'
3268  vcoordname='mid cld bot'
3269  l=1
3270  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3271  ,l,nrec,fldsize,spval,tmp &
3272  ,recname,reclevtyp,reclev,varname,vcoordname &
3273  ,pbotm)
3274 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa)
3275 
3276 ! retrieve time averaged middle cloud top temperature using nemsio
3277  varname='tmp_ave'
3278  vcoordname='mid cld top'
3279  l=1
3280  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3281  ,l,nrec,fldsize,spval,tmp &
3282  ,recname,reclevtyp,reclev,varname,vcoordname &
3283  ,ttopm)
3284 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa)
3285 
3286 ! retrieve time averaged high cloud top pressure using nemsio *********
3287  varname='pres_ave'
3288  vcoordname='high cld top'
3289  l=1
3290  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3291  ,l,nrec,fldsize,spval,tmp &
3292  ,recname,reclevtyp,reclev,varname,vcoordname &
3293  ,ptoph)
3294 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa)
3295 
3296 ! retrieve time averaged high cloud bottom pressure using nemsio
3297  varname='pres_ave'
3298  vcoordname='high cld bot'
3299  l=1
3300  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3301  ,l,nrec,fldsize,spval,tmp &
3302  ,recname,reclevtyp,reclev,varname,vcoordname &
3303  ,pboth)
3304 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa)
3305 
3306 ! retrieve time averaged high cloud top temperature using nemsio
3307  varname='tmp_ave'
3308  vcoordname='high cld top'
3309  l=1
3310  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3311  ,l,nrec,fldsize,spval,tmp &
3312  ,recname,reclevtyp,reclev,varname,vcoordname &
3313  ,ttoph)
3314 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa)
3315 
3316 ! retrieve boundary layer cloud cover using nemsio
3317  varname='tcdc_ave'
3318  vcoordname='bndary-layer cld'
3319  l=1
3320  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3321  ,l,nrec,fldsize,spval,tmp &
3322  ,recname,reclevtyp,reclev,varname,vcoordname &
3323  ,pblcfr)
3324 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa)
3325 ! where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction
3326 !$omp parallel do private(i,j)
3327  do j = jsta_2l, jend_2u
3328  do i=1,im
3329  if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01
3330  enddo
3331  enddo
3332 
3333 ! retrieve cloud work function using nemsio
3334  varname='cwork_ave'
3335  vcoordname='atmos col'
3336  l=1
3337  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3338  ,l,nrec,fldsize,spval,tmp &
3339  ,recname,reclevtyp,reclev,varname,vcoordname &
3340  ,cldwork)
3341 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa)
3342 
3343 ! retrieve water runoff using nemsio
3344  varname='watr_acc'
3345  vcoordname='sfc'
3346  l=1
3347  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3348  ,l,nrec,fldsize,spval,tmp &
3349  ,recname,reclevtyp,reclev,varname,vcoordname &
3350  ,runoff)
3351 ! mask water areas
3352 !$omp parallel do private(i,j)
3353  do j=jsta,jend
3354  do i=1,im
3355  if (sm(i,j) /= 0.0) runoff(i,j) = spval
3356  enddo
3357  enddo
3358 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa)
3359 
3360 ! retrieve shelter max temperature using nemsio
3361  varname='tmax_max'
3362  vcoordname='2 m above gnd'
3363  l=1
3364  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3365  ,l,nrec,fldsize,spval,tmp &
3366  ,recname,reclevtyp,reclev,varname,vcoordname &
3367  ,maxtshltr)
3368 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,maxtshltr(isa,jsa)
3369 
3370 ! retrieve shelter min temperature using nemsio
3371  varname='tmin_min'
3372  vcoordname='2 m above gnd'
3373  l=1
3374  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3375  ,l,nrec,fldsize,spval,tmp &
3376  ,recname,reclevtyp,reclev,varname,vcoordname &
3377  ,mintshltr)
3378 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
3379 ! 1,mintshltr(im/2,(jsta+jend)/2)
3380 
3381 !$omp parallel do private(i,j)
3382  do j=jsta_2l,jend_2u
3383  do i=1,im
3384  maxrhshltr(i,j) = spval
3385  minrhshltr(i,j) = spval
3386  enddo
3387  enddo
3388 
3389 ! retrieve ice thickness using nemsio
3390  varname='icetk'
3391  vcoordname='sfc'
3392  l=1
3393  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3394  ,l,nrec,fldsize,spval,tmp &
3395  ,recname,reclevtyp,reclev,varname,vcoordname &
3396  ,dzice)
3397 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa)
3398 
3399 ! retrieve wilting point using nemsio
3400  varname='wilt'
3401  vcoordname='sfc'
3402  l=1
3403  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3404  ,l,nrec,fldsize,spval,tmp &
3405  ,recname,reclevtyp,reclev,varname,vcoordname &
3406  ,smcwlt)
3407 ! mask water areas
3408 !$omp parallel do private(i,j)
3409  do j=jsta,jend
3410  do i=1,im
3411  if (sm(i,j) /= 0.0) smcwlt(i,j) = spval
3412  enddo
3413  enddo
3414 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,smcwlt(isa,jsa)
3415 
3416 ! retrieve sunshine duration using nemsio
3417  varname='sunsd_acc'
3418  vcoordname='sfc'
3419  l=1
3420  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3421  ,l,nrec,fldsize,spval,tmp &
3422  ,recname,reclevtyp,reclev,varname,vcoordname &
3423  ,suntime)
3424 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,suntime(isa,jsa)
3425 
3426 ! retrieve field capacity using nemsio
3427  varname='fldcp'
3428  vcoordname='sfc'
3429  l=1
3430  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3431  ,l,nrec,fldsize,spval,tmp &
3432  ,recname,reclevtyp,reclev,varname,vcoordname &
3433  ,fieldcapa)
3434 ! mask water areas
3435 !$omp parallel do private(i,j)
3436  do j=jsta,jend
3437  do i=1,im
3438  if (sm(i,j) /= 0.0) fieldcapa(i,j) = spval
3439  enddo
3440  enddo
3441 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,fieldcapa(isa,jsa)
3442 
3443 ! retrieve time averaged surface visible beam downward solar flux
3444  varname='vbdsf_ave'
3445  vcoordname='sfc'
3446  l=1
3447  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3448  ,l,nrec,fldsize,spval,tmp &
3449  ,recname,reclevtyp,reclev,varname,vcoordname &
3450  ,avisbeamswin)
3451 
3452 ! retrieve time averaged surface visible diffuse downward solar flux
3453  varname='vddsf_ave'
3454  vcoordname='sfc'
3455  l=1
3456  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3457  ,l,nrec,fldsize,spval,tmp &
3458  ,recname,reclevtyp,reclev,varname,vcoordname &
3459  ,avisdiffswin)
3460 
3461 ! retrieve time averaged surface near IR beam downward solar flux
3462  varname='nbdsf_ave'
3463  vcoordname='sfc'
3464  l=1
3465  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3466  ,l,nrec,fldsize,spval,tmp &
3467  ,recname,reclevtyp,reclev,varname,vcoordname &
3468  ,airbeamswin)
3469 
3470 ! retrieve time averaged surface near IR diffuse downward solar flux
3471  varname='nddsf_ave'
3472  vcoordname='sfc'
3473  l=1
3474  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3475  ,l,nrec,fldsize,spval,tmp &
3476  ,recname,reclevtyp,reclev,varname,vcoordname &
3477  ,airdiffswin)
3478 
3479 ! retrieve time averaged surface clear sky outgoing LW
3480  varname='csulf'
3481  vcoordname='sfc'
3482  l=1
3483  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3484  ,l,nrec,fldsize,spval,tmp &
3485  ,recname,reclevtyp,reclev,varname,vcoordname &
3486  ,alwoutc)
3487 
3488 ! retrieve time averaged TOA clear sky outgoing LW
3489  varname='csulf'
3490  vcoordname='nom. top'
3491  l=1
3492  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3493  ,l,nrec,fldsize,spval,tmp &
3494  ,recname,reclevtyp,reclev,varname,vcoordname &
3495  ,alwtoac)
3496 
3497 ! retrieve time averaged surface clear sky outgoing SW
3498  varname='csusf'
3499  vcoordname='sfc'
3500  l=1
3501  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3502  ,l,nrec,fldsize,spval,tmp &
3503  ,recname,reclevtyp,reclev,varname,vcoordname &
3504  ,aswoutc)
3505 
3506 ! retrieve time averaged TOA clear sky outgoing LW
3507  varname='csusf'
3508  vcoordname='nom. top'
3509  l=1
3510  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3511  ,l,nrec,fldsize,spval,tmp &
3512  ,recname,reclevtyp,reclev,varname,vcoordname &
3513  ,aswtoac)
3514 
3515 ! retrieve time averaged surface clear sky incoming LW
3516  varname='csdlf'
3517  vcoordname='sfc'
3518  l=1
3519  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3520  ,l,nrec,fldsize,spval,tmp &
3521  ,recname,reclevtyp,reclev,varname,vcoordname &
3522  ,alwinc)
3523 
3524 ! retrieve time averaged surface clear sky incoming SW
3525  varname='csdsf'
3526  vcoordname='sfc'
3527  l=1
3528  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3529  ,l,nrec,fldsize,spval,tmp &
3530  ,recname,reclevtyp,reclev,varname,vcoordname &
3531  ,aswinc)
3532 
3533 ! retrieve shelter max specific humidity using nemsio
3534  varname='spfhmax_max'
3535  vcoordname='2 m above gnd'
3536  l=1
3537  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3538  ,l,nrec,fldsize,spval,tmp &
3539  ,recname,reclevtyp,reclev,varname,vcoordname &
3540  ,maxqshltr)
3541 ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',
3542 ! 1,maxtshltr(isa,jsa)
3543 
3544 ! retrieve shelter min temperature using nemsio
3545  varname='spfhmin_min'
3546  vcoordname='2 m above gnd'
3547  l=1
3548  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3549  ,l,nrec,fldsize,spval,tmp &
3550  ,recname,reclevtyp,reclev,varname,vcoordname &
3551  ,minqshltr)
3552 
3553 ! retrieve storm runoff using nemsio
3554  varname='ssrun_acc'
3555  vcoordname='sfc'
3556  l=1
3557  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3558  ,l,nrec,fldsize,spval,tmp &
3559  ,recname,reclevtyp,reclev,varname,vcoordname &
3560  ,ssroff)
3561 ! mask water areas
3562 !$omp parallel do private(i,j)
3563  do j=jsta,jend
3564  do i=1,im
3565  if (sm(i,j) /= 0.0) ssroff(i,j) = spval
3566  enddo
3567  enddo
3568 
3569 ! retrieve direct soil evaporation
3570  varname='evbs_ave'
3571  vcoordname='sfc'
3572  l=1
3573  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3574  ,l,nrec,fldsize,spval,tmp &
3575  ,recname,reclevtyp,reclev,varname,vcoordname &
3576  ,avgedir)
3577 ! mask water areas
3578 !$omp parallel do private(i,j)
3579  do j=jsta,jend
3580  do i=1,im
3581  if (sm(i,j) /= 0.0) avgedir(i,j) = spval
3582  enddo
3583  enddo
3584 
3585 ! retrieve CANOPY WATER EVAP
3586  varname='evcw_ave'
3587  vcoordname='sfc'
3588  l=1
3589  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3590  ,l,nrec,fldsize,spval,tmp &
3591  ,recname,reclevtyp,reclev,varname,vcoordname &
3592  ,avgecan)
3593 ! mask water areas
3594 !$omp parallel do private(i,j)
3595  do j=jsta,jend
3596  do i=1,im
3597  if (sm(i,j) /= 0.0) avgecan(i,j) = spval
3598  enddo
3599  enddo
3600 
3601 ! retrieve PLANT TRANSPIRATION
3602  varname='trans_ave'
3603  vcoordname='sfc'
3604  l=1
3605  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3606  ,l,nrec,fldsize,spval,tmp &
3607  ,recname,reclevtyp,reclev,varname,vcoordname &
3608  ,avgetrans)
3609 ! mask water areas
3610 !$omp parallel do private(i,j)
3611  do j=jsta,jend
3612  do i=1,im
3613  if (sm(i,j) /= 0.0) avgetrans(i,j) = spval
3614  enddo
3615  enddo
3616 
3617 ! retrieve snow sublimation
3618  varname='sbsno_ave'
3619  vcoordname='sfc'
3620  l=1
3621  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3622  ,l,nrec,fldsize,spval,tmp &
3623  ,recname,reclevtyp,reclev,varname,vcoordname &
3624  ,avgesnow)
3625 ! mask water areas
3626 !$omp parallel do private(i,j)
3627  do j=jsta,jend
3628  do i=1,im
3629  if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval
3630  enddo
3631  enddo
3632 
3633 ! retrive total soil moisture
3634  varname='soilm'
3635  vcoordname='0-200 cm down'
3636  l=1
3637  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3638  ,l,nrec,fldsize,spval,tmp &
3639  ,recname,reclevtyp,reclev,varname,vcoordname &
3640  ,smstot)
3641 ! mask water areas
3642 !$omp parallel do private(i,j)
3643  do j=jsta,jend
3644  do i=1,im
3645  if (sm(i,j) /= 0.0) smstot(i,j) = spval
3646  enddo
3647  enddo
3648 
3649 ! retrieve snow phase change heat flux
3650  varname='snohf'
3651  vcoordname='sfc'
3652  l=1
3653  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3654  ,l,nrec,fldsize,spval,tmp &
3655  ,recname,reclevtyp,reclev,varname,vcoordname &
3656  ,snopcx)
3657 ! mask water areas
3658 !$omp parallel do private(i,j)
3659  do j=jsta,jend
3660  do i=1,im
3661  if (sm(i,j) /= 0.0) snopcx(i,j) = spval
3662  enddo
3663  enddo
3664 
3665 ! GFS does not have deep convective cloud top and bottom fields
3666 
3667 !$omp parallel do private(i,j)
3668  do j=jsta,jend
3669  do i=1,im
3670  htopd(i,j) = spval
3671  hbotd(i,j) = spval
3672  htops(i,j) = spval
3673  hbots(i,j) = spval
3674  cuppt(i,j) = spval
3675  enddo
3676  enddo
3677 
3678 ! retrieve dust emission fluxes
3679  do k = 1, nbin_du
3680  if ( k == 1) varname='duem001'
3681  if ( k == 2) varname='duem002'
3682  if ( k == 3) varname='duem003'
3683  if ( k == 4) varname='duem004'
3684  if ( k == 5) varname='duem005'
3685  vcoordname='atmos sfc'
3686  l=1
3687  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3688  ,l,nrec,fldsize,spval,tmp &
3689  ,recname,reclevtyp,reclev,varname,vcoordname&
3690  ,duem(1,jsta_2l,k))
3691 ! if(debugprint)print*,'sample ',VarName,' = ',duem(isa,jsa,k)
3692  enddo
3693 
3694 ! retrieve dust sedimentation fluxes
3695  do k = 1, nbin_du
3696  if ( k == 1) varname='dust1sd'
3697  if ( k == 2) varname='dust2sd'
3698  if ( k == 3) varname='dust3sd'
3699  if ( k == 4) varname='dust4sd'
3700  if ( k == 5) varname='dsut5sd'
3701  vcoordname='atmos sfc'
3702  l=1
3703  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3704  ,l,nrec,fldsize,spval,tmp &
3705  ,recname,reclevtyp,reclev,varname,vcoordname&
3706  ,dusd(1,jsta_2l,k))
3707 ! if(debugprint)print*,'sample ',VarName,' = ',dusd(isa,jsa,k)
3708  enddo
3709 
3710 ! retrieve dust dry deposition fluxes
3711  do k = 1, nbin_du
3712  if ( k == 1) varname='dust1dp'
3713  if ( k == 2) varname='dust2dp'
3714  if ( k == 3) varname='dust3dp'
3715  if ( k == 4) varname='dust4dp'
3716  if ( k == 5) varname='dust5dp'
3717  vcoordname='atmos sfc'
3718  l=1
3719  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3720  ,l,nrec,fldsize,spval,tmp &
3721  ,recname,reclevtyp,reclev,varname,vcoordname&
3722  ,dudp(1,jsta_2l,k))
3723  print *,'dudp,ck=',maxval(dudp(1:im,jsta:jend,k)), &
3724  minval(dudp(1:im,jsta:jend,k))
3725 ! if(debugprint)print*,'sample ',VarName,' = ',dudp(isa,jsa,k)
3726  enddo
3727 
3728 ! retrieve dust wet deposition fluxes
3729  do k = 1, nbin_du
3730  if ( k == 1) varname='dust1wtl'
3731  if ( k == 2) varname='dust2wtl'
3732  if ( k == 3) varname='dust3wtl'
3733  if ( k == 4) varname='dust4wtl'
3734  if ( k == 5) varname='dust5wtl'
3735  vcoordname='atmos sfc'
3736  l=1
3737  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3738  ,l,nrec,fldsize,spval,tmp &
3739  ,recname,reclevtyp,reclev,varname,vcoordname&
3740  ,duwt(1,jsta_2l,k))
3741  enddo
3742 ! retrieve dust scavenging fluxes
3743  do k = 1, nbin_du
3744  if ( k == 1) varname='dust1wtc'
3745  if ( k == 2) varname='dust2wtc'
3746  if ( k == 3) varname='dust3wtc'
3747  if ( k == 4) varname='dust4wtc'
3748  if ( k == 5) varname='dust5wtc'
3749  vcoordname='atmos sfc'
3750  l=1
3751  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3752  ,l,nrec,fldsize,spval,tmp &
3753  ,recname,reclevtyp,reclev,varname,vcoordname&
3754  ,dusv(1,jsta_2l,k))
3755  enddo
3756 
3757 ! retrieve seasalt emission fluxes
3758  do k = 1, nbin_ss
3759  if ( k == 1) varname='ssem001'
3760  if ( k == 2) varname='ssem002'
3761  if ( k == 3) varname='ssem003'
3762  if ( k == 4) varname='ssem004'
3763  if ( k == 5) varname='ssem005'
3764  vcoordname='atmos sfc'
3765  l=1
3766  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3767  ,l,nrec,fldsize,spval,tmp &
3768  ,recname,reclevtyp,reclev,varname,vcoordname&
3769  ,ssem(1,jsta_2l,k))
3770  enddo
3771 
3772 ! retrieve seasalt emission fluxes
3773  do k = 1, nbin_ss
3774  if ( k == 1) varname='seas1sd'
3775  if ( k == 2) varname='seas2sd'
3776  if ( k == 3) varname='seas3sd'
3777  if ( k == 4) varname='seas4sd'
3778  if ( k == 5) varname='seas5sd'
3779  vcoordname='sfc'
3780  l=1
3781  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3782  ,l,nrec,fldsize,spval,tmp &
3783  ,recname,reclevtyp,reclev,varname,vcoordname&
3784  ,sssd(1,jsta_2l,k))
3785  enddo
3786 
3787 
3788 ! retrieve seasalt dry deposition fluxes
3789  do k = 1, nbin_ss
3790  if ( k == 1) varname='seas1dp'
3791  if ( k == 2) varname='seas2dp'
3792  if ( k == 3) varname='seas3dp'
3793  if ( k == 4) varname='seas4dp'
3794  if ( k == 5) varname='seas5dp'
3795  vcoordname='atmos sfc'
3796  l=1
3797  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3798  ,l,nrec,fldsize,spval,tmp &
3799  ,recname,reclevtyp,reclev,varname,vcoordname&
3800  ,ssdp(1,jsta_2l,k))
3801  enddo
3802 
3803 ! retrieve seasalt wet deposition fluxes
3804  do k = 1, nbin_ss
3805  if ( k == 1) varname='seas1wtl'
3806  if ( k == 2) varname='seas2wtl'
3807  if ( k == 3) varname='seas3wtl'
3808  if ( k == 4) varname='seas4wtl'
3809  if ( k == 5) varname='seas5wtl'
3810  vcoordname='atmos sfc'
3811  l=1
3812  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3813  ,l,nrec,fldsize,spval,tmp &
3814  ,recname,reclevtyp,reclev,varname,vcoordname&
3815  ,sswt(1,jsta_2l,k))
3816  enddo
3817 
3818 ! retrieve seasalt scavenging fluxes
3819  do k = 1, nbin_ss
3820  if ( k == 1) varname='seas1wtc'
3821  if ( k == 2) varname='seas1wtc'
3822  if ( k == 3) varname='seas1wtc'
3823  if ( k == 4) varname='seas1wtc'
3824  if ( k == 5) varname='seas1wtc'
3825  vcoordname='atmos sfc'
3826  l=1
3827  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3828  ,l,nrec,fldsize,spval,tmp &
3829  ,recname,reclevtyp,reclev,varname,vcoordname&
3830  ,sssv(1,jsta_2l,k))
3831  enddo
3832 
3833 ! retrieve bc emission fluxes
3834  do k = 1, nbin_bc
3835  if ( k == 1) varname='bceman'
3836  if ( k == 2) varname='bcembb'
3837  vcoordname='atmos sfc'
3838  l=1
3839  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3840  ,l,nrec,fldsize,spval,tmp &
3841  ,recname,reclevtyp,reclev,varname,vcoordname&
3842  ,bcem(1,jsta_2l,k))
3843  enddo
3844 
3845 ! retrieve bc sedimentation fluxes
3846  do k = 1, nbin_bc
3847  if ( k == 1) varname='bc1sd'
3848  if ( k == 2) varname='bc2sd'
3849  vcoordname='atmos sfc'
3850  l=1
3851  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3852  ,l,nrec,fldsize,spval,tmp &
3853  ,recname,reclevtyp,reclev,varname,vcoordname&
3854  ,bcsd(1,jsta_2l,k))
3855  enddo
3856 
3857 ! retrieve bc dry deposition fluxes
3858  do k = 1, nbin_bc
3859  if ( k == 1) varname='bc1dp'
3860  if ( k == 2) varname='bc2dp'
3861  vcoordname='atmos sfc'
3862  l=1
3863  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3864  ,l,nrec,fldsize,spval,tmp &
3865  ,recname,reclevtyp,reclev,varname,vcoordname&
3866  ,bcdp(1,jsta_2l,k))
3867  enddo
3868 
3869 ! retrieve bc large wet deposition fluxes
3870  do k = 1, nbin_bc
3871  if ( k == 1) varname='bc1wtl'
3872  if ( k == 2) varname='bc2wtl'
3873  vcoordname='atmos sfc'
3874  l=1
3875  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3876  ,l,nrec,fldsize,spval,tmp &
3877  ,recname,reclevtyp,reclev,varname,vcoordname&
3878  ,bcwt(1,jsta_2l,k))
3879  enddo
3880 
3881 ! retrieve bc convective wet deposition fluxes
3882  do k = 1, nbin_bc
3883  if ( k == 1) varname='bc1wtc'
3884  if ( k == 2) varname='bc2wtc'
3885  vcoordname='atmos sfc'
3886  l=1
3887  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3888  ,l,nrec,fldsize,spval,tmp &
3889  ,recname,reclevtyp,reclev,varname,vcoordname&
3890  ,bcsv(1,jsta_2l,k))
3891  enddo
3892 
3893 ! retrieve oc emission fluxes
3894  do k = 1, nbin_oc
3895  if ( k == 1) varname='oceman'
3896  if ( k == 2) varname='ocembb'
3897  vcoordname='atmos sfc'
3898  l=1
3899  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3900  ,l,nrec,fldsize,spval,tmp &
3901  ,recname,reclevtyp,reclev,varname,vcoordname&
3902  ,ocem(1,jsta_2l,k))
3903  enddo
3904 
3905 ! retrieve oc sedimentation fluxes
3906  do k = 1, nbin_oc
3907  if ( k == 1) varname='oc1sd'
3908  if ( k == 2) varname='oc2sd'
3909  vcoordname='atmos sfc'
3910  l=1
3911  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3912  ,l,nrec,fldsize,spval,tmp &
3913  ,recname,reclevtyp,reclev,varname,vcoordname&
3914  ,ocsd(1,jsta_2l,k))
3915  enddo
3916 
3917 ! retrieve oc dry deposition fluxes
3918  do k = 1, nbin_oc
3919  if ( k == 1) varname='oc1dp'
3920  if ( k == 2) varname='oc2dp'
3921  vcoordname='atmos sfc'
3922  l=1
3923  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3924  ,l,nrec,fldsize,spval,tmp &
3925  ,recname,reclevtyp,reclev,varname,vcoordname&
3926  ,ocdp(1,jsta_2l,k))
3927  enddo
3928 
3929 ! retrieve oc large wet deposition fluxes
3930  do k = 1, nbin_oc
3931  if ( k == 1) varname='oc1wtl'
3932  if ( k == 2) varname='oc2wtl'
3933  vcoordname='atmos sfc'
3934  l=1
3935  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3936  ,l,nrec,fldsize,spval,tmp &
3937  ,recname,reclevtyp,reclev,varname,vcoordname&
3938  ,ocwt(1,jsta_2l,k))
3939  enddo
3940 
3941 ! retrieve oc convective wet deposition fluxes
3942  do k = 1, nbin_oc
3943  if ( k == 1) varname='oc1wtc'
3944  if ( k == 2) varname='oc2wtc'
3945  vcoordname='atmos sfc'
3946  l=1
3947  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3948  ,l,nrec,fldsize,spval,tmp &
3949  ,recname,reclevtyp,reclev,varname,vcoordname&
3950  ,ocsv(1,jsta_2l,k))
3951  enddo
3952 
3953 ! retrieve MIE AOD
3954  varname='maod'
3955  vcoordname='sfc'
3956  l=1
3957  call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3958  ,l,nrec,fldsize,spval,tmp &
3959  ,recname,reclevtyp,reclev,varname,vcoordname&
3960  ,maod(1,jsta_2l))
3961 
3962 
3963 ! done with flux file, close it for now
3964  call nemsio_close(ffile,iret=status)
3965  deallocate(tmp,recname,reclevtyp,reclev)
3966 
3967 !lzhang
3968 !! retrieve sfc mass concentration
3969 ! VarName='DUSMASS'
3970 ! VcoordName='atmos col'
3971 ! l=1
3972 ! call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3973 ! ,l,nrec,fldsize,spval,tmp &
3974 ! ,recname,reclevtyp,reclev,VarName,VcoordName &
3975 ! ,dusmass)
3976 ! if(debugprint)print*,'sample ',VarName,' = ',dusmass(isa,jsa)
3977 
3978 !lzhang
3979 !! retrieve col mass density
3980 ! VarName='DUCMASS'
3981 ! VcoordName='atmos col'
3982 ! l=1
3983 ! call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3984 ! ,l,nrec,fldsize,spval,tmp &
3985 ! ,recname,reclevtyp,reclev,VarName,VcoordName &
3986 ! ,ducmass)
3987 !! if(debugprint)print*,'sample ',VarName,' = ',ducmass(isa,jsa)
3988 
3989 !lzhang
3990 !! retrieve sfc mass concentration (pm2.5)
3991 ! VarName='DUSMASS25'
3992 ! VcoordName='atmos col'
3993 ! l=1
3994 ! call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
3995 ! ,l,nrec,fldsize,spval,tmp &
3996 ! ,recname,reclevtyp,reclev,VarName,VcoordName &
3997 ! ,dusmass25)
3998 ! if(debugprint)print*,'sample ',VarName,' = ',dusmass25(isa,jsa)
3999 
4000 !lzhang
4001 !! retrieve col mass density (pm2.5)
4002 ! VarName='DUCMASS25'
4003 ! VcoordName='atmos col'
4004 ! l=1
4005 ! call assignnemsiovar(im,jsta,jend,jsta_2l,jend_2u &
4006 ! ,l,nrec,fldsize,spval,tmp &
4007 ! ,recname,reclevtyp,reclev,VarName,VcoordName &
4008 ! ,ducmass25)
4009 ! if(debugprint)print*,'sample ',VarName,' = ',ducmass25(isa,jsa)
4010 
4011 ! pos east
4012  call collect_loc(gdlat,dummy)
4013  if(me == 0)then
4014  latstart = nint(dummy(1,1)*gdsdegr)
4015  latlast = nint(dummy(im,jm)*gdsdegr)
4016  print*,'laststart,latlast B bcast= ',latstart,latlast,'gdsdegr=',gdsdegr,&
4017  'dummy(1,1)=',dummy(1,1),dummy(im,jm),'gdlat=',gdlat(1,1)
4018  end if
4019  call mpi_bcast(latstart,1,mpi_integer,0,mpi_comm_comp,irtn)
4020  call mpi_bcast(latlast,1,mpi_integer,0,mpi_comm_comp,irtn)
4021  write(6,*) 'laststart,latlast,me A calling bcast=',latstart,latlast,me
4022  call collect_loc(gdlon,dummy)
4023  if(me == 0)then
4024  lonstart = nint(dummy(1,1)*gdsdegr)
4025  lonlast = nint(dummy(im,jm)*gdsdegr)
4026  end if
4027  call mpi_bcast(lonstart,1,mpi_integer,0,mpi_comm_comp,irtn)
4028  call mpi_bcast(lonlast, 1,mpi_integer,0,mpi_comm_comp,irtn)
4029 
4030  write(6,*)'lonstart,lonlast A calling bcast=',lonstart,lonlast
4031 !
4032 
4033 ! generate look up table for lifted parcel calculations
4034 
4035  thl = 210.
4036  plq = 70000.
4037  pt_tbl = 10000. ! this is for 100 hPa added by Moorthi
4038 
4039  CALL table(ptbl,ttbl,pt_tbl, &
4040  rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
4041 
4042  CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
4043 
4044 !
4045 !
4046  IF(me == 0)THEN
4047  WRITE(6,*)' SPL (POSTED PRESSURE LEVELS) BELOW: '
4048  WRITE(6,51) (spl(l),l=1,lsm)
4049  50 FORMAT(14(f4.1,1x))
4050  51 FORMAT(8(f8.1,1x))
4051  ENDIF
4052 !
4053 !$omp parallel do private(l)
4054  DO l = 1,lsm
4055  alsl(l) = log(spl(l))
4056  END DO
4057 !
4058 !HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
4059  if(me == 0)then
4060  print*,'writing out igds'
4061  igdout = 110
4062 ! open(igdout,file='griddef.out',form='unformatted'
4063 ! + ,status='unknown')
4064  if(maptype == 1)THEN ! Lambert conformal
4065  WRITE(igdout)3
4066  WRITE(6,*)'igd(1)=',3
4067  WRITE(igdout)im
4068  WRITE(igdout)jm
4069  WRITE(igdout)latstart
4070  WRITE(igdout)lonstart
4071  WRITE(igdout)8
4072  WRITE(igdout)cenlon
4073  WRITE(igdout)dxval
4074  WRITE(igdout)dyval
4075  WRITE(igdout)0
4076  WRITE(igdout)64
4077  WRITE(igdout)truelat2
4078  WRITE(igdout)truelat1
4079  WRITE(igdout)255
4080  ELSE IF(maptype == 2)THEN !Polar stereographic
4081  WRITE(igdout)5
4082  WRITE(igdout)im
4083  WRITE(igdout)jm
4084  WRITE(igdout)latstart
4085  WRITE(igdout)lonstart
4086  WRITE(igdout)8
4087  WRITE(igdout)cenlon
4088  WRITE(igdout)dxval
4089  WRITE(igdout)dyval
4090  WRITE(igdout)0
4091  WRITE(igdout)64
4092  WRITE(igdout)truelat2 !Assume projection at +-90
4093  WRITE(igdout)truelat1
4094  WRITE(igdout)255
4095  ! Note: The calculation of the map scale factor at the standard
4096  ! lat/lon and the PSMAPF
4097  ! Get map factor at 60 degrees (N or S) for PS projection, which will
4098  ! be needed to correctly define the DX and DY values in the GRIB GDS
4099  if (truelat1 < 0.) THEN
4100  lat = -60.
4101  else
4102  lat = 60.
4103  end if
4104 
4105  CALL msfps(lat,truelat1*0.001,psmapf)
4106 
4107  ELSE IF(maptype == 3) THEN !Mercator
4108  WRITE(igdout)1
4109  WRITE(igdout)im
4110  WRITE(igdout)jm
4111  WRITE(igdout)latstart
4112  WRITE(igdout)lonstart
4113  WRITE(igdout)8
4114  WRITE(igdout)latlast
4115  WRITE(igdout)lonlast
4116  WRITE(igdout)truelat1
4117  WRITE(igdout)0
4118  WRITE(igdout)64
4119  WRITE(igdout)dxval
4120  WRITE(igdout)dyval
4121  WRITE(igdout)255
4122  ELSE IF(maptype == 0 .OR. maptype == 203)THEN !A STAGGERED E-GRID
4123  WRITE(igdout)203
4124  WRITE(igdout)im
4125  WRITE(igdout)jm
4126  WRITE(igdout)latstart
4127  WRITE(igdout)lonstart
4128  WRITE(igdout)136
4129  WRITE(igdout)cenlat
4130  WRITE(igdout)cenlon
4131  WRITE(igdout)dxval
4132  WRITE(igdout)dyval
4133  WRITE(igdout)64
4134  WRITE(igdout)0
4135  WRITE(igdout)0
4136  WRITE(igdout)0
4137  END IF
4138  end if
4139 !
4140 !
4141 
4142  RETURN
4143  END
4144  subroutine rg2gg(im,jm,numi,a)
4145 !
4146  implicit none
4147 !
4148  integer,intent(in):: im,jm,numi(jm)
4149  real,intent(inout):: a(im,jm)
4150  integer j,ir,ig
4151  real r,t(im)
4152  do j=1,jm
4153  r = real(numi(j))/real(im)
4154  do ig=1,im
4155  ir = mod(nint((ig-1)*r),numi(j)) + 1
4156  t(ig) = a(ir,j)
4157  enddo
4158  do ig=1,im
4159  a(ig,j) = t(ig)
4160  enddo
4161  enddo
4162  end subroutine rg2gg
4163  subroutine gg2rg(im,jm,numi,a)
4164 !
4165  implicit none
4166 !
4167  integer,intent(in):: im,jm,numi(jm)
4168  real,intent(inout):: a(im,jm)
4169  integer j,ir,ig
4170  real r,t(im)
4171  do j=1,jm
4172  r = real(numi(j))/real(im)
4173  do ir=1,numi(j)
4174  ig = nint((ir-1)/r) + 1
4175  t(ir) = a(ig,j)
4176  enddo
4177  do ir=1,numi(j)
4178  a(ir,j) = t(ir)
4179  enddo
4180  enddo
4181  end subroutine gg2rg
4182 
4183  subroutine uninterpred(iord,kmsk,lonsperlat,lonr,latr,fi,f)
4184 !!
4185  implicit none
4186 !!
4187  integer, intent(in) :: iord, lonr, latr
4188  integer, intent(in) :: kmsk(lonr,latr), lonsperlat(latr)
4189  real, intent(in) :: fi(lonr,latr)
4190  real, intent(out) :: f(lonr,latr)
4191  integer j,lons
4192 !!
4193 !!$omp parallel do private(j,lons)
4194  do j=1,latr
4195  lons = lonsperlat(j)
4196  if(lons /= lonr) then
4197  call intlon(iord,1,lons,lonr,kmsk(1,j),fi(1,j),f(1,j))
4198  else
4199  f(:,j) = fi(:,j)
4200  endif
4201  enddo
4202  end subroutine
4203  subroutine intlon(iord,imsk,m1,m2,k1,f1,f2)
4204  implicit none
4205  integer,intent(in) :: iord,imsk,m1,m2
4206  integer,intent(in) :: k1(m1)
4207  real, intent(in) :: f1(m1)
4208  real, intent(out):: f2(m2)
4209  integer i2,in,il,ir
4210  real r,x1
4211  r = real(m1)/real(m2)
4212  do i2=1,m2
4213  x1 = (i2-1)*r
4214  il = int(x1)+1
4215  ir = mod(il,m1)+1
4216  if(iord == 2 .and. (imsk == 0 .or. k1(il) == k1(ir))) then
4217  f2(i2) = f1(il)*(il-x1) + f1(ir)*(x1-il+1)
4218  else
4219  in = mod(nint(x1),m1) + 1
4220  f2(i2) = f1(in)
4221  endif
4222  enddo
4223  end subroutine intlon
masks
Definition: MASKS_mod.f:1
physcons_post
Definition: physcons.f:1
params_mod
Definition: params.F:1
modstuff2
subroutine modstuff2(im, ix, km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om, me)
modstuff2() computes model coordinate dependent functions.
Definition: GFSPOSTSIG.F:333
soil
Definition: SOIL_mod.f:1
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:341
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27
lookup_mod
Definition: LOOKUP.f:1