subroutine inter(yrdoy,flat,flon,pteran,pcloud,icat,isno) c c -- author: c -- lucy liu / 26-feb-96 c -- date c -- 26-feb-96 c -- purpose: c -- given date and latitude, longitude, find terrain, c -- cloud pressures, snow/ice, and surface category c -- from appropriate data bases c save c logical frscal/.true./, leap integer lusur, lucat, bindoy(12), yr, doy, yrdoy real flat, flon, pteran, pcloud character afile*70 data bindoy /15,46,74,105,135,166,196,227,258,288,319,349/ c c --- leap is a statement function. leap (iyr) is true if iyr is a leap c --- year (can be divided by 4 and not 100 or by 400). c leap (iyr) = ((mod (iyr,4).eq.0) .and. (mod (iyr,100).ne.0)) .or. + (mod (iyr,400).eq.0) c if(frscal) then c afile='/usr/people/cgw/source/overpass/TERPRS.DAT' call get_lun(lusur) if (lusur .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' TERPRS.DAT' stop endif open(lusur,file=afile,form='unformatted',status='old',err=900) c afile='/usr/people/cgw/source/overpass/SURFCAT.DAT' call get_lun(lucat) if (lucat .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' SURFCAT.DAT' stop endif open(lucat,file=afile,form='formatted',status='old',err=900) c frscal=.false. c endif c call getsur(lusur,lucat, .false.) c write(6,*) 'enter yr doy and lat lon in degree:' c open (42,file='pres_in',status='unknown',form='unformatted') c open (43,file='pres_out',status='unknown',form='unformatted') c read (42) yrdoy, flat, flon c close (42) yr = int(yrdoy/1000) doy = mod(yrdoy,1000) c write (6,*) yrdoy,yr,doy if (leap (yr )) then do 50 mon = 3, 12 bindoy(mon) = bindoy(mon) + 1 50 continue endif do 200 mon = 1, 12 if (doy.le.bindoy(mon)) go to 250 200 continue mon = 1 250 month1 = mon month2 = mon - 1 if (month2.eq.0) month2 = 12 if (month1.ne.monpre) then c c --- read new snow-ice cover data when there is month change c print*,'calling getcov' call getcov (month1,month2, lfatal) monpre = month1 endif c lat=nint(flat*3.14159/180.*1.e4) lon=nint(flon*3.14159/180.*1.e4) idif = bindoy(month1) - doy if (month2.eq.12 .and. doy.gt.bindoy(12) .and. + doy.le.bindoy(12)+16) idif = 31 - (doy - bindoy(12)) call insurf (lat,lon, iter,icat) call incovr (idif,lat,lon, isno,icld) pteran=iter/1013. pcloud=icld/1013. c write (43) iter/1013., icat, isno, icld/1013. c close (43) write(6,*) 'inter',iter/1013., icat, isno, icld/1013. c c close(lusur) c close(lucat) return 900 write(6,*) 'error in opening ',afile end c subroutine getcov (month1,month2, lfatal) c ********************************************************************** c c module name: getcov c c version : 1.0 c c programmer : lucy liu hstx. 1996 c modified by: c c purpose : read snow-ice cover data for current & previous months. c read cloud pressure data for current & previous months. c c note : month n's snow-ice data cover/cloud pressure data from c middle of month n-1 to the middle of month n. c c method : c 1. open the data set/file. c 2. read two months data and store them in arrays map1 and map2, c cld1 and cld2 in common area covpre. c c calling routine: main c c calling sequence: getcov (month1,month2, lfatal) c c subroutines called: none c c common blocks used: covpre c c arguments and local variables: c c name type i/o descriptions c ------ ----- --- ----------------------------------------------- c month1 i*4 i current month index c month2 i*4 i previous month index c lfatal l*4 o flag for fatal error c lucov i*4 l unit number for reading snow-ice data file c tmp r*4 l buffer for each record read c (363) c mname c*3 l member name of a partitioned data set in ibmmvs, c (12) or extention of file name in other environments c dsn1 c*44 l data set name for current month c dsn2 c*44 l data set name for previous month c ********************************************************************** c --- common areas ----------------------------------------------------- c integer monpre real map1(363,180), map2(363,180), cld1(144,72),cld2(144,72) common /covpre/ monpre, map1, map2, cld1, cld2 c c --- arguments -------------------------------------------------------- c integer month1, month2 logical lfatal c c --- local variables -------------------------------------------------- c real*4 tmpsno(363), tmpcld(144) integer lucov character mname(12)*3, afile*45 c data mname /'JAN','FEB','MAR','APR','MAY','JUN', + 'JUL','AUG','SEP','OCT','NOV','DEC'/ c lfatal = .false. c afile = '/usr/people/cgw/source/overpass/SNOICE.'//mname(month1) c c --- read current month data c call get_lun(lucov) if (lucov .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' SNOICE.'//mname(month1) stop endif open (lucov,file=afile, status='old', form='unformatted', + action='read', err=910) do 250 j = 1, 180 read (lucov, end=920, err=930) tmpsno do 200 k = 1, 363 map1(k,j) = tmpsno(k) 200 continue 250 continue close(lucov) c c --- read previous month data c afile = '/usr/people/cgw/source/overpass/SNOICE.'//mname(month2) call get_lun(lucov) if (lucov .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' SNOICE.'//mname(month2) stop endif open (lucov, file=afile, status='old', form='unformatted', + action='read', err=910) do 350 j = 1, 180 read (lucov, end=920, err=930) tmpsno do 300 k = 1, 363 map2(k,j) = tmpsno(k) 300 continue 350 continue close(lucov) c afile = '/usr/people/cgw/source/overpass/CLDPRES.'//mname(month1) c c --- read current month data c call get_lun(lucld) if (lucld .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' CLDPRES.'//mname(month1) stop endif open (lucld,file=afile, status='old', form='unformatted', + action='read', err=910) do 450 j = 1, 72 read (lucld, end=920, err=930) tmpcld do 400 k = 1, 144 cld1(k,j) = tmpcld(k) 400 continue 450 continue close(lucld) c c --- read previous month data c afile = '/usr/people/cgw/source/overpass/CLDPRES.'//mname(month2) call get_lun(lucld) if (lucld .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' CLDPRES.'//mname(month2) stop endif open (lucld, file=afile, status='old', form='unformatted', + action='read', err=910) do 550 j = 1, 72 read (lucld, end=920, err=930) tmpcld do 500 k = 1, 144 cld2(k,j) = tmpcld(k) 500 continue 550 continue close(lucld) c return c c --- error found c 910 write(6,*) 'cannot open ',afile go to 990 920 write(6,*) 'unexpected end of ',afile go to 990 930 write(6,*) 'error in reading ',afile 990 lfatal = .true. c 999 return end c subroutine getsur (lusur,lucat, lfatal) c ********************************************************************** c c module name: getsur c c version : 1.0 c c programmer : lucy liu, hstx, 1996 c modified by: c c purpose : to read terrain height pressure data c to read surface category code c c method : c read 360 data reords (each having 720 elements) and store them in ps c read 360*6 data reords (each having 120 elements) and store them in cat c c calling routine: main c c calling sequence: getsur (lusur,lucat, lfatal) c c subroutines called: none c c common blocks used: surpre c c arguments and local variables: c c name type i/o descriptions c ------ ----- --- ----------------------------------------------- c lusur i*4 i logical unit number for reading terrain height c pressure data set/file c lucat i*4 i logical unit number for reading surface category code c data set/file c lfatal l*4 o flag of fatal error c*********************************************************************** c --- common areas ----------------------------------------------------- c real*4 ps(360,720) integer cat(360,720) common /surpre/ ps, cat c c --- arguments -------------------------------------------------------- c integer lusur, lucat c c --- local variables -------------------------------------------------- c logical lfatal character tmpch*120 c c --- read terrain height pressure data and store in array ps c lfatal = .false. read (lusur, end=920, err=940) ps c c --- read surface catgegory code and store in array cat c do 910 i=1,360 do 910 k=1,6 read (lucat, '(a120)', end=950, err=960) tmpch do 910 j=1,120 cat(i,(k-1)*120+j)=ichar(tmpch(j:j))-ichar('0') 910 continue return c c --- error found c 920 write(6,*) 'unexpected eof of terrain height pressure.' go to 990 940 write(6,*) 'error in reading terrain height pressure.' go to 990 950 write(6,*) 'unexpected eof of surface category code.' go to 990 960 write(6,*) 'error in reading surface category code.' 990 lfatal = .true. c 999 return end c subroutine incovr (idif,lat,lon, snoice,cldprs) c ********************************************************************** c c module name: incovr c c version : 1.0 c c programmer : david k. lee, stx (sasc), for nimbus-7, 1981 c rewritten : lucy liu, hstx 1996 c c purpose : compute the snow-ice thickness, covered on the earth, c in tenths of inches. c compute the cloud pressure. c c method : use a 1 x 1 degree snow-ice map for interpolation. c use a 2.5 x 2.5 degree cloud-pressure map for interpolation. c c calling routine: main c c calling sequence: incovr (idif,lat,lon, snoice,cldprs) c c subroutines called: none c c common blocks used: covpre c c arguments and local variables: c c name type i/o descriptions c ------ ----- --- ----------------------------------------------- c idif i*4 l difference c lat i*4 i latitude in radians * 10000 c lon i*4 i longitude in radians * 10000 c snoice i*4 o snow-ice depth at given latitude and longitude c cldprs i*4 o cloud-pressure at given latitude and longitude c raddeg r*4 l radians to degrees c ilat i*4 l c ilon i*4 l c flat r*4 l latitude in degrees c flon r*4 l longitude in degrees c ********************************************************************** c --- common areas ----------------------------------------------------- c integer monpre real map1(363,180), map2(363,180), cld1(144,72),cld2(144,72) common /covpre/ monpre, map1, map2, cld1, cld2 c c --- arguments -------------------------------------------------------- c integer idif, lat, lon integer snoice, cldprs c c --- local variables -------------------------------------------------- c real raddeg, flat, flon, grid integer ilat, ilon c data raddeg /57.295779/, grid/2.5/ c c --- fill data c snoice = -7777 cldprs = -7777 c ilat = nint (90. - float (lat) * raddeg / 1.e4) + 1 ilon = nint (amod (float (lon) * raddeg / 1.e4 + 360.,360.)) + 1 if (ilat.gt.180 ) ilat = 180 if (ilon.gt.360) ilon = 360 c if (map1(ilon+3,ilat).ge.0. .and. map2(ilon+3,ilat).ge.0.) then snoice = nint ((map1(ilon+3,ilat) + ((map2(ilon+3,ilat) - + map1(ilon+3,ilat)) / 30.5) * idif) / 10.) endif c flat = (-1.*float(lat)*raddeg/1.e4-grid/2.+90.) / grid if (lon.ge.0) then flon=(float(lon)*raddeg/1.e4-grid/2.) / grid else flon=(float(lon)*raddeg/1.e4-grid/2.+360.) / grid end if if (flat.le.0) then ilat = 0 else if (flat.ge.(180./grid-1)) then ilat = 180./grid-1 else ilat = nint(flat) end if if (flon.lt.0) flon = flon + 360. / grid ilon = nint(flon) ilat = ilat + 1 ilon = ilon + 1 if (ilat.gt.72) ilat=72 if (ilon.gt.144) ilon = 144 cldprs = nint (cld1(ilon,ilat) + (cld2(ilon,ilat) - + cld1(ilon,ilat)) / 30.5 * idif) return end c subroutine insurf (lat,lon, terpre,catcod) c ********************************************************************** c c module name: insurf c c version : 1.0 c c programmer : david k. lee, stx (sasc), for nimbus-7, 1981 c modified by: lucy liu, hstx, 1996 c c purpose : compute the terrain height pressure at a given c latitude and longitude on the surface of the earth. c find surface category code c c method : use a bilinear interpolation between the four c nearest grid points. c c calling routine: main c c calling sequence: insurf (lat,lon, terpre,catcod) c c subroutines called: none c c common blocks used: surpre c c arguments and local variables: c c name type i/o descriptions c ------ ----- --- ----------------------------------------------- c lat i*4 i latitude in radians * 10000 c lon i*4 i longitude in radians * 10000 c terpre i*4 o terrain height pressure at given c latitude and longitude c catcod i*4 o surface category code at given c latitude and longitude c raddeg r*4 l radians to degrees c flat r*4 l latitude in degree c flon r*4 l latitude in degree c ps1 r*4 l c ps2 r*4 l c ********************************************************************** c --- common areas ----------------------------------------------------- real*4 ps(360,720) integer cat(360,720) common /surpre/ ps, cat c c --- arguments -------------------------------------------------------- c integer lat,lon integer terpre,catcod c c --- local variables -------------------------------------------------- c real raddeg, epsil, grid real flat, flon, ps1, ps2 integer lat1, lat2, lon1, lon2 c data raddeg /57.295779/, epsil/1.e-4/, grid/0.5/ c c --- fill data c terpre = -7777 catcod = 9 c flat = (-1.*float(lat)*raddeg/1.e4-grid/2.+90.) / grid if (lon.ge.0) then flon=(float(lon)*raddeg/1.e4-grid/2.) / grid else flon=(float(lon)*raddeg/1.e4-grid/2.+360.) / grid end if if (flat.le.0) then flat = 0. lat1 = int(flat) lat2 = lat1 else if (flat.ge.(180./grid-1)) then flat = 180./grid-1 lat1 = int(flat) lat2 = lat1 else lat1 = int(flat) lat2 = lat1 + 1 end if if (flon.lt.0) flon = flon + 360. / grid lon1 = int(flon) lon2 = lon1 + 1 if (lon1.ge.int(360./grid)) lon1 = lon1-int(360./grid) if (lon2.ge.int(360./grid)) lon2 = lon2-int(360./grid) flat = flat + 1 flon = flon + 1 lat1 = lat1 + 1 lon1 = lon1 + 1 lat2 = lat2 + 1 lon2 = lon2 + 1 ps1 = ps(lat1,lon1) + (flat-lat1) * (ps(lat2,lon1)-ps(lat1,lon1)) ps2 = ps(lat1,lon2) + (flat-lat1) * (ps(lat2,lon2)-ps(lat1,lon2)) terpre = nint(ps1 + (flon-lon1) * (ps2 - ps1)) c if (abs(flat-lat1).le.epsil .and. abs(flon-lon1).le.epsil .or. + abs(flat-lat1).le.epsil .and. abs(flon-lon2).le.epsil .or. + abs(flat-lat2).le.epsil .and. abs(flon-lon1).le.epsil .or. + abs(flat-lat2).le.epsil .and. abs(flon-lon2).le.epsil) then catcod = cat(int(flat),int(flon)) else if (abs(flat-lat1).le.epsil) then if (cat(lat1,lon1).eq.cat(lat1,lon2)) then catcod = cat(lat1,lon1) else if (cat(lat1,lon1).eq.0 .and. cat(lat1,lon2).eq.1 .or. + cat(lat1,lon1).eq.1 .and. cat(lat1,lon2).eq.0) then catcod = 3 else if (cat(lat1,lon1).eq.1 .and. cat(lat1,lon2).eq.2 .or. + cat(lat1,lon1).eq.2 .and. cat(lat1,lon2).eq.1) then catcod = 4 end if else if (abs(flat-lat2).le.epsil) then if (cat(lat2,lon1).eq.cat(lat2,lon2)) then catcod = cat(lat2,lon1) else if (cat(lat2,lon1).eq.0 .and. cat(lat2,lon2).eq.1 .or. + cat(lat2,lon1).eq.1 .and. cat(lat2,lon2).eq.0) then catcod = 3 else if (cat(lat2,lon1).eq.1 .and. cat(lat2,lon2).eq.2 .or. + cat(lat2,lon1).eq.2 .and. cat(lat2,lon2).eq.1) then catcod = 4 end if else if (abs(flon-lon1).le.epsil) then if (cat(lat1,lon1).eq.cat(lat2,lon1)) then catcod = cat(lat1,lon1) else if (cat(lat1,lon1).eq.0 .and. cat(lat2,lon1).eq.1 .or. + cat(lat1,lon1).eq.1 .and. cat(lat2,lon1).eq.0) then catcod = 3 else if (cat(lat1,lon1).eq.1 .and. cat(lat2,lon1).eq.2 .or. + cat(lat1,lon1).eq.2 .and. cat(lat2,lon1).eq.1) then catcod = 4 end if else if (abs(flon-lon2).le.epsil) then if (cat(lat1,lon2).eq.cat(lat2,lon2)) then catcod = cat(lat1,lon2) else if (cat(lat1,lon2).eq.0 .and. cat(lat2,lon2).eq.1 .or. + cat(lat1,lon2).eq.1 .and. cat(lat2,lon2).eq.0) then catcod = 3 else if (cat(lat1,lon2).eq.1 .and. cat(lat2,lon2).eq.2 .or. + cat(lat1,lon2).eq.2 .and. cat(lat2,lon2).eq.1) then catcod = 4 end if else if (cat(lat1,lon1).eq.cat(lat1,lon2) .and. + cat(lat1,lon1).eq.cat(lat2,lon1) .and. + cat(lat1,lon1).eq.cat(lat2,lon2)) then catcod = cat(lat1,lon1) else if (cat(lat1,lon1).ne.2 .and. cat(lat1,lon2).ne.2 .and. + cat(lat2,lon1).ne.2 .and. cat(lat2,lon2).ne.2) then catcod = 3 else if (cat(lat1,lon1).ne.0 .and. cat(lat1,lon2).ne.0 .and. + cat(lat2,lon1).ne.0 .and. cat(lat2,lon2).ne.0) then catcod = 4 else write(6,*)'no surface category code at lat=',lat,' lon=',lon end if return end