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