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