c GrADS routines for ATM562 (Fortran 77 version) c http://www.atmos.albany.edu/facstaff/rfovell/ATM562/grads_routines_augmented.f c this version also carries water vapor 2D field (qv) and base state (qb) c version of 10/9/2018 c ======= OVERVIEW =================================================== c subroutine dumpgrads expects the following to be passed as input c igradscnt - a counter c nx, nz - array dimensions c u, v, w, th, pi - 2D arrays dimensioned (nx, nz) c rhou, tb, pib, ub, qb - 1D base state arrays dimensioned (nz) c cpd - specific heat at constant pressure c c with this input, this routine writes c full and perturbation u, th, and pi; full v; and c perturbation pressure (in mb) c subroutine write_to_grads_ctl writes out the ctl file c EXTRA/OPTIONAL: c to add more fields to the GrADS output: c (1) alter dumpgrads to pass more fields, or compute them in dumpgrads c (2) add calls to write_to_grads_dat in dumpgrads for new fields c (3) update ngradsvars (# of variables written out) in write_to_grads_ctl c (4) add code in write_to_grads_ctl to describe the new field c -------------------------------------------------------------------- c subroutine dumpgrads c -------------------------------------------------------------------- subroutine dumpgrads(igradscnt,u,v,w,th,pi,qv,rhou,tb,pib,ub, 1 qb,nx,nz,cpd) dimension u(nx,nz),v(nx,nz),w(nx,nz),th(nx,nz),pi(nx,nz), 1 qv(nx,nz) dimension rhou(nz),tb(nz),pib(nz),ub(nz),qb(nz) dimension temp(nx,nz) dimension zero(nz) nxm=nx nzm=nz do k=1,nzm zero(k)=0. enddo c some calls pass a non-zero base state array to augment to the 2D field call write_to_grads_dat(u,nxm,nzm,nx,nz,zero,0.,1) ! full U is predicted in model call write_to_grads_dat(u,nxm,nzm,nx,nz,ub,-1.,1) ! pert u made by subtracting ub call write_to_grads_dat(v,nxm,nzm,nx,nz,zero,0.,0) ! v component - at scalar point call write_to_grads_dat(w,nxm,nzm,nx,nz,zero,0.,2) call write_to_grads_dat(th,nxm,nzm,nx,nz,tb,1.,0) ! full pot temp call write_to_grads_dat(th,nxm,nzm,nx,nz,zero,0.,0) ! pert pot temp call write_to_grads_dat(qv,nxm,nzm,nx,nz,qb,1.,0) ! full qv call write_to_grads_dat(qv,nxm,nzm,nx,nz,zero,0.,0) ! pert qv call write_to_grads_dat(pi,nxm,nzm,nx,nz,pib,1.,0) ! full ndim pressure call write_to_grads_dat(pi,nxm,nzm,nx,nz,zero,0.,0) ! pert ndim pressure do i=1,nxm do k=1,nzm temp(i,k)=0.01*rhou(k)*cpd*tb(k)*pi(i,k) enddo enddo call write_to_grads_dat(temp,nxm,nzm,nx,nz,zero,0.,0) ! pert pressure in millibars igradscnt=igradscnt+1 print *,' done writing grads write number ',igradscnt return end c -------------------------------------------------------------------- c subroutine write_to_grads_dat c -------------------------------------------------------------------- subroutine write_to_grads_dat(array,nxx,nzx,nxg,nzg,zmean,zfactor, 1 iavg) dimension array(nxx,nzx),zmean(nzx) dimension dummy(nxg,1,nzg) do i=1,nxg do k=1,nzg dummy(i,1,k)=0. enddo enddo if(iavg.eq.0)then ! no averaging needed do i=1,nxg do k=1,nzg dummy(i,1,k)=array(i,k)+zfactor*zmean(k) enddo enddo else if(iavg.eq.1)then ! average in x-direction do k=1,nzg do i=1,nxg-1 dummy(i,1,k)=0.5*(array(i+1,k)+array(i,k))+zfactor*zmean(k) enddo enddo else if(iavg.eq.2)then ! average in z-direction do i=1,nxg do k=1,nzg-1 dummy(i,1,k)=0.5*(array(i,k+1)+zfactor*zmean(k+1) 1 +array(i, k )+zfactor*zmean( k )) enddo enddo else if(iavg.eq.3)then ! average from strfcn point do i=2,nxg-1 do k=2,nzg-1 top_pair=0.5*(array(i+1,k+1)+array(i,k+1))+zfactor*zmean(k+1) bot_pair=0.5*(array(i+1, k )+array(i, k ))+zfactor*zmean( k ) dummy(i,1,k)=0.5*(top_pair+bot_pair) enddo enddo endif ! iavg do k=2,nzg-1 write(66) ((dummy(i,j,k),i=2,nxg-1),j=1,1) enddo return end c -------------------------------------------------------------------- c subroutine write_to_grads_ctl c -------------------------------------------------------------------- subroutine write_to_grads_ctl(casename,nxg,nzg,dt,dx,dz, 1 igradscnt,ipltint,byteswap) character*4 ctl,dat character*80 casename,ctlfile,datfile close(44) close(66) ngradsvars=11 print *,' writing grads control file' c if(icall.eq.0)then ! first call c create names for ctl, dat files knx=index(casename,' ') write(ctl,'(a4)') '.ctl' write(dat,'(a4)') '.dat' ctlfile=casename(1:knx-1)//ctl datfile=casename(1:knx-1)//dat write(6,*) 'casename is ',casename write(6,*) 'ctlfile is ',ctlfile write(6,*) 'datfile is ',datfile open(44,file=ctlfile,status='unknown') open(66,file=datfile,status='unknown', 1 form='unformatted') c endif ! icall c gradstinc=float(ipltint)*dt/60. ! minutes igradstinc=ifix(gradstinc) ! if not integer number of min, time is counted wrong print *,' nxg= ',nxg,' nzg= ',nzg,' gradstinc ',gradstinc write(44,900) datfile ! grads does not like underscores in ctl, dat file names 900 format('DSET ^',a80) write(44,901) 901 format('TITLE ATM 562 model',/, 1 'OPTIONS sequential',/,'UNDEF -9.99E33') if(byteswap.eq.1) write(44,902) 902 format('OPTIONS byteswapped') c write(44,100) nxg-2,dx/2/1000.,dx/1000. ! dx in km write(44,100) nxg-2,-1*float(nxg-2)*dx/2/1000.+0.5*dx/1000., 1 dx/1000. ! dx in km 100 format('XDEF ',i3,' LINEAR ',2f10.5,/,'YDEF 1 LINEAR 1.0 1.0') write(44,101) nzg-2 101 format('ZDEF ',i3,' levels') do k=2,nzg-1 write(44,102) (float(k)-1.5)*dz/1000. 102 format(f10.3) enddo if(igradstinc.lt.1)then print *,' *** WARNING - since plotting interval less ', 1 'than one minute, GrADS will not report time ', 1 'correctly.' print *,' *** Plot interval recorded as 1 min in ', 1 'control file' igradstinc=1 endif write(44,103) igradscnt,igradstinc 103 format('TDEF ',i5,' LINEAR 00:00Z01JAN2000 ',i3,'mn') write(44,104) ngradsvars 104 format('VARS ',i4) write(44,105) nzg-2 105 format('u ',i3,' 00 horizontal velocity') write(44,106) nzg-2 106 format('upr',i3,' 00 pert horizontal velocity') write(44,199) nzg-2 199 format('v ',i3,' 00 north-south velocity') write(44,107) nzg-2 107 format('w ',i3,' 00 vertical velocity') write(44,108) nzg-2 108 format('th ',i3,' 00 potential temperature') write(44,109) nzg-2 109 format('thpr',i3,' 00 pert potential temperature') write(44,118) nzg-2 118 format('qv ',i3,' 00 vapor mixing ratio') write(44,119) nzg-2 119 format('qvpr',i3,' 00 pert vapor mixing ratio') write(44,110) nzg-2 110 format('pi ',i3,' 00 ndim pressure') write(44,111) nzg-2 111 format('pipr',i3,' 00 pert ndim pressure') write(44,112) nzg-2 112 format('pprmb',i3,' 00 pert pressure in millibars') write(44,999) 999 format('ENDVARS') close(44) return end