subroutine inter(yrdoy,flat,flon,pteran,pcloud,icat,isno) c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c implicit none c save c logical frscal/.true./, leap, lfatal integer lusur, lucat, bindoy(12), yr, doy, yrdoy, icat, isno, 1 iyr, mon, month1, month2, monpre, lat, lon, idif, iter, icld real flat, flon, pteran, pcloud character*60 terprs, surfcat character*60 afile integer bindoy1(12),bindoy2(12) include 'flnams.com' equivalence (fname(17),terprs),(fname(18),surfcat) data bindoy1 /15,46,74,105,135,166,196,227,258,288,319,349/ data bindoy2 /15,46,75,106,136,167,197,228,259,289,320,350/ 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=terprs 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=surfcat 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 call getsur(lusur,lucat, .false.) c frscal=.false. c endif c yr = int(yrdoy/1000) doy = mod(yrdoy,1000) c write (6,*) yrdoy,yr,doy if (leap (yr )) then do mon = 1, 12 bindoy(mon) = bindoy1(mon) enddo else do mon = 1, 12 bindoy(mon) = bindoy2(mon) enddo 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 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,flat,flon, isno,icld) pteran=iter/1013. pcloud=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(360,180), map2(360,180), cld1(360,180),cld2(360,180) common /covpre/ monpre, map1, map2, cld1, cld2 c c --- arguments -------------------------------------------------------- c integer month1, month2 logical lfatal c c --- local variables -------------------------------------------------- c real*4 tmpsno(18), tmpcld(18) integer lucov character mname(12)*3, afile*60 include 'flnams.com' character*60 snowice, cldpres equivalence (fname(19),snowice),(fname(20),cldpres) c data mname /'jan','feb','mar','apr','may','jun', + 'jul','aug','sep','oct','nov','dec'/ c lfatal = .false. c afile = snowice(1:index(snowice,' ')-1)//'.'//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 ' v8snowice.'//mname(month1) stop endif open (lucov,file=afile, status='old', + err=910) do 250 j = 1, 180 do 225 k=1,20 read (lucov, '(18f8.2)', end=920, err=930) tmpsno do 200 m = 1, 18 n = (k-1)*18 + m map1(n,j) = tmpsno(m) 200 continue 225 continue 250 continue close(lucov) c c --- read previous month data c afile = snowice(1:index(snowice,' ')-1)//'.'//mname(month2) call get_lun(lucov) if (lucov .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' v8snowice.'//mname(month2) stop endif open (lucov, file=afile, status='old', + err=910) do 350 j = 1, 180 do 325 k=1,20 read (lucov, '(18f8.2)', end=920, err=930) tmpsno do 300 m = 1, 18 n = (k-1)*18 + m map2(n,j) = tmpsno(m) 300 continue 325 continue 350 continue close(lucov) c afile = cldpres(1:index(cldpres,' ')-1)//'.'//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 ' v8cldpres.'//mname(month1) stop endif open (lucld,file=afile, status='old', + err=910) do 450 j = 1, 180 do 425 k=1,20 read (lucld, '(18f8.2)', end=920, err=930) tmpcld do 400 m = 1, 18 n = (k-1)*18 + m cld1(n,j) = tmpcld(m) 400 continue 425 continue 450 continue close(lucld) c c --- read previous month data c afile = cldpres(1:index(cldpres,' ')-1)//'.'//mname(month2) call get_lun(lucld) if (lucld .lt. 0) then print*,'error in getting a logical unit number for the', 1 ' v8cldpres.'//mname(month2) stop endif open (lucld, file=afile, status='old', + err=910) do 550 j = 1, 180 do 525 k=1,20 read (lucld, '(18f8.2)', end=920, err=930) tmpcld do 500 m = 1, 18 n = (k-1)*18 + m cld2(n,j) = tmpcld(m) 500 continue 525 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,flat,flon,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,flat,flon, 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 flat i*4 i latitude in degrees c flon i*4 i longitude in degrees 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 ilat i*4 l latitude index c ilon i*4 l longitude index c ********************************************************************** c --- common areas ----------------------------------------------------- c integer monpre real map1(360,180), map2(360,180), cld1(360,180),cld2(360,180) common /covpre/ monpre, map1, map2, cld1, cld2 c c --- arguments -------------------------------------------------------- c integer idif integer snoice, cldprs real flat, flon c c --- local variables -------------------------------------------------- c integer ilat, ilon c c --- fill data c snoice = -7777 cldprs = -7777 c ilat = nint (flat + 90.0) + 1 ilon = nint (flon + 180.0) + 1 if (ilat.gt.180 ) ilat = 180 if (ilon.gt.360) ilon = 360 c if (map1(ilon,ilat).ge.0. .and. map2(ilon,ilat).ge.0.) then snoice = nint ((map1(ilon,ilat) + ((map2(ilon,ilat) - + map1(ilon,ilat)) / 30.5) * idif) / 10.) endif c if (cld1(ilon,ilat).ge.0. .and. cld2(ilon,ilat).ge.0.) then cldprs = nint (cld1(ilon,ilat) + (cld2(ilon,ilat) - + cld1(ilon,ilat)) / 30.5 * idif) endif c 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