subroutine snohires (si,sm,glat1,glon1,list,lusaf,im,jm,nm) ! ! inputs: sm (eta land/sea mask) ! glat (latitude array of eta grid) ! glon (longitude array of eta grid) ! list (unit number of printout) ! lusaf (logical - true:use ims snow and usaf snow, ! false: use ims snow only) ! dtr (degrees to radians conversion factor) ! ! outputs: si (snowdepth and sea-ice on eta grid, see calling ! subroutine cnsts for more details on si output) ! implicit real (a-h, o-z) ! real,intent(out) :: si(0:im+1,0:jm+1,nm) real,intent(in) :: sm(0:im+1,0:jm+1,nm) real,intent(in) :: glat1(0:im+1,0:jm+1,nm), glon1(0:im+1,0:jm+1,nm) logical,intent(inout) :: lusaf integer,intent(in) :: im,jm,nm ! integer::i,j,n ! real afsno (512,512) integer mskaf (512,512) ! real scvh (1024,1024) integer mskscvh (1024,1024) real,parameter::dtr=3.1415926/180. ! real::glat(0:im+1,0:jm+1,nm), glon(0:im+1,0:jm+1,nm) ! ! input units (45-nesdis daily snow/ice, ! 41-usaf daily snow/ice, ! 43-nesdis 1/16-bedient land/sea mask, ! 42-usaf 1/8-bedient land/sea mask) ! data insnoims/45/, insnoaf/41/, inmskaf/42/, inmskims/43/ ! ! set lower threshold for snodepth in tenths of inches (mininum ! snow depth is 1.5" where there is snow coverage). ! data idepth /15/ ! !************************** begin execution *********************** ! glon=-glon1 glat=glat1 ! specify the unit number of the printer. ioutuprt = list ! ! specify parameters of the n.h. 1024x1024 ims 1/16-mesh grid ! ! location of the pole: xpnmc8 = 513. ypnmc8 = 512. ! n -- the number of grid intervals from the pole to the equator: ennmc8 = 16.0e0 * 31.2043316e0 ! the longitudinal rotation angle of the grid alnmc8 = 10.e0 ! the orientation west longitude of the grid: orient8 = 80.e0 ! ! specify parameters of the n.h. 512 x 512 grid type 88 (usaf grid). ! ! location of pole: xpnmcaf = 257.e0 ypnmcaf = 256.e0 ! grid mesh length at 60n in km (fyi, but not needed) !.... xmeshl = 47.625e0 ! number of grid intervals from pole to equator ennmcaf = 8.0 * 31.2043316e0 ! the longitude rotation angle of the grid alnmcaf = 10.e0 ! the orientation west longitude of the grid orientaf = 80.e0 ! ! call to read nesdis/ims 1/16-bedient daily n.h. snow/ice via sno16get ! write(6,*) 'calling sno16get' call sno16get(scvh,iyear,imonth,iday,insnoims,list) write(6,*) 'back from sno16get' !test !test write(60,7701) iyear, imonth, iday !test 7701 format('ims data : ', 3i2.2) !test ! ! at present time the ims file has a 2-digit year. add century: ! now iyear from sno16get is 4-digits so we don't need to do this ! if (iyear .gt. 90) then ! iyear = iyear + 1900 ! else ! iyear = iyear + 2000 ! endif ! ! calculate julian date: nday = iw3jdn(iyear,imonth,iday) write (ioutuprt, 1) insnoims, iyear, imonth, iday, nday 1 format(1h ,' daily nesdis snow read in via unit no.=',i3/ & 1h ,' for valid date of y-m-d =',i4,'-',i2,'-',i2/ & 1h ,' which is julian day =', i10) ! ! read the nesdis-ims land-water mask ! (sea=0,land=1) do 10 j = 1, 1024 read(inmskims,'(80i1)') (mskscvh(i,j),i=1,1024) 10 continue ! call printims (scvh, mskscvh) ! ! print specified usaf snow processing flag write(ioutuprt,2211) lusaf 2211 format(//1h ,5x,'subroutine snohires '/ & 1h ,1x,' will try to process usaf snowcover: lusaf = ',l2) ! if (lusaf) then ! ! read air force 1/8-bedient daily n.h. snow/ice via sno8get ! write(6,*) 'calling sno8get' call sno8get(afsno,iyear,imonth,iday,insnoaf,list,lusaf) !test write(60,7702) iyear, imonth, iday 7702 format('usaf data: ', 3i2.2) !test ! ! if i/o error encountered in sno8get, then lusaf is returned false ! if (.not. lusaf) then write(ioutuprt,56) 56 format(1h ,'warning: file err in usaf snow anal') else ! the usaf snow header has a 2-digit year. add century: ! now iyear from sno8get is 4-digits so we don't need to do this ! if (iyear .gt. 90) then ! iyear = iyear + 1900 ! else ! iyear = iyear + 2000 ! endif ! ! calculate julian date: nafday = iw3jdn(iyear,imonth,iday) ! write(ioutuprt,2) insnoaf, iyear, imonth, iday, nafday 2 format(1h ,' daily usaf snow read in via unit no.=',i3/ & 1h ,' for valid date of y-m-d ='i4,'-',i2,'-',i2/ & 1h ,' which is julian day =', i10) ! if (nafday .lt. nday-4) then lusaf = .false. write(ioutuprt,2644) nday,nafday 2644 format(/1h ,'****** warning ****** warning ******'/ & 1h ,'date of usaf snow anal more than 4 days old'/ & 1h ,' - will fall-back to nesdis ims snow cover'/ & 1h ,43x,' nesdis snow anal julday=',i7/ & 1h ,43x,' usaf snow anal julday=',i7) ! endif endif endif ! ! note: upon return from call sno8get above, usaf snow/ice field ! has following physical ranges: ! - values over sea points are 0.0 or 11.0 (sea-ice flag) ! - values over land/coast are 0.0 or pos depth in meters ! - snowdepth over land is actual, not water equivalent yet ! ! read the usaf afgwc land/coast/sea mask ! (sea=1,land=2,coastal-land=4,offworld=9) read(inmskaf) mskaf ! !---------- i/o of primary input fields is complete ------------- ! init radians to degrees rtd = 1./dtr ! ! idepth is criterion of snowdepth threshold in tenths of ! inches below which usaf snowdepth will be assumed zero. ! here convert idepth to meters ! depth = float(idepth) * 2.54e-3 ! if ( lusaf) then write(ioutuprt,2321) lusaf,idepth 2321 format(1h //' usaf snodep anal will be used, lusaf=',l2/ & 1h , 35x,'snodepth threshold (tenths of inches) =',i3) call printaf (afsno,mskaf) else write(ioutuprt,2322) lusaf 2322 format(1h //' usaf snodep analysis will be ignored'/ & 1h , 35x,'logical flag lusaf=',l2) endif ! !----------initialize snow/ice arrays to zero on eta grid------------ ! si = 0.0 !-------------------------------------------------------------------- ! ! ****** now begin major loop over all eta grids and points ******* ! do n=1,nm do j=1,jm do i=1,im ! !--------------- determine lat/lon of eta grid point ------------- ! (here long will be east long) ! yylat = glat(i,j,n)*rtd xlong = 360. - glon(i,j,n)*rtd !!!!!!TEST!!!!!! ! if (xlong.gt.360)then ! print *,i,j,n,xlong,yylat ! endif !!!!!!!!!!!!!TEST!!!!!!!!!!!! ! where eta domain south of 22 n lat (including any s.h.), ! we keep default zero snow/ice ! ! print *,i,j,n,xlong,yylat,glon(i,j,n),glat(i,j,n),glon1(i,j,n) & ! ,glat1(i,j,n) if (yylat.lt.22.0e0) go to 4300 !---------------------------------------------------------------- ! ! determine location of eta point on the 1024 x 1024 nesdis/ims grid ! rm= ennmc8*cos(yylat * dtr) / (1.e0 + sin(yylat * dtr)) rad = (xlong - alnmc8) * dtr x = xpnmc8 + rm * cos(rad) y = ypnmc8 + rm * sin(rad) ! is = int(x) ip1 = is + 1 js = int(y) jp1 = js + 1 ! !-------------------------------------------------------------------- !.......first interpolate nesdis snow/ice as primary default...... ! ! the value of snow cover or sea ice on the 1024 x 1024 ! nesdis/ims grid is 1 for snow/ice points, 0 for snow/ice-free. ! we utilize an 1024 x 1024 land mask ! from ims to distinguish sea ice from snow points ! (actually, the data we got from ims have different values for ! snow and for ice. but since we want to use as much of the original ! program [written for the 1/2-mesh sab snow/ice] as possible, we give ! snow cover and sea ice the same value (=1) and let the mask do the ! job. ! ! now use eta and ims land-sea masks to ensure only land ! points are interpolated to land points (to determine snow) ! and only sea points are interpolated to sea points (for ice) ! (nesdis/ims land mask: sea=0,land=1, while the eta mask is the ! other way round). ! iland = 1 if( sm(i,j,n) .gt. 0.9 ) iland=0 ! ipoint = nint(x) jpoint = nint(y) ! print *,i,j,n,ipoint,jpoint if ( mskscvh(ipoint,jpoint) .eq. iland ) then si(i,j,n) = scvh(ipoint,jpoint) go to 3351 endif ! ! nearest neighbor not same sfc type, so use all 4 surrounding points ! kount = 0 ! xratio = x - real(is) yratio = y - real(js) ! area11 = (1.0e0 - xratio) * (1.0e0 - yratio) area21 = xratio * (1.0e0 - yratio) area12 = (1.0e0 - xratio) * yratio area22 = xratio * yratio ! if( mskscvh(is, js) .eq. iland) then kount = kount + 1 area = area11 ipoint = is jpoint = js end if ! if( mskscvh(is, jp1) .eq. iland ) then kount = kount +1 if (kount .eq. 1) then ipoint = is jpoint = jp1 elseif (area12 .gt. area) then area = area12 ipoint = is jpoint = jp1 end if end if ! if( mskscvh(ip1, js) .eq. iland ) then kount = kount + 1 if (kount .eq. 1) then area = area21 ipoint = ip1 jpoint = js elseif (area21 .gt. area) then area = area21 ipoint = ip1 jpoint = js end if end if ! if( mskscvh(ip1, jp1) .eq. iland ) then kount = kount + 1 if (kount .eq. 1) then area = area22 ipoint = ip1 jpoint = jp1 elseif (area22 .gt. area) then area = area22 ipoint = ip1 jpoint = jp1 end if end if ! ! determine sno/ice using nearest neighbor with same sfc type ! if(kount .gt. 0) then si(i,j,n) = scvh(ipoint,jpoint) else ! ! no immediately surrounding points in the 1024 x 1024 field of ! snow/ice have the same land-sea type as the eta point. the ! eta point may be small island or lake or small bay or pennin. ! (invariably a small lake in eta grid) ! so expand search radius and take first sfc type match ! ipoint = nint(x) jpoint = nint(y) ! !!! do 3346 ll=1,4 do 3346 ll=1,40 !DRAGAN jpe = min (1024, jpoint+ll) jpb = max (1 , jpoint-ll) ipe = min (1024, ipoint+ll) ipb = max (1 , ipoint-ll) ! do 2346 mk=jpb,jpe do 2346 nk=ipb,ipe if (mskscvh(nk,mk) .eq. iland) then si(i,j,n) = scvh(nk,mk) go to 3351 endif 2346 continue 3346 continue ! ! no land/sea mask matches found, so ! a) north of 55n, we assign snow/ice irrespective of sfc type ! b) south of 55n, we keep a priori zero default ! (the "b" option best for warmer lats of u.s., where this condition ! is virtually always a small eta lake with no counterpart water ! nearby in the nesdis/ims grid, e.g., salt lake, where we must ! avoid getting sea-ice owing to surrounding snow cover) ! write (ioutuprt, 3347) i,j,yylat,xlong,ipoint,jpoint,iland 3347 format(1h ,'**no ims msk match,', & ! 'eta-i,j,elat,elon,ims-i,ims-j,lnd:',2i6,2f7.2,2i3,i2) 'eta-i,j,elat,elon,ims-i,ims-j,lnd:',2i6,2f7.2,2i5,i2) if (yylat .ge. 55.0 ) then si(i,j,n) = scvh(ipoint,jpoint) else si(i,j,n) = 0.0 endif endif ! 3351 continue ! if (.not. lusaf) go to 4300 ! ! !-------------- begin usaf snow/ice interpolation------------------ ! if current usaf snow/ice anal was successfully read (lusaf=true), ! add the usaf information as follows: ! ! the quality of nesdis/ims sea-ice cover and snow coverage (especially ! in areas with small amounts of snow) are better than the usaf data. ! we do not use usaf sea-ice cover data at all. over land, the ! presence/absence of snow is determined by the ims snow coverage, i.e., ! at a given location, if the ims data indicate no snow, then we assume ! there is no snow, no matter what the usaf data say. if the ims data ! show snow but the usaf data have no snow or less than 1.5" of snow at ! the location, we assume there is a 1.5" of snow. if both the ims data ! and usaf data indicate snow, and the usaf snow depth is more than 1.5", ! then we use the usaf snow depth. ! ! on the 512 x 512 usaf grid, the data values are 0.0 for ! no snow or ice, 11.0 for ice points, and positive but ! less than 11.0 for snow points. ! !--------------------------------------if eta sea point, skip usaf anl if ( sm(i,j,n) .gt. 0.9 ) go to 4300 ! ! if the ims data indicate no snow on this point, also skip usaf anl: if (si(i,j,n).lt.0.001) go to 4300 ! ! otherwise, set snow depth to be the lower threshold, 1.5": si(i,j,n) = depth ! !------------------------------------------------------------------- ! this is an eta land point, so apply usaf snowdepth anal ! ! determine location of eta point on the 512 x 512 usaf grid ! rm= ennmcaf*cos(yylat * dtr) / (1.e0 + sin(yylat * dtr)) rad = (xlong - alnmcaf) * dtr x = xpnmcaf + rm * cos(rad) y = ypnmcaf + rm * sin(rad) ! is = int(x) ip1 = is + 1 js = int(y) jp1 = js + 1 ! !-----if outside of usaf grid domain (i.e. s.h.) we keep zero default-- ! if ((is .lt. 1) .or. (is .gt. 511) .or. (js .lt. 1) & .or. (js .gt. 511)) then go to 4300 endif !-------------------------------------------------------------------- ! ! now use eta and usaf land-sea mask to ensure only land points are ! interpolated to land points (to determine snow) ! (usaf land mask: sea=1,land=2,coastal-land=4,offworld=9) ! note: in reaching this stage, we have already insured we are on ! eta land point and in n.h., i.e. not mskaf=9 (i.e. not offworld) ! iland = 2 ! ipoint = nint(x) jpoint = nint(y) if ( mskaf(ipoint,jpoint) .ge. iland) then si(i,j,n) = afsno(ipoint,jpoint) si(i,j,n) = amax1(depth,si(i,j,n)) go to 4300 endif ! kount = 0 ! xratio = x - real(is) yratio = y - real(js) ! area11 = (1.0e0 - xratio) * (1.0e0 - yratio) area21 = xratio * (1.0e0 - yratio) area12 = (1.0e0 - xratio) * yratio area22 = xratio * yratio ! if( mskaf(is, js) .ge. iland) then kount = kount + 1 area = area11 ipoint = is jpoint = js end if ! if( mskaf(is, jp1) .ge. iland ) then kount = kount +1 if (kount .eq. 1) then ipoint = is jpoint = jp1 elseif (area12 .gt. area) then area = area12 ipoint = is jpoint = jp1 end if end if ! if( mskaf(ip1, js) .ge. iland ) then kount = kount + 1 if (kount .eq. 1) then area = area21 ipoint = ip1 jpoint = js elseif (area21 .gt. area) then area = area21 ipoint = ip1 jpoint = js end if end if ! if( mskaf(ip1, jp1) .ge. iland ) then kount = kount + 1 if (kount .eq. 1) then area = area22 ipoint = ip1 jpoint = jp1 elseif (area22 .gt. area) then area = area22 ipoint = ip1 jpoint = jp1 end if end if ! ! determine sno/ice considering the number of points surrounding ! eta grid point with the same land-sea mask flag ! if (kount .gt. 0) then si(i,j,n) = afsno(ipoint,jpoint) si(i,j,n) = amax1(depth,si(i,j,n)) ! else ! ! no immediately surrounding points in the 512 x 512 field of ! sno/ice have the same land-sea mask as the eta point. ! the eta point may be an island or lake or bay or penninsula. ! so expand search radius and take first mask flag match ! ipoint = nint(x) jpoint = nint(y) ! !!!!! do 7346 ll=1,7 do 7346 ll=1,70 !DRAGAN jpe = min (512, jpoint+ll) jpb = max (1 , jpoint-ll) ipe = min (512, ipoint+ll) ipb = max (1 , ipoint-ll) ! do 6346 mk=jpb,jpe do 6346 nk=ipb,ipe if (mskaf(nk,mk) .ge. iland) then si(i,j,n) = afsno(nk,mk) si(i,j,n) = amax1(depth,si(i,j,n)) go to 4300 endif 6346 continue 7346 continue ! ! no land mask match found, so we print warning and stay ! with earlier value determined from nesdis/ims anal, which ! we convert from cover flag (0,1) to default depth of 1.5". ! si(i,j,n) = depth write (ioutuprt, 7347) i,j,yylat,xlong,is,js,si(i,j,n) 7347 format(1h ,'*** warning ***..no usaf land msk match ', & ' at eta-i,j,elat,elon,usaf-i,j:',2i6,2f7.2,2i3/ & 1h ,' default to nesdis value of ',f5.1,' times .10 m') ! endif ! ! !******************** end major eta grid point loop **************** ! 4300 continue enddo enddo enddo ! ! print sample of snow/ice on eta grid ! call printeta (si,sm,im,jm,nm) !test call printyl(si,sm) return end