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