program adjbges
c
c adjust bges
c
c take observed pentad precip - analyzed precip
c and use it to correct surface moisture
c
c
c unit 10 = input bges
c unit 11 = input observed precip (same grid) (prev pentad)
c unit 12 = input model precip (same grid) (prev pentad)
c (optional) land mask
c (optional) runoff
c unit 51 = output bges
c unit 52 = output grib file
c
c assume that bges and observed precip are correct
c use idate(4) from bges file to get date of month
c
c v1.3 .. used by R2 2001-2004 (over the counter)
c
c v1.4 5/2004 for NCO operations
c read decoded grib files in real*4 rather than real*8
c to eliminate double precision version of wgrib
c do pentad averaging .. eliminate pentad averging program.
c
c v1.5 12/2012 for wcoss .. some f95 notation Wesley Ebisuzaki
parameter (lonf= 192 )
parameter (latg= 94 )
parameter (lsoil= 2 )
c parameter (lonf=192)
c parameter (latg=94)
c parameter (lsoil=2)
parameter (ijdim= lonf * latg)
parameter (lugi=10,lugi1=11,lugi2=12)
parameter (lugo1=51)
parameter (lugo2=52)
parameter (nhour= 6 ) ! analysis cycle time
parameter (undef=9.999e20) ! wgrib undefined number
parameter (undeflow=9.998e20) ! wgrib undefined number
c rlev1 = 0 .. 10 cm below ground
c rlev2 = 10 .. 200 cm below ground
parameter(rlev1 = 112 * 256 * 256 + 0 * 256 + 10)
parameter(rlev2 = 112 * 256 * 256 + 10 * 256 + 200)
character*32 label
real fhour ! forecast hour
integer idate(4) ! initial date code
real sfctmp(ijdim) ! sfc temp
real soilm(ijdim,lsoil) ! soil moisture
real snow(ijdim) ! sfc temp
real temp(ijdim), tempv(ijdim,lsoil) ! temp array
real obspre(ijdim) ! observed precip
real mdlpre(ijdim) ! model precip
real offfset, bias(ijdim,lsoil) ! correction term in cm
real land(ijdim) ! land mask, 1=land 0=water
real runoff(ijdim) ! runoff
real*4 tmp1(ijdim), tmp2(ijdim), tmp3(ijdim)
real soilcap(lsoil) ! soil moisture capacity for each layer (cm)
integer year, month, day, hour
real factor, rtmp(ijdim)
logical dorun, flag1(ijdim), flag2(ijdim), flag3(ijdim)
data soilcap /100.,1900./ ! layer depth in mm
dorun = .true.
c initialize bias
bias = 0.0
c read observed precip
c read real*4, convert to real*8
read(lugi1,end=100,err=100) tmp1
obspre = tmp1
goto 110
100 continue
c no obs precip
write(*,*) "WARNING!!!!"
write(*,*) "Obs precip read error or end of record"
write(*,*) "ZERO SOIL WETNESS CORRECTION"
do i = 1, ijdim
obspre(i) = undef
enddo
110 continue
call maxmin(obspre,ijdim,rmin,rmax)
write(*,*) 'obspre min/max (mm/s) =', rmin, rmax
c
c read model precip
c
mdlpre = 0.0
land = 0.0
runoff = 0.0
flag1 = .true.
flag2 = .true.
flag3 = .true.
n = 0
120 continue
read(lugi2,end=140,err=130) tmp1
read(lugi2,end=130,err=130) tmp2
read(lugi2,end=130,err=130) tmp3
rtmp = tmp1
call maxmin(rtmp,ijdim,rmin,rmax)
write(*,*) 'n mdlpre min/max=', n, rmin, rmax
do i = 1, ijdim
if (tmp1(i).gt.undeflow) flag1(i) = .false.
if (tmp2(i).gt.undeflow) flag2(i) = .false.
if (tmp3(i).gt.undeflow) flag3(i) = .false.
mdlpre(i) = mdlpre(i) + tmp1(i)
land(i) = land(i) + tmp2(i)
runoff(i) = runoff(i) + tmp3(i)
enddo
n = n + 1
goto 120
130 continue
write(*,*) 'WARNING!!!!'
write(*,*) 'READ ERROR - MODEL precip - precipadj'
write(*,*) 'WARNING!!!!'
n = 0
140 continue
if (n.ne.20.and.n.ne.24) then
write(*,*) "WARNING!!!!"
write(*,*) "Bad model precip file"
write(*,*) "running without soil moisture adjustment"
mdlpre = undef
land = 1.0
runoff = 0.0
n = 0
dorun = .false.
stop 8
else
c Average model parameters
factor = 1.0 / n
do i = 1, ijdim
if (flag1(i)) then
mdlpre(i) = factor*mdlpre(i)
else
mdlpre(i) = undef
endif
if (flag2(i)) then
land(i) = factor*land(i)
else
land(i) = 1.0
endif
if (flag3(i)) then
runoff(i) = factor*runoff(i)
else
runoff(i) = 0.0
land(i) = 0.0
endif
enddo
endif
c fixup values, runoff in mm (per 6 hours)
if (dorun) then
do i = 1, ijdim
if (abs((runoff(i)-undef)/undef).lt.1e-5) runoff(i) = 0.0
enddo
endif
close(lugi1)
close(lugi2)
write(*,*) 'obspre (pnt 1 100):',obspre(1), obspre(100)
write(*,*) 'mdlpre:(pnt i 100)',mdlpre(1), mdlpre(100)
write(*,*) 'runoff is ',dorun
if (dorun) write(*,*) 'runoff:',runoff(1), runoff(100)
call maxmin(obspre,ijdim,rmin,rmax)
write(*,*) 'obspre min/max=', rmin, rmax
call maxmin(mdlpre,ijdim,rmin,rmax)
write(*,*) 'mdlpre min/max=', rmin, rmax
call maxmin(runoff,ijdim,rmin,rmax)
write(*,*) 'runoff min/max=', rmin, rmax
call maxmin(land,ijdim,rmin,rmax)
write(*,*) 'land min/max=', rmin, rmax
read(lugi) label
read(lugi) fhour, idate
write(lugo1) label
write(lugo1) fhour,idate
write(6,*) 'bguess date ',fhour,' yr=',idate(4),' mo=',
1 idate(2),' day=',idate(3),' hr=',idate(1)
year = idate(4)
c Y2K, y2k problem
if (year.lt.200) year = year + 1900
month = idate(2)
day = idate(3)
hour = idate(1)
call pdays(year,month,day,factor)
c PRATE = kg/m/m/s
c x 86400 to get to mm/day
c x factor which is a correction for leap year
c x (nhour/24) to get correction for a n-hour cycle
c net result .. correction in mm
write(*,*) ' time factor= ',factor, ' should be approx 1'
factor = 86400.0 * factor * (nhour/24.0)
write(*,*) ' scaling factor= ',factor
read(lugi) sfctmp
read(lugi) soilm
read(lugi) snow
call maxmin(sfctmp,ijdim,rmin,rmax)
write(*,*) 'sfctmp min/max=', rmin, rmax
call maxmin(snow,ijdim,rmin,rmax)
write(*,*) 'snow min/max=', rmin, rmax
call maxmin(soilm(1,1),ijdim,rmin,rmax)
write(*,*) 'soilm(:,1) min/max=', rmin, rmax
call maxmin(soilm(1,2),ijdim,rmin,rmax)
write(*,*) 'soilm(:,2) min/max=', rmin, rmax
c now to do the correction
c
c v1.4
c
c case1:
c obspre .ge. mdlpre-runoff .and. runoff .gt. 0
c => assume extra precip would have gone into runoff
c => do nothing
c obspre .lt. mdlpre-runoff .and. runoff .gt. 0
c => apply precip correction
c obspre .ne. mdlpre .and. runoff .eq. 0
c => assume no runoff either way
c => apply precip correction
c
c surface temp < 0C .. do nothing .. frozen ground
c
c note:
c runoff == 0 => v1.0 scheme
c temperature correction
c correction is applied to top layer first
c
c write diagnostic file
c soilm(ijdim,lsoil)
c
cif (lsoil.eq.2) then
c write(*,*) '>>writing soilm'
c call setctr(7, 1, 80)
c call ezgbw(soilm(1,1),lonf,latg,.false.,.false.,
c 2 144,.false.,rlev1,year,month,day,hour,lugo2,
c 3 0,0,'hour','uninit anl')
c call ezgbw(soilm(1,2),lonf,latg,.false.,.false.,
c 2 144,.false.,rlev2,year,month,day,hour,lugo2,
c 3 0,0,'hour','uninit anl')
c write(*,*) '<<writing soilm'
cendif
do i = 1, ijdim
c check for water point
if (land(i).eq.0) goto 200
c check for undefined values of precip
if (abs((obspre(i)-undef)/undef).lt.1e-5) goto 200
c should be defined but check anyways
if (abs((mdlpre(i)-undef)/undef).lt.1e-5) goto 200
c check for surface temperature
if (sfctmp(i).lt.273.16) goto 200
c figure out the precip surplus
if (runoff(i).gt.0.0) then
offset = (obspre(i) - (mdlpre(i)-runoff(i)))*factor
if (offset.gt.0.0) offset = 0.0
else
offset = (obspre(i) - mdlpre(i))*factor
endif
if (offset.eq.0.0) goto 200
c
c do k = lsoil, 1, -1 => bottom to top
c do k = 1, lsoil => top to bottom
do k = 1, lsoil
smoist = soilm(i,k)*soilcap(k) + offset
if (smoist.lt.0.0) then
bias(i,k) = offset - smoist
offset = smoist
soilm(i,k) = 0.0
else if (smoist.gt.soilcap(k)) then
bias(i,k) = offset - (smoist - soilcap(k))
offset = smoist - soilcap(k)
soilm(i,k) = 1.0
else
bias(i,k) = offset
soilm(i,k) = smoist / soilcap(k)
goto 200
endif
enddo
200 continue
enddo
c
c write diagnostic file
c soilm(ijdim,lsoil)
c bias(ijdim,lsoil)
c
if (lsoil.eq.2) then
c write(*,*) '>>new diag'
c call ezgbw(soilm(1,1),lonf,latg,.false.,.false.,
c 2 144,.false.,rlev1,year,month,day,hour,lugo2,
c 3 0,0,'hour','init anl')
c call ezgbw(soilm(1,2),lonf,latg,.false.,.false.,
c 2 144,.false.,rlev2,year,month,day,hour,lugo2,
c 3 0,0,'hour','init anl')
c
call ezgbw1(bias(1,1),lonf,latg,234,.false.,rlev1,
2 year,month,day,hour,lugo2)
call ezgbw1(bias(1,2),lonf,latg,234,.false.,rlev2,
2 year,month,day,hour,lugo2)
call maxmin(bias(1,1),ijdim,rmin,rmax)
write(*,*) 'bias(:,1) min/max=', rmin, rmax
call maxmin(bias(1,2),ijdim,rmin,rmax)
write(*,*) 'bias(:,2) min/max=', rmin, rmax
call maxmin(soilm(1,1),ijdim,rmin,rmax)
write(*,*) 'new soilm(:,1) min/max=', rmin, rmax
call maxmin(soilm(1,2),ijdim,rmin,rmax)
write(*,*) 'new soilm(:,2) min/max=', rmin, rmax
endif
write(lugo1) sfctmp
write(lugo1) soilm
write(lugo1) snow
c copy rest of the file
read(lugi) tempv
write(lugo1) tempv
1000 read(lugi,end=1200) temp
write(lugo1) temp
goto 1000
1200 continue
stop
end
subroutine maxmin(a,idim,rmin,rmax)
parameter (undef=9.999e20) ! wgrib undefined number
parameter (undeflow=9.998e20) ! wgrib undefined number
real a(idim), rmax, rmin
rmax=-1e20
rmin= 1e20
do i = 1, idim
if (a(i).lt.undeflow) then
if (rmin.gt.a(i)) rmin = a(i)
if (rmax.lt.a(i)) rmax = a(i)
endif
enddo
return
end
subroutine pdays(year,month,day,factor)
c
c returns the scaling factor for 6 vs 5 days
c factor = (no. of days in preceeding pentad) /
c no of days in current pentad
c
integer year,month,day
real factor
factor = 1.0
if (month.eq.1 .or. month.ge.4) return
if (mod(year,4).ne.0) return
if (mod(year,100).eq.0 .and. mod(year,400).ne.0) return
c must be leap year
if (month.eq.2) then
if (day.ge.25) then
factor = 5.0 / 6.0
return
endif
return
endif
if (day.eq.1) then
factor = 5.0 / 6.0
return
endif
if (day.le.6) then
factor = 6.0 / 5.0
return
endif
return
end
c
c added variance as a type
c
c=================================================================
c
c change center/model id (in common block)
c
subroutine setctr(ictr, isctr, imodel)
c
integer ictr, isctr, imodel
integer center, subcen, model
common /gribc/center, subcen, model
center = ictr
subcen = isctr
model = imodel
return
end
c=================================================================
c
block data gribitc
integer center, subcen, model
common /gribc/center, subcen, model
c nmc
data center/7/
c reanalysis
data subcen/3/
c who knows
data model/195/
end
c
c=================================================================
subroutine ezgbw1(data,nx,ny,iparm,islola,level,
2 iyear,imonth,iday,ihour,iunit)
c
c
c ez grib write 1 W. Ebisuzaki
c
c sets many defaults - so you don't have to
c lola - grid, gaussian (future)
c writes only 1 field that is for an instantaneous time
c
c
c real data(nx,ny) - data to written as grib record
c integer iparm - grib parameter number ex. 33 -> U wind
c logical islola - .true. -> lola .false. -> gaussian grid
c real level - 0.0 -> 1.0 (sigma), 1.+ ->2000 (prs)
c -1 -> surface
c
c
real data(nx,ny), level
integer iparm,iyear,imonth,iday,ihour,iunit
logical islola, ldummy(1)
call ezgbw(data,nx,ny,ldummy,.false.,iparm,islola,level,
2 iyear,imonth,iday,ihour,iunit,0,0,'hour','inst')
return
end
c=================================================================
subroutine ezgbw2(data,nx,ny,iparm,islola,level,
2 iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
c
c ez grib write 2 W. Ebisuzaki
c
c sets many defaults - so you don't have to
c lola - grid, gaussian (future)
c writes only 1 field that is for a time average/accumulation/inst
c
c
real data(nx,ny), level
integer iparm,iyear,imonth,iday,ihour,iunit
logical islola, ldummy(1)
integer itime1,itime2
character*(*) timeun,type
call ezgbw(data,nx,ny,ldummy,.false., iparm,islola,level,
2 iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
return
end
c=================================================================
subroutine ezgbw(data,nx,ny,bit,hasbit,iparm,islola,level,
2 iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
c
c
c ez grib write W. Ebisuzaki
c
c sets many defaults - so you don't have to
c lola - grid, gaussian (future)
c
c real data(nx,ny) - data to written as grib record
c logical bit(*), hasbit - bitmap of defined values
c integer iparm - grib parameter number ex. 33 -> U wind
c logical islola - .true. -> lola .false. -> gaussian grid
c real level - 0.0 -> 1.0 (sigma), 1.+ ->2000 (prs)
c -1 -> surface
c integer itime1, itime2 - time average from time1 to time2
c character*(*) timeun - 'hour', 'day', 'month'
c character*(*) type - 'ave', 'acc', 'inst', 'climo', 'var'
c - 'init anl', 'uninit anl'
c note: type='climo'
c itime1 = length of mean
c itime2 = number of years of data
c
parameter (ilpds=28)
parameter (mxbit=16)
parameter (iptv=132)
parameter (mxnxny=384*190)
parameter (ngrib = 100 + ilpds + mxnxny*(mxbit+1)/8+10)
real data(nx,ny), grib(ngrib), level
integer iyear, imonth, iday, ihour
character*(*) timeun, type
logical hasbit, bit(*)
integer levelt, level1, level2, itime1, itime2
integer iftu, time1, time2, timerg, nave, nmiss
integer dscale(255)
integer gridty, lengrb
logical init, islola
integer center, subcen, model
common /gribc/center, subcen, model
data init/.false./
save init, dscale
if (nx*ny.gt.mxnxny) then
write(*,*) '**************'
write(*,*) '*** ezgbwr ***'
write(*,*) 'array too small: edit, recompile'
write(*,*) '*** Fatal Error ***'
call exit(99)
stop
endif
if (.not.init) then
call idsdef(iptv,dscale)
init = .true.
endif
nave = 0
nmiss = 0
c gridtyp 0 -> lola 4 -> gaussian grid
gridty = 0
if (.not.islola) gridty = 4
c set colatitude
if (.not.islola) then
if (ny.eq.94) then
c T62 Gaussian lat
colat1 = (90.0-88.542)*3.141592654/180.0
elseif (ny.eq.190) then
c T126 Gaussian lat
colat1 = (90.0-89.227)*3.141592654/180.0
else
write(*,*) '*** Missing Gaussian colat in gribit ***'
write(*,*) '*** please add ***'
write(*,*) '*** 0 used ***'
colat1 = 0.0
endif
endif
c set level
level1 = 0
if (level.gt.256*256) then
levelt = level / (256*256)
level1 = mod(int(level)/256, 256)
level2 = mod(int(level), 256)
elseif (level.ge.2000.0) then
c meters above ground - 2000
level2 = level - 2000.0
levelt = 105
elseif (level.gt.1.0) then
c must be pressure level
level2 = level
levelt = 100
elseif (level.ge.0.0 .and. level.le.1.0) then
c must be sigma level
level2 = 10000.0 * level
levelt = 107
elseif (level.eq.-1.0) then
c must be the surface
level2 = 0
levelt = 1
elseif (level.gt.-256.0.and.level.lt.-1.0) then
rlevel = -level
levelt = int(rlevel)
if (levelt.le.100.or.levelt.eq.102.or.levelt.eq.103 .or.
1 levelt.eq.105.or.levelt.eq.107.or.levelt.eq.109 .or.
1 levelt.eq.111.or.levelt.eq.113.or.levelt.eq.125 .or.
1 levelt.eq.160.or.levelt.eq.200.or.levelt.eq.201) then
level1 = 0
level2 = int(1000*(rlevel - int(rlevel)) + 0.5)
else
level2 = int(1000*(rlevel - int(rlevel)) + 0.5)
level1 = level2 / 256
level2 = mod(level2,256)
endif
else if (level.eq.-257.0) then
c sigma 0..1
levelt = 108
level2 = 100
else
c ????
write(*,*) 'ezgbw2: unknown level?', level
call exit(99)
stop
endif
c for an average/accum.
if (timeun.eq.'hour'.or.timeun.eq.'HOUR') then
iftu = 1
elseif (timeun.eq.'day'.or.timeun.eq.'DAY') then
iftu = 2
elseif (timeun.eq.'month'.or.timeun.eq.'MONTH') then
iftu = 3
elseif (timeun.eq.'year'.or.timeun.eq.'YEAR') then
iftu = 4
write(*,*) 'ezwgb: check year parameter!!'
stop
else
write(*,*) 'ezwgb bad timeun:',timeun
call exit(99)
endif
write(*,*) 'type = ',type
if (type.eq.'init anl' .or. type.eq.'INIT ANL') then
timerg = 1
time1 = 0
time2 = 0
else if (type.eq.'uninit anl' .or. type.eq.'UNINIT ANL') then
timerg = 0
time1 = 0
time2 = 0
else if (type.eq.'inst'.or.type.eq.'INST'.or.itime1.eq.
1 itime2) then
timerg = 10
time1 = 0
time2 = min0(itime1, itime2)
else if (type.eq.'acc'.or.type.eq.'ACC') then
timerg = 4
time1 = min0(itime1, itime2)
time2 = max0(itime1, itime2)
else if (type.eq.'ave'.or.type.eq.'AVE') then
timerg = 3
time1 = min0(itime1, itime2)
time2 = max0(itime1, itime2)
else if (type.eq.'climo'.or.type.eq.'CLIMO') then
c time1 = length of average
c time2 = number of years of data
timerg = 51
time1 = 0
time2 = itime1
nave = itime2
else if (type.eq.'var'.or.type.eq.'VAR') then
timerg = 118
time1 = min0(itime1, itime2)
time2 = max0(itime1, itime2)
else
write(*,*) 'ezwgb2 type:',type
call exit(99)
endif
ids = dscale(iparm)
if (type.eq.'var'.or.type.eq.'VAR') ids = 2*ids
c bitmap?
ibms = 0
if (hasbit) ibms = 1
call gribit(data,bit,gridty,nx,ny,mxbit,colat1,ilpds,
1 iptv,center,subcen,model,ibms,iparm,
2 levelt, level1, level2,
3 iyear, imonth, iday, ihour,
4 iftu, time1, time2, timerg,
5 nave, nmiss, ids, grib, lengrb, ierr)
if (ierr.eq.0) call binwrt(iunit, grib, lengrb)
return
end
c
c=====================================================================
c
subroutine ezgbwk(data,lbm,kpds,kgds,iunit)
c
c write grib record, kpds, kgds must be defined by caller
c
parameter (mxbit=16)
parameter (ilpds = 28)
parameter (mxnxny=384*190)
parameter (ngrib = 100 + ilpds + mxnxny*(mxbit+1)/8+10)
real data(*)
logical lbm(*)
integer kpds(*), kgds(*)
integer center, subcen, model
real grib(ngrib)
c data(nx,ny) -- data to write
c lbm(nx,ny) -- optional bitmap
c lola/gaussian
idrt = kgds(1)
c dim of data(nx,ny)
nx = kgds(2)
ny = kgds(3)
c maximum number of bits
c colatitude
c colat1 = (90e3 - (kpds(6)))*acos(-1.0)/(180e3)
colat1 = (90e3 - (kgds(4)))*acos(-1.0)/(180e3)
c length of pds (usually 28)
c ilpds
c parameter table (usually 1)
iptv = kpds(19)
c forecast center
center = kpds(1)
c subcenter
subcen = kpds(23)
c generating model
model = kpds(2)
c has bitmap?
ibms = mod(kpds(4)/64,2)
write(*,*) 'ezgbwk1 ibms=',ibms
c input parameter
ipu = kpds(5)
c level indicator
itl=kpds(6)
c old code
c if(itl.LE.100.OR.itl.EQ.102.OR.itl.EQ.103.OR.
c & itl.EQ.105.OR.itl.EQ.107.OR.itl.EQ.111.OR.
c & itl.EQ.160) then
c mod 2/28/95
if (itl.le.100 .or. itl.eq.102 .or. itl.eq.103 .or.
1 itl.eq.105 .or. itl.eq.107 .or. itl.eq.109 .or.
1 itl.eq.111 .or. itl.eq.113 .or. itl.eq.125 .or.
1 itl.eq.160 .or. itl.eq.200 .or. itl.eq.201) then
il1=0
il2=kpds(7)
else
il1=mod(kpds(7)/256,256)
il2=mod(kpds(7),256)
endif
c date
c iyr = kpds(8)
iyr = kpds(8) + (kpds(21)-1)*100
imo = kpds(9)
idy = kpds(10)
ihr = kpds(11)
c forecast unit, and times, and units
iftu = kpds(13)
ip1 = kpds(14)
ip2 = kpds(15)
itr = kpds(16)
c number in/missing from average
ina = kpds(17)
inm = kpds(20)
c integer decimal scaling
ids = kpds(22)
c write(*,*) 'idrt=',idrt,' nx/ny=',nx,ny,' mxbit=',mxbit
c write(*,*) 'colat1=',colat1,' ilpds=',ilpds,' iptv=',iptv
c write(*,*) 'icen=',icen,' igen=',igen,' ibms=',ibms,' param=',ipu
c write(*,*) 'times:',itl,il1,il2
c write(*,*) 'date:',iyr,imo,idy,ihr
c write(*,*) 'forcast times:',iftu,ip1,ip2,itr
c write(*,*) 'no. for ave=',ina,' missing=',inm,' dec scale=',ids
call gribit(data,lbm,idrt,nx,ny,mxbit,colat1,ilpds,iptv,
a center,subcen,model,
1 ibms,ipu,itl,il1,il2,iyr,imo,idy,ihr,iftu,ip1,ip2,itr,
2 ina,inm,ids,grib,lengrb,ierr)
if (ierr.eq.0) call binwrt(iunit, grib, lengrb)
return
end
c=====================================================================
subroutine binwrt(iunit,string,n)
c
c binary write
c
character*1 string(n)
write(iunit) string
return
end
c=====================================================================
CFPP$ NOCONCUR R
SUBROUTINE GRIBIT(F,LBM,IDRT,IM,JM,MXBIT,COLAT1,
& ILPDS,IPTV,ICEN,ISUBCN,IGEN,IBMS,IPU,ITL,
& IL1,IL2,
& IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
& GRIB,LGRIB,IERR)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: GRIBIT CREATE GRIB MESSAGE
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31
C
C ABSTRACT: CREATE A GRIB MESSAGE FROM A FULL FIELD.
C AT PRESENT, ONLY GLOBAL LATLON GRIDS AND GAUSSIAN GRIDS ARE ALLOWED.
C
C PROGRAM HISTORY LOG:
C 92-10-31 IREDELL
C
C USAGE: CALL GRIBIT(F,LBM,IDRT,IM,JM,MXBIT,COLAT1,
C & ILPDS,IPTV,ICEN,IGEN,IBMS,IPU,ITL,IL1,IL2,
C & IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
C & GRIB,LGRIB,IERR)
C INPUT ARGUMENT LIST:
C F - REAL (IM*JM) FIELD DATA TO PACK INTO GRIB MESSAGE
C LBM - LOGICAL (IM*JM) BITMAP TO USE IF IBMS=1
C IDRT - INTEGER DATA REPRESENTATION TYPE
C (0 FOR LATLON OR 4 FOR GAUSSIAN)
C IM - INTEGER LONGITUDINAL DIMENSION
C JM - INTEGER LATITUDINAL DIMENSION
C MXBIT - INTEGER MAXIMUM NUMBER OF BITS TO USE (0 FOR NO LIMIT)
C COLAT1 - REAL FIRST COLATITUDE OF GAUSSIAN GRID IF IDRT=4
C ILPDS - INTEGER LENGTH OF THE PDS (USUALLY 28)
C IPTV - INTEGER PARAMETER TABLE VERSION (USUALLY 1)
C ICEN - INTEGER FORECAST CENTER (USUALLY 7)
C ISUBCN - INTEGER SUBCENTER
C IGEN - INTEGER MODEL GENERATING CODE
C IBMS - INTEGER BITMAP FLAG (0 FOR NO BITMAP)
C IPU - INTEGER PARAMETER AND UNIT INDICATOR
C ITL - INTEGER TYPE OF LEVEL INDICATOR
C IL1 - INTEGER FIRST LEVEL VALUE (0 FOR SINGLE LEVEL)
C IL2 - INTEGER SECOND LEVEL VALUE
C & IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
C & GRIB,LGRIB,IERR)
C IYR - INTEGER YEAR (YYYY)
C IMO - INTEGER MONTH
C IDY - INTEGER DAY
C IHR - INTEGER HOUR
C IFTU - INTEGER FORECAST TIME UNIT (1 FOR HOUR)
C IP1 - INTEGER FIRST TIME PERIOD
C IP2 - INTEGER SECOND TIME PERIOD (0 FOR SINGLE PERIOD)
C ITR - INTEGER TIME RANGE INDICATOR (10 FOR SINGLE PERIOD)
C INA - INTEGER NUMBER INCLUDED IN AVERAGE
C INM - INTEGER NUMBER MISSING FROM AVERAGE
C IDS - INTEGER DECIMAL SCALING
C
C OUTPUT ARGUMENT LIST:
C GRIB - CHARACTER (LGRIB) GRIB MESSAGE
C LGRIB - INTEGER LENGTH OF GRIB MESSAGE
C (NO MORE THAN 100+ILPDS+IM*JM*(MXBIT+1)/8)
C IERR - INTEGER ERROR CODE (0 FOR SUCCESS)
C
C SUBPROGRAMS CALLED:
C GTBITS - COMPUTE NUMBER OF BITS AND ROUND DATA APPROPRIATELY
C W3FI72 - ENGRIB DATA INTO A GRIB1 MESSAGE
C
C ATTRIBUTES:
C LANGUAGE: CRAY FORTRAN
C
C$$$
REAL F(IM*JM)
LOGICAL LBM(IM*JM)
CHARACTER GRIB(*)
INTEGER IPDS(25),IGDS(18),IBDS(9)
C
PARAMETER(JLPDS=28)
C INTEGER IBM(IM*JM*IBMS+1-IBMS)
C REAL FR(IM*JM)
INTEGER IBM(192 * 94)
REAL FR(192 * 94)
CHARACTER PDS(JLPDS)
C
c write(*,*) '>>gribit'
NF=IM*JM
LONI=NINT(360.E3/IM)
IF(IDRT.EQ.0) THEN
LAT1=NINT(90.E3)
LATI=NINT(180.E3/(JM-1))
IF(IM.EQ.144.AND.JM.EQ.73) THEN
IGRID=2
ELSEIF(IM.EQ.360.AND.JM.EQ.181) THEN
IGRID=3
ELSE
IGRID=255
ENDIF
ELSEIF(IDRT.EQ.4) THEN
LAT1=NINT(90.E3-180.E3/ACOS(-1.)*COLAT1)
LATI=JM/2
IF(IM.EQ.192.AND.JM.EQ.94) THEN
IGRID=98
ELSE
IGRID=255
ENDIF
ELSE
IERR=40
RETURN
ENDIF
IPDS(01)=ILPDS ! LENGTH OF PDS
IPDS(02)=IPTV ! PARAMETER TABLE VERSION ID
IPDS(03)=ICEN ! CENTER ID
IPDS(04)=IGEN ! GENERATING MODEL ID
IPDS(05)=IGRID ! GRID ID
IPDS(06)=1 ! GDS FLAG
IPDS(07)=IBMS ! BMS FLAG
IPDS(08)=IPU ! PARAMETER UNIT ID
IPDS(09)=ITL ! TYPE OF LEVEL ID
IPDS(10)=IL1 ! LEVEL 1 OR 0
IPDS(11)=IL2 ! LEVEL 2
c IPDS(12)=IYR ! YEAR
IPDS(12)=mod((IYR-1),100) + 1
IPDS(13)=IMO ! MONTH
IPDS(14)=IDY ! DAY
IPDS(15)=IHR ! HOUR
IPDS(16)=0 ! MINUTE
IPDS(17)=IFTU ! FORECAST TIME UNIT ID
IPDS(18)=IP1 ! TIME PERIOD 1
IPDS(19)=IP2 ! TIME PERIOD 2 OR 0
IPDS(20)=ITR ! TIME RANGE INDICATOR
IPDS(21)=INA ! NUMBER IN AVERAGE
IPDS(22)=INM ! NUMBER MISSING
c IPDS(23)=20 ! CENTURY
IPDS(23)=((IYR-IPDS(12))/100) + 1
IPDS(24)=ISUBCN ! SUBCENTER
IPDS(25)=IDS ! DECIMAL SCALING
IGDS(01)=0 ! NUMBER OF VERTICAL COORDS
IGDS(02)=255 ! VERTICAL COORD FLAG
IGDS(03)=IDRT ! DATA REPRESENTATION TYPE
IGDS(04)=IM ! EAST-WEST POINTS
IGDS(05)=JM ! NORTH-SOUTH POINTS
IGDS(06)=LAT1 ! LATITUDE OF ORIGIN
IGDS(07)=0 ! LONGITUDE OF ORIGIN
IGDS(08)=128 ! RESOLUTION FLAG
IGDS(09)=-LAT1 ! LATITUDE OF END
IGDS(10)=-LONI ! LONGITUDE OF END
IGDS(11)=LATI ! LAT INCREMENT OR GAUSSIAN LATS
IGDS(12)=LONI ! LONGITUDE INCREMENT
IGDS(13)=0 ! SCANNING MODE FLAGS
DO I=14,18
IGDS(I)=0 ! NOT USED
ENDDO
DO I=1,9
IBDS(I)=0 ! BDS FLAGS
ENDDO
NBM=NF
c write(*,*) 'gribit:ibms=',ibms
IF(IBMS.NE.0) THEN
NBM=0
DO I=1,NF
IF(LBM(I)) THEN
IBM(I)=1
NBM=NBM+1
ELSE
IBM(I)=0
ENDIF
ENDDO
IF(NBM.EQ.NF) IPDS(7)=0
ENDIF
c write(*,*) 'gribit2 nbm=',nbm
IF(NBM.EQ.0) THEN
DO I=1,NF
FR(I)=0.
ENDDO
NBIT=0
ELSE
fmin = f(1)
fmax = fmin
do ii = 1, nf
fmin = amin1(fmin,f(ii))
fmax = amax1(fmax,f(ii))
enddo
c write(*,*) 'fmin/fmax=',fmin,fmax
CALL GTBITS(IPDS(7),IDS,NF,IBM,F,FR,FMIN,FMAX,NBIT)
c write(*,*) 'fmin/max=',fmin,fmax,' nbit=',nbit
IF(MXBIT.GT.0) NBIT=MIN(NBIT,MXBIT)
ENDIF
c write(*,*)'gribit4'
CALL W3FI72(0,FR,0,NBIT,0,IPDS,PDS,
& 1,255,IGDS,0,0,IBM,NF,IBDS,
& NFO,GRIB,LGRIB,IERR)
c write(*,*)'gribit5'
RETURN
END
CFPP$ NOCONCUR R
SUBROUTINE GTBITS(IBM,IDS,LEN,MG,G,GROUND,GMIN,GMAX,NBIT)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: GTBITS COMPUTE NUMBER OF BITS AND ROUND FIELD.
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31
C
C ABSTRACT: THE NUMBER OF BITS REQUIRED TO PACK A GIVEN FIELD
C AT A PARTICULAR DECIMAL SCALING IS COMPUTED USING THE FIELD RANGE.
C THE FIELD IS ROUNDED OFF TO THE DECIMAL SCALING FOR PACKING.
C THE MINIMUM AND MAXIMUM ROUNDED FIELD VALUES ARE ALSO RETURNED.
C GRIB BITMAP MASKING FOR VALID DATA IS OPTIONALLY USED.
C
C PROGRAM HISTORY LOG:
C 92-10-31 IREDELL
C
C USAGE: CALL GTBITS(IBM,IDS,LEN,MG,G,GMIN,GMAX,NBIT)
C INPUT ARGUMENT LIST:
C IBM - INTEGER BITMAP FLAG (=0 FOR NO BITMAP)
C IDS - INTEGER DECIMAL SCALING
C (E.G. IDS=3 TO ROUND FIELD TO NEAREST MILLI-VALUE)
C LEN - INTEGER LENGTH OF THE FIELD AND BITMAP
C MG - INTEGER (LEN) BITMAP IF IBM=1 (0 TO SKIP, 1 TO KEEP)
C G - REAL (LEN) FIELD
C
C OUTPUT ARGUMENT LIST:
C GROUND - REAL (LEN) FIELD ROUNDED TO DECIMAL SCALING
C GMIN - REAL MINIMUM VALID ROUNDED FIELD VALUE
C GMAX - REAL MAXIMUM VALID ROUNDED FIELD VALUE
C NBIT - INTEGER NUMBER OF BITS TO PACK
C
C SUBPROGRAMS CALLED:
C ISRCHNE - FIND FIRST VALUE IN AN ARRAY NOT EQUAL TO TARGET VALUE
C
C ATTRIBUTES:
C LANGUAGE: CRAY FORTRAN
C
C$$$
DIMENSION MG(LEN),G(LEN),GROUND(LEN)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON
DS=10.**IDS
IF(IBM.EQ.0) THEN
GROUND(1)=NINT(G(1)*DS)/DS
GMAX=GROUND(1)
GMIN=GROUND(1)
DO I=2,LEN
GROUND(I)=NINT(G(I)*DS)/DS
GMAX=MAX(GMAX,GROUND(I))
GMIN=MIN(GMIN,GROUND(I))
ENDDO
ELSE
I1=ISRCHNE(LEN,MG,1,0)
IF(I1.GT.0.AND.I1.LE.LEN) THEN
GROUND(I1)=NINT(G(I1)*DS)/DS
GMAX=GROUND(I1)
GMIN=GROUND(I1)
DO I=I1+1,LEN
IF(MG(I).NE.0) THEN
GROUND(I)=NINT(G(I)*DS)/DS
GMAX=MAX(GMAX,GROUND(I))
GMIN=MIN(GMIN,GROUND(I))
ENDIF
ENDDO
ELSE
GMAX=0.
GMIN=0.
ENDIF
ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C COMPUTE NUMBER OF BITS
NBIT=LOG((GMAX-GMIN)*DS+0.9)/LOG(2.)+1.
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
RETURN
END
SUBROUTINE IDSDEF(IPTV,IDS)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: IDSDEF SETS DEFAULT DECIMAL SCALINGS
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31
C
C ABSTRACT: SETS DECIMAL SCALINGS DEFAULTS FOR VARIOUS PARAMETERS.
C A DECIMAL SCALING OF -3 MEANS DATA IS PACKED IN KILO-SI UNITS.
C
C PROGRAM HISTORY LOG:
C 92-10-31 IREDELL
C
C USAGE: CALL IDSDEF(IPTV,IDS)
C INPUT ARGUMENTS:
C IPTV PARAMTER TABLE VERSION (ONLY 1 OR 2 IS RECOGNIZED)
C OUTPUT ARGUMENTS:
C IDS INTEGER (255) DECIMAL SCALINGS
C (UNKNOWN DECIMAL SCALINGS WILL NOT BE SET)
C
C ATTRIBUTES:
C LANGUAGE: CRAY FORTRAN
C
C$$$
DIMENSION IDS(255)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IF(IPTV.EQ.1.OR.IPTV.EQ.2.or.IPTV.eq.132) THEN
IDS(001)=-1 ! PRESSURE (PA)
IDS(002)=-1 ! SEA-LEVEL PRESSURE (PA)
IDS(003)=5 ! PRESSURE TENDENCY (PA/S)
!
!
IDS(006)=-1 ! GEOPOTENTIAL (M2/S2)
IDS(007)=0 ! GEOPOTENTIAL HEIGHT (M)
IDS(008)=0 ! GEOMETRIC HEIGHT (M)
IDS(009)=0 ! STANDARD DEVIATION OF HEIGHT (M)
!
IDS(011)=1 ! TEMPERATURE (K)
IDS(012)=1 ! VIRTUAL TEMPERATURE (K)
IDS(013)=1 ! POTENTIAL TEMPERATURE (K)
IDS(014)=1 ! PSEUDO-ADIABATIC POTENTIAL TEMPERATURE (K)
IDS(015)=1 ! MAXIMUM TEMPERATURE (K)
IDS(016)=1 ! MINIMUM TEMPERATURE (K)
IDS(017)=1 ! DEWPOINT TEMPERATURE (K)
IDS(018)=1 ! DEWPOINT DEPRESSION (K)
IDS(019)=4 ! TEMPERATURE LAPSE RATE (K/M)
IDS(020)=0 ! VISIBILITY (M)
! RADAR SPECTRA 1 ()
! RADAR SPECTRA 2 ()
! RADAR SPECTRA 3 ()
!
IDS(025)=1 ! TEMPERATURE ANOMALY (K)
IDS(026)=-1 ! PRESSURE ANOMALY (PA)
IDS(027)=0 ! GEOPOTENTIAL HEIGHT ANOMALY (M)
! WAVE SPECTRA 1 ()
! WAVE SPECTRA 2 ()
! WAVE SPECTRA 3 ()
IDS(031)=0 ! WIND DIRECTION (DEGREES)
IDS(032)=1 ! WIND SPEED (M/S)
IDS(033)=1 ! ZONAL WIND (M/S)
IDS(034)=1 ! MERIDIONAL WIND (M/S)
IDS(035)=-4 ! STREAMFUNCTION (M2/S) WNE
IDS(036)=-4 ! VELOCITY POTENTIAL (M2/S) WNE
IDS(037)=-1 ! MONTGOMERY STREAM FUNCTION (M2/S2)
IDS(038)=8 ! SIGMA VERTICAL VELOCITY (1/S)
IDS(039)=3 ! PRESSURE VERTICAL VELOCITY (PA/S)
IDS(040)=4 ! GEOMETRIC VERTICAL VELOCITY (M/S)
IDS(041)=6 ! ABSOLUTE VORTICITY (1/S)
IDS(042)=6 ! ABSOLUTE DIVERGENCE (1/S)
IDS(043)=6 ! RELATIVE VORTICITY (1/S)
IDS(044)=6 ! RELATIVE DIVERGENCE (1/S)
IDS(045)=4 ! VERTICAL U SHEAR (1/S)
IDS(046)=4 ! VERTICAL V SHEAR (1/S)
IDS(047)=0 ! DIRECTION OF CURRENT (DEGREES)
! SPEED OF CURRENT (M/S)
! U OF CURRENT (M/S)
! V OF CURRENT (M/S)
IDS(051)=4 ! SPECIFIC HUMIDITY (KG/KG)
IDS(052)=0 ! RELATIVE HUMIDITY (PERCENT)
IDS(053)=4 ! HUMIDITY MIXING RATIO (KG/KG)
IDS(054)=1 ! PRECIPITABLE WATER (KG/M2)
IDS(055)=-1 ! VAPOR PRESSURE (PA)
IDS(056)=-1 ! SATURATION DEFICIT (PA)
IDS(057)=1 ! EVAPORATION (KG/M2)
IDS(058)=1 ! CLOUD ICE (KG/M2)
IDS(059)=6 ! PRECIPITATION RATE (KG/M2/S)
IDS(060)=0 ! THUNDERSTORM PROBABILITY (PERCENT)
IDS(061)=1 ! TOTAL PRECIPITATION (KG/M2)
IDS(062)=1 ! LARGE-SCALE PRECIPITATION (KG/M2)
IDS(063)=1 ! CONVECTIVE PRECIPITATION (KG/M2)
IDS(064)=6 ! WATER EQUIVALENT SNOWFALL RATE (KG/M2/S)
IDS(065)=0 ! WATER EQUIVALENT OF SNOW DEPTH (KG/M2)
IDS(066)=2 ! SNOW DEPTH (M)
! MIXED-LAYER DEPTH (M)
! TRANSIENT THERMOCLINE DEPTH (M)
! MAIN THERMOCLINE DEPTH (M)
! MAIN THERMOCLINE ANOMALY (M)
IDS(071)=0 ! TOTAL CLOUD COVER (PERCENT)
IDS(072)=0 ! CONVECTIVE CLOUD COVER (PERCENT)
IDS(073)=0 ! LOW CLOUD COVER (PERCENT)
IDS(074)=0 ! MIDDLE CLOUD COVER (PERCENT)
IDS(075)=0 ! HIGH CLOUD COVER (PERCENT)
IDS(076)=1 ! CLOUD WATER (KG/M2)
!
IDS(078)=1 ! CONVECTIVE SNOW (KG/M2)
IDS(079)=1 ! LARGE SCALE SNOW (KG/M2)
IDS(080)=1 ! WATER TEMPERATURE (K)
IDS(081)=0 ! SEA-LAND MASK ()
! DEVIATION OF SEA LEVEL FROM MEAN (M)
IDS(083)=5 ! ROUGHNESS (M)
IDS(084)=0 ! ALBEDO (PERCENT)
IDS(085)=1 ! SOIL TEMPERATURE (K)
IDS(086)=0 ! SOIL WETNESS (KG/M2)
IDS(087)=0 ! VEGETATION (PERCENT)
! SALINITY (KG/KG)
IDS(089)=4 ! DENSITY (KG/M3)
IDS(090)=1 ! RUNOFF (KG/M2)
IDS(091)=0 ! ICE CONCENTRATION ()
! ICE THICKNESS (M)
IDS(093)=0 ! DIRECTION OF ICE DRIFT (DEGREES)
! SPEED OF ICE DRIFT (M/S)
! U OF ICE DRIFT (M/S)
! V OF ICE DRIFT (M/S)
! ICE GROWTH (M)
! ICE DIVERGENCE (1/S)
IDS(099)=1 ! SNOW MELT (KG/M2)
! SIG HEIGHT OF WAVES AND SWELL (M)
IDS(101)=0 ! DIRECTION OF WIND WAVES (DEGREES)
! SIG HEIGHT OF WIND WAVES (M)
! MEAN PERIOD OF WIND WAVES (S)
IDS(104)=0 ! DIRECTION OF SWELL WAVES (DEGREES)
! SIG HEIGHT OF SWELL WAVES (M)
! MEAN PERIOD OF SWELL WAVES (S)
IDS(107)=0 ! PRIMARY WAVE DIRECTION (DEGREES)
! PRIMARY WAVE MEAN PERIOD (S)
IDS(109)=0 ! SECONDARY WAVE DIRECTION (DEGREES)
! SECONDARY WAVE MEAN PERIOD (S)
IDS(111)=0 ! NET SOLAR RADIATIVE FLUX AT SURFACE (W/M2)
IDS(112)=0 ! NET LONGWAVE RADIATIVE FLUX AT SURFACE (W/M2)
IDS(113)=0 ! NET SOLAR RADIATIVE FLUX AT TOP (W/M2)
IDS(114)=0 ! NET LONGWAVE RADIATIVE FLUX AT TOP (W/M2)
IDS(115)=0 ! NET LONGWAVE RADIATIVE FLUX (W/M2)
IDS(116)=0 ! NET SOLAR RADIATIVE FLUX (W/M2)
IDS(117)=0 ! TOTAL RADIATIVE FLUX (W/M2)
!
!
!
IDS(121)=0 ! LATENT HEAT FLUX (W/M2)
IDS(122)=0 ! SENSIBLE HEAT FLUX (W/M2)
IDS(123)=0 ! BOUNDARY LAYER DISSIPATION (W/M2)
IDS(124)=3 ! U WIND STRESS (N/M2)
IDS(125)=3 ! V WIND STRESS (N/M2)
! WIND MIXING ENERGY (J)
! IMAGE DATA ()
IDS(128)=-1 ! MEAN SEA-LEVEL PRESSURE (STDATM) (PA)
IDS(129)=-1 ! MEAN SEA-LEVEL PRESSURE (MAPS) (PA)
IDS(130)=-1 ! MEAN SEA-LEVEL PRESSURE (ETA) (PA)
IDS(131)=1 ! SURFACE LIFTED INDEX (K)
IDS(132)=1 ! BEST LIFTED INDEX (K)
IDS(133)=1 ! K INDEX (K)
IDS(134)=1 ! SWEAT INDEX (K)
IDS(135)=10 ! HORIZONTAL MOISTURE DIVERGENCE (KG/KG/S)
IDS(136)=4 ! SPEED SHEAR (1/S)
IDS(137)=5 ! 3-HR PRESSURE TENDENCY (PA/S)
IDS(138)=6 ! BRUNT-VAISALA FREQUENCY SQUARED (1/S2)
IDS(139)=11 ! POTENTIAL VORTICITY (MASS-WEIGHTED) (1/S/M)
IDS(140)=0 ! RAIN MASK ()
IDS(141)=0 ! FREEZING RAIN MASK ()
IDS(142)=0 ! ICE PELLETS MASK ()
IDS(143)=0 ! SNOW MASK ()
IDS(144)=4 ! SOILW
!
!
!
!
!
!
! COVARIANCE BETWEEN V AND U (M2/S2)
! COVARIANCE BETWEEN U AND T (K*M/S)
! COVARIANCE BETWEEN V AND T (K*M/S)
!
!
IDS(155)=0 ! GROUND HEAT FLUX (W/M2)
IDS(156)=0 ! CONVECTIVE INHIBITION (W/M2)
! CONVECTIVE APE (J/KG)
! TURBULENT KE (J/KG)
! CONDENSATION PRESSURE OF LIFTED PARCEL (PA)
IDS(160)=0 ! CLEAR SKY UPWARD SOLAR FLUX (W/M2)
IDS(161)=0 ! CLEAR SKY DOWNWARD SOLAR FLUX (W/M2)
IDS(162)=0 ! CLEAR SKY UPWARD LONGWAVE FLUX (W/M2)
IDS(163)=0 ! CLEAR SKY DOWNWARD LONGWAVE FLUX (W/M2)
IDS(164)=0 ! CLOUD FORCING NET SOLAR FLUX (W/M2)
IDS(165)=0 ! CLOUD FORCING NET LONGWAVE FLUX (W/M2)
IDS(166)=0 ! VISIBLE BEAM DOWNWARD SOLAR FLUX (W/M2)
IDS(167)=0 ! VISIBLE DIFFUSE DOWNWARD SOLAR FLUX (W/M2)
IDS(168)=0 ! NEAR IR BEAM DOWNWARD SOLAR FLUX (W/M2)
IDS(169)=0 ! NEAR IR DIFFUSE DOWNWARD SOLAR FLUX (W/M2)
!
IDS(170)=3 ! old (but current as of 1/94) u-stress (n/m**2)
IDS(171)=3 ! old (but current as of 1/94) v-stress (n/m**2)
!
!
IDS(172)=3 ! MOMENTUM FLUX (N/M2)
IDS(173)=0 ! MASS POINT MODEL SURFACE ()
IDS(174)=0 ! VELOCITY POINT MODEL SURFACE ()
IDS(175)=0 ! SIGMA LAYER NUMBER ()
IDS(176)=2 ! LATITUDE (DEGREES)
IDS(177)=2 ! EAST LONGITUDE (DEGREES)
IDS(178)=2 ! UMAS wne
IDS(179)=2 ! VMAS wne
!
!
!
IDS(181)=9 ! X-GRADIENT LOG PRESSURE (1/M)
IDS(182)=9 ! Y-GRADIENT LOG PRESSURE (1/M)
IDS(183)=5 ! X-GRADIENT HEIGHT (M/M)
IDS(184)=5 ! Y-GRADIENT HEIGHT (M/M)
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
IDS(201)=0 ! ICE-FREE WATER SURCACE (PERCENT)
!
!
IDS(204)=0 ! DOWNWARD SOLAR RADIATIVE FLUX (W/M2)
IDS(205)=0 ! DOWNWARD LONGWAVE RADIATIVE FLUX (W/M2)
!
IDS(207)=0 ! MOISTURE AVAILABILITY (PERCENT)
! EXCHANGE COEFFICIENT (KG/M2/S)
IDS(209)=0 ! NUMBER OF MIXED LAYER NEXT TO SFC ()
!
IDS(211)=0 ! UPWARD SOLAR RADIATIVE FLUX (W/M2)
IDS(212)=0 ! UPWARD LONGWAVE RADIATIVE FLUX (W/M2)
IDS(213)=0 ! NON-CONVECTIVE CLOUD COVER (PERCENT)
IDS(214)=6 ! CONVECTIVE PRECIPITATION RATE (KG/M2/S)
IDS(215)=7 ! TOTAL DIABATIC HEATING RATE (K/S)
IDS(216)=7 ! TOTAL RADIATIVE HEATING RATE (K/S)
IDS(217)=7 ! TOTAL DIABATIC NONRADIATIVE HEATING RATE (K/S)
IDS(218)=2 ! PRECIPITATION INDEX (FRACTION)
IDS(219)=1 ! STD DEV OF IR T OVER 1X1 DEG AREA (K)
IDS(220)=4 ! NATURAL LOG OF SURFACE PRESSURE OVER 1 KPA ()
!
IDS(222)=0 ! 5-WAVE GEOPOTENTIAL HEIGHT (M)
IDS(223)=1 ! PLANT CANOPY SURFACE WATER (KG/M2)
!
!
! BLACKADARS MIXING LENGTH (M)
! ASYMPTOTIC MIXING LENGTH (M)
IDS(228)=1 ! POTENTIAL EVAPORATION (KG/M2)
IDS(229)=0 ! SNOW PHASE-CHANGE HEAT FLUX (W/M2)
!
IDS(231)=3 ! CONVECTIVE CLOUD MASS FLUX (PA/S)
IDS(232)=0 ! DOWNWARD TOTAL RADIATION FLUX (W/M2)
IDS(233)=0 ! UPWARD TOTAL RADIATION FLUX (W/M2)
IDS(224)=1 ! BASEFLOW-GROUNDWATER RUNOFF (KG/M2)
IDS(225)=1 ! STORM SURFACE RUNOFF (KG/M2)
IDS(234)=1 ! BASEFLOW-GROUNDWATER RUNOFF (KG/M2)
!
!
IDS(238)=0 ! SNOW COVER (PERCENT)
IDS(239)=1 ! SNOW TEMPERATURE (K)
!
IDS(241)=7 ! LARGE SCALE CONDENSATION HEATING RATE (K/S)
IDS(242)=7 ! DEEP CONVECTIVE HEATING RATE (K/S)
IDS(243)=10 ! DEEP CONVECTIVE MOISTENING RATE (KG/KG/S)
IDS(244)=7 ! SHALLOW CONVECTIVE HEATING RATE (K/S)
IDS(245)=10 ! SHALLOW CONVECTIVE MOISTENING RATE (KG/KG/S)
IDS(246)=7 ! VERTICAL DIFFUSION HEATING RATE (KG/KG/S)
IDS(247)=7 ! VERTICAL DIFFUSION ZONAL ACCELERATION (M/S/S)
IDS(248)=7 ! VERTICAL DIFFUSION MERID ACCELERATION (M/S/S)
IDS(249)=10 ! VERTICAL DIFFUSION MOISTENING RATE (KG/KG/S)
IDS(250)=7 ! SOLAR RADIATIVE HEATING RATE (K/S)
IDS(251)=7 ! LONGWAVE RADIATIVE HEATING RATE (K/S)
! DRAG COEFFICIENT ()
! FRICTION VELOCITY (M/S)
! RICHARDSON NUMBER ()
!
ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
RETURN
END
FUNCTION ISRCHNE(N,IX,INCX,ITARGET)
INTEGER IX(*),ITARGET
J=1
ISRCHNE=0
IF(N.LE.0) RETURN
IF(INCX.LT.0) J=1-(N-1)*INCX
DO I=1,N
IF(IX(J).NE.ITARGET) THEN
ISRCHNE=I
RETURN
ENDIF
J=J+INCX
ENDDO
RETURN
END