module deter_sfc_mod !$$$ module documentation block ! . . . . ! module: deter_sfc_mod ! prgmmr: lueken ! ! abstract: subroutine used to determine land surface type ! ! program history log: ! 2011-08-01 lueken - Moved all land surface type subroutines to new module ! 2013-01-23 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! ! subroutines included: ! sub deter_sfc ! sub deter_sfc_type ! sub deter_sfc2 ! sub deter_sfc_fov ! sub deter_sfc_amsre_low ! sub deter_sfc_gmi ! sub deter_zsfc_model ! sub reduce2full ! sub init_sfc ! sub time_int_sfc ! sub accum_sfc ! sub calc_sfc ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind use satthin, only: sno_full,isli_full,sst_full,soil_moi_full, & soil_temp_full,soil_type_full,veg_frac_full,veg_type_full, & fact10_full,zs_full,sfc_rough_full,zs_full_gfs use constants, only: zero,one,two,one_tenth,deg2rad,rad2deg use gridmod, only: nlat,nlon,regional,tll2xy,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc, & rlats,rlons,dx_gfs,txy2ll,lpl_gfs use guess_grids, only: nfldsfc,hrdifsfc,ntguessfc use calc_fov_crosstrk, only: npoly, fov_ellipse_crosstrk, inside_fov_crosstrk use calc_fov_conical, only: fov_ellipse_conical, inside_fov_conical implicit none ! Set default to private private ! Set passed variables to public public deter_sfc public deter_sfc_type public deter_sfc2 public deter_sfc_fov public deter_sfc_amsre_low public deter_sfc_gmi public deter_zsfc_model contains subroutine deter_sfc(alat,alon,dlat_earth,dlon_earth,obstime,isflg, & idomsfc,sfcpct,ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc determine land surface type ! prgmmr: derber org: np2 date: 2005-01-27 ! ! abstract: determines land surface type based on surrounding land ! surface types ! ! program history log: ! 2005-01-27 derber ! 2005-03-03 treadon - add implicit none, define zero ! 2006-02-01 parrish - change names of sno,isli,sst ! ! input argument list: ! alat ! alon ! obstime- observation time relative to analysis time ! dlat_earth ! dlon_earth ! ! output argument list: ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! sfcpct(0:3)- percentage of 4 surface types ! (0) - sea percentage ! (1) - land percentage ! (2) - sea ice percentage ! (3) - snow percentage ! tsavg - sea surface temperature ! idomsfc ! ts ! dfcr ! vty,vfr,sty,stp,sm,sn,zz,ff10 ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth,obstime,alat,alon integer(i_kind) ,intent( out) :: isflg,idomsfc real(r_kind),dimension(0:3),intent( out) :: sfcpct real(r_kind),dimension(0:3),intent( out) :: ts real(r_kind) ,intent( out) :: tsavg,sfcr real(r_kind) ,intent( out) :: vty,vfr,sty,stp,sm,sn,zz,ff10 real(r_kind),parameter:: minsnow=one_tenth integer(i_kind) istyp00,istyp01,istyp10,istyp11 integer(i_kind):: itsfc,itsfcp integer(i_kind):: ix,iy,ixp,iyp,j real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11,dtsfc,dtsfcp,wgtmin real(r_kind):: sno00,sno01,sno10,sno11,dlat,dlon real(r_kind):: sst00,sst01,sst10,sst11 real(r_kind),dimension(0:3)::wgtavg ! First do surface field since it is on model grid iy=int(alon); ix=int(alat) dy =alon-iy; dx =alat-ix dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=dx*dy ix=min(max(1,ix),nlat); iy=min(max(0,iy),nlon) ixp=min(nlat,ix+1); iyp=iy+1 if(iy==0) iy=nlon if(iyp==nlon+1) iyp=1 ! Interpolate fields which only vary in space (no time component) ! zz = surface height zz = zs_full(ix,iy) *w00 + zs_full(ixp,iy) *w10 + & zs_full(ix,iyp)*w01 + zs_full(ixp,iyp)*w11 if(regional)then ! call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) dlat=alat dlon=alon else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats_sfc,nlat_sfc,1) call grdcrd1(dlon,rlons_sfc,nlon_sfc,1) end if iy=int(dlon); ix=int(dlat) dy =dlon-iy; dx =dlat-ix dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=one-w00-w10-w01 ix=min(max(1,ix),nlat_sfc); iy=min(max(0,iy),nlon_sfc) ixp=min(nlat_sfc,ix+1); iyp=iy+1 if(iy==0) iy=nlon_sfc if(iyp==nlon_sfc+1) iyp=1 ! Get time interpolation factors for surface files if(obstime > hrdifsfc(1) .and. obstime <= hrdifsfc(nfldsfc))then do j=1,nfldsfc-1 if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then itsfc=j itsfcp=j+1 dtsfc=(hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j)) end if end do else if(obstime <=hrdifsfc(1))then itsfc=1 itsfcp=1 dtsfc=one else itsfc=nfldsfc itsfcp=nfldsfc dtsfc=one end if dtsfcp=one-dtsfc ! Set surface type flag. Begin by assuming obs over ice-free water istyp00 = isli_full(ix ,iy ) istyp10 = isli_full(ixp,iy ) istyp01 = isli_full(ix ,iyp) istyp11 = isli_full(ixp,iyp) sno00= sno_full(ix ,iy ,itsfc)*dtsfc+sno_full(ix ,iy ,itsfcp)*dtsfcp sno01= sno_full(ix ,iyp,itsfc)*dtsfc+sno_full(ix ,iyp,itsfcp)*dtsfcp sno10= sno_full(ixp,iy ,itsfc)*dtsfc+sno_full(ixp,iy ,itsfcp)*dtsfcp sno11= sno_full(ixp,iyp,itsfc)*dtsfc+sno_full(ixp,iyp,itsfcp)*dtsfcp sst00= sst_full(ix ,iy ,itsfc)*dtsfc+sst_full(ix ,iy ,itsfcp)*dtsfcp sst01= sst_full(ix ,iyp,itsfc)*dtsfc+sst_full(ix ,iyp,itsfcp)*dtsfcp sst10= sst_full(ixp,iy ,itsfc)*dtsfc+sst_full(ixp,iy ,itsfcp)*dtsfcp sst11= sst_full(ixp,iyp,itsfc)*dtsfc+sst_full(ixp,iyp,itsfcp)*dtsfcp ! Interpolate sst to obs location tsavg=sst00*w00+sst10*w10+sst01*w01+sst11*w11 if(istyp00 >=1 .and. sno00 > minsnow)istyp00 = 3 if(istyp01 >=1 .and. sno01 > minsnow)istyp01 = 3 if(istyp10 >=1 .and. sno10 > minsnow)istyp10 = 3 if(istyp11 >=1 .and. sno11 > minsnow)istyp11 = 3 sfcpct = zero sfcpct(istyp00)=sfcpct(istyp00)+w00 sfcpct(istyp01)=sfcpct(istyp01)+w01 sfcpct(istyp10)=sfcpct(istyp10)+w10 sfcpct(istyp11)=sfcpct(istyp11)+w11 isflg = 0 if(sfcpct(0) > 0.99_r_kind)then isflg = 0 else if(sfcpct(1) > 0.99_r_kind)then isflg = 1 else if(sfcpct(2) > 0.99_r_kind)then isflg = 2 else if(sfcpct(3) > 0.99_r_kind)then isflg = 3 else isflg = 4 end if ! vty = vegetation type ! sty = soil type ts(0:3)=zero wgtavg(0:3)=zero vfr=zero stp=zero sty=zero vty=zero sm=zero sn=zero idomsfc=isli_full(ix ,iy ) wgtmin = w00 if(istyp00 == 1)then vty = veg_type_full(ix ,iy) sty = soil_type_full(ix ,iy) wgtavg(1) = wgtavg(1) + w00 ts(1)=ts(1)+w00*sst00 vfr =vfr +w00*(veg_frac_full(ix ,iy ,itsfc ) *dtsfc+ & veg_frac_full(ix ,iy ,itsfcp) *dtsfcp) stp =stp +w00*(soil_temp_full(ix ,iy ,itsfc )*dtsfc+ & soil_temp_full(ix ,iy ,itsfcp)*dtsfcp) sm =sm +w00*(soil_moi_full(ix ,iy ,itsfc ) *dtsfc+ & soil_moi_full(ix ,iy ,itsfcp) *dtsfcp) else if(istyp00 == 2)then wgtavg(2) = wgtavg(2) + w00 ts(2)=ts(2)+w00*sst00 else if(istyp00 == 3)then wgtavg(3) = wgtavg(3) + w00 ts(3)=ts(3)+w00*sst00 sn = sn + w00*sno00 else wgtavg(0) = wgtavg(0) + w00 ts(0)=ts(0)+w00*sst00 end if if(istyp01 == 1)then if(wgtmin < w01 .or. (vty == zero .and. sty == zero))then vty = veg_type_full(ix ,iyp) sty = soil_type_full(ix ,iyp) end if wgtavg(1) = wgtavg(1) + w01 ts(1)=ts(1)+w01*sst01 vfr =vfr +w01*(veg_frac_full(ix ,iyp,itsfc ) *dtsfc+ & veg_frac_full(ix ,iyp,itsfcp) *dtsfcp) stp =stp +w01*(soil_temp_full(ix ,iyp,itsfc )*dtsfc+ & soil_temp_full(ix ,iyp,itsfcp)*dtsfcp) sm =sm +w01*(soil_moi_full(ix ,iyp,itsfc ) *dtsfc+ & soil_moi_full(ix ,iyp,itsfcp) *dtsfcp) else if(istyp01 == 2)then wgtavg(2) = wgtavg(2) + w01 ts(2)=ts(2)+w01*sst01 else if(istyp01 == 3)then wgtavg(3) = wgtavg(3) + w01 ts(3)=ts(3)+w01*sst01 sn = sn + w01*sno01 else wgtavg(0) = wgtavg(0) + w01 ts(0)=ts(0)+w01*sst01 end if if(wgtmin < w01)then idomsfc=isli_full(ix ,iyp) wgtmin = w01 end if if(istyp10 == 1)then if(wgtmin < w10 .or. (vty == zero .and. sty == zero))then vty = veg_type_full(ixp,iy) sty = soil_type_full(ixp,iy) end if wgtavg(1) = wgtavg(1) + w10 ts(1)=ts(1)+w10*sst10 vfr =vfr +w10*(veg_frac_full(ixp,iy ,itsfc ) *dtsfc+ & veg_frac_full(ixp,iy ,itsfcp) *dtsfcp) stp =stp +w10*(soil_temp_full(ixp,iy ,itsfc )*dtsfc+ & soil_temp_full(ixp,iy ,itsfcp)*dtsfcp) sm =sm +w10*(soil_moi_full(ixp,iy ,itsfc ) *dtsfc+ & soil_moi_full(ixp,iy ,itsfcp) *dtsfcp) else if(istyp10 == 2)then wgtavg(2) = wgtavg(2) + w10 ts(2)=ts(2)+w10*sst10 else if(istyp10 == 3)then wgtavg(3) = wgtavg(3) + w10 ts(3)=ts(3)+w10*sst10 sn = sn + w10*sno10 else wgtavg(0) = wgtavg(0) + w10 ts(0)=ts(0)+w10*sst10 end if if(wgtmin < w10)then idomsfc=isli_full(ixp,iy ) wgtmin = w10 end if if(istyp11 == 1)then if(wgtmin < w11 .or. (vty == zero .and. sty == zero))then vty = veg_type_full(ixp,iyp) sty = soil_type_full(ixp,iyp) endif wgtavg(1) = wgtavg(1) + w11 ts(1)=ts(1)+w11*sst11 vfr =vfr +w11*(veg_frac_full(ixp,iyp,itsfc ) *dtsfc+ & veg_frac_full(ixp,iyp,itsfcp) *dtsfcp) stp =stp +w11*(soil_temp_full(ixp,iyp,itsfc )*dtsfc+ & soil_temp_full(ixp,iyp,itsfcp)*dtsfcp) sm =sm +w11*(soil_moi_full(ixp,iyp,itsfc ) *dtsfc+ & soil_moi_full(ixp,iyp,itsfcp) *dtsfcp) else if(istyp11 == 2)then wgtavg(2) = wgtavg(2) + w11 ts(2)=ts(2)+w11*sst11 else if(istyp11 == 3)then wgtavg(3) = wgtavg(3) + w11 ts(3)=ts(3)+w11*sst11 sn = sn + w11*sno11 else wgtavg(0) = wgtavg(0) + w11 ts(0)=ts(0)+w11*sst11 end if if(wgtmin < w11)then idomsfc=isli_full(ixp,iyp) wgtmin = w11 end if if(wgtavg(0) > zero)then ts(0) = ts(0)/wgtavg(0) else ts(0) = tsavg end if if(wgtavg(1) > zero)then ts(1) = ts(1)/wgtavg(1) sm = sm/wgtavg(1) vfr = vfr/wgtavg(1) stp = stp/wgtavg(1) else ts(1) = tsavg sm=one end if if(wgtavg(2) > zero)then ts(2) = ts(2)/wgtavg(2) else ts(2) = tsavg end if if(wgtavg(3) > zero)then ts(3) = ts(3)/wgtavg(3) sn = sn/wgtavg(3) else ts(3) = tsavg end if ! ts(0)=max(ts(0),270._r_kind) ! ts(2)=min(ts(2),280._r_kind) ! ts(3)=min(ts(3),280._r_kind) ! Space-time interpolation of fields from surface wind speed ff10=(fact10_full(ix ,iy ,itsfc )*w00+ & fact10_full(ixp,iy ,itsfc )*w10+ & fact10_full(ix ,iyp,itsfc )*w01+ & fact10_full(ixp,iyp,itsfc )*w11)*dtsfc + & (fact10_full(ix ,iy ,itsfcp)*w00+ & fact10_full(ixp,iy ,itsfcp)*w10+ & fact10_full(ix ,iyp,itsfcp)*w01+ & fact10_full(ixp,iyp,itsfcp)*w11)*dtsfcp sfcr=(sfc_rough_full(ix ,iy ,itsfc )*w00+ & sfc_rough_full(ixp,iy ,itsfc )*w10+ & sfc_rough_full(ix ,iyp,itsfc )*w01+ & sfc_rough_full(ixp,iyp,itsfc )*w11)*dtsfc + & (sfc_rough_full(ix ,iy ,itsfcp)*w00+ & sfc_rough_full(ixp,iy ,itsfcp)*w10+ & sfc_rough_full(ix ,iyp,itsfcp)*w01+ & sfc_rough_full(ixp,iyp,itsfcp)*w11)*dtsfcp return end subroutine deter_sfc subroutine deter_sfc_type(dlat_earth,dlon_earth,obstime,isflg,tsavg) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc determine land surface type ! prgmmr: derber org: np2 date: 2005-01-27 ! ! abstract: determines land surface type based on surrounding land ! surface types ! ! program history log: ! 2005-01-27 derber ! 2005-03-03 treadon - add implicit none, define zero ! 2006-02-01 parrish - change names of sno,isli,sst ! ! input argument list: ! dlat ! dlon ! obstime ! ! output argument list: ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! tsavg ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth,obstime integer(i_kind),intent( out) :: isflg real(r_kind) ,intent( out) :: tsavg logical outside integer(i_kind) istyp00,istyp01,istyp10,istyp11 integer(i_kind):: ix,iy,ixp,iyp,j,itsfc,itsfcp real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11,dtsfc real(r_kind):: dtsfcp,dlat,dlon real(r_kind):: sst00,sst01,sst10,sst11 real(r_kind):: sno00,sno01,sno10,sno11 real(r_kind),parameter:: minsnow=one_tenth real(r_kind),dimension(0:3):: sfcpct if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats_sfc,nlat_sfc,1) call grdcrd1(dlon,rlons_sfc,nlon_sfc,1) end if iy=int(dlon); ix=int(dlat) dy =dlon-iy; dx =dlat-ix dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=one-w00-w10-w01 ix=min(max(1,ix),nlat_sfc); iy=min(max(0,iy),nlon_sfc) ixp=min(nlat_sfc,ix+1); iyp=iy+1 if(iy==0) iy=nlon_sfc if(iyp==nlon_sfc+1) iyp=1 ! Get time interpolation factors for surface files if(obstime > hrdifsfc(1) .and. obstime <= hrdifsfc(nfldsfc))then do j=1,nfldsfc-1 if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then itsfc=j itsfcp=j+1 dtsfc=(hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j)) end if end do else if(obstime <=hrdifsfc(1))then itsfc=1 itsfcp=1 dtsfc=one else itsfc=nfldsfc itsfcp=nfldsfc dtsfc=one end if dtsfcp=one-dtsfc ! Set surface type flag. Begin by assuming obs over ice-free water istyp00 = isli_full(ix ,iy ) istyp10 = isli_full(ixp,iy ) istyp01 = isli_full(ix ,iyp) istyp11 = isli_full(ixp,iyp) sno00= sno_full(ix ,iy ,itsfc)*dtsfc+sno_full(ix ,iy ,itsfcp)*dtsfcp sno01= sno_full(ix ,iyp,itsfc)*dtsfc+sno_full(ix ,iyp,itsfcp)*dtsfcp sno10= sno_full(ixp,iy ,itsfc)*dtsfc+sno_full(ixp,iy ,itsfcp)*dtsfcp sno11= sno_full(ixp,iyp,itsfc)*dtsfc+sno_full(ixp,iyp,itsfcp)*dtsfcp sst00= sst_full(ix ,iy ,itsfc)*dtsfc+sst_full(ix ,iy ,itsfcp)*dtsfcp sst01= sst_full(ix ,iyp,itsfc)*dtsfc+sst_full(ix ,iyp,itsfcp)*dtsfcp sst10= sst_full(ixp,iy ,itsfc)*dtsfc+sst_full(ixp,iy ,itsfcp)*dtsfcp sst11= sst_full(ixp,iyp,itsfc)*dtsfc+sst_full(ixp,iyp,itsfcp)*dtsfcp ! Interpolate sst to obs location tsavg=sst00*w00+sst10*w10+sst01*w01+sst11*w11 if(istyp00 >=1 .and. sno00 > minsnow)istyp00 = 3 if(istyp01 >=1 .and. sno01 > minsnow)istyp01 = 3 if(istyp10 >=1 .and. sno10 > minsnow)istyp10 = 3 if(istyp11 >=1 .and. sno11 > minsnow)istyp11 = 3 sfcpct = zero sfcpct(istyp00)=sfcpct(istyp00)+w00 sfcpct(istyp01)=sfcpct(istyp01)+w01 sfcpct(istyp10)=sfcpct(istyp10)+w10 sfcpct(istyp11)=sfcpct(istyp11)+w11 isflg = 0 if(sfcpct(0) > 0.99_r_kind)then isflg = 0 else if(sfcpct(1) > 0.99_r_kind)then isflg = 1 else if(sfcpct(2) > 0.99_r_kind)then isflg = 2 else if(sfcpct(3) > 0.99_r_kind)then isflg = 3 else isflg = 4 end if return end subroutine deter_sfc_type subroutine deter_sfc2(dlat_earth,dlon_earth,obstime,idomsfc,tsavg,ff10,sfcr,zz) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc determine land surface type ! prgmmr: derber org: np2 date: 2005-01-27 ! ! abstract: determines land surface type based on surrounding land ! surface types ! ! program history log: ! 2005-01-27 derber ! 2005-03-03 treadon - add implicit none, define zero ! 2006-02-01 parrish - change names of sno,isli,sst ! 2010-12-05 pondeca - add local terrain elevation zz as ! optional output variable ! ! input argument list: ! dlat_earth ! dlon_earth ! obstime- observation time relative to analysis time ! ! output argument list: ! tsavg - sea surface temperature ! idomsfc ! sfcr ! ff10 ! zz ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth,obstime integer(i_kind),intent( out) :: idomsfc real(r_kind) ,intent( out) :: tsavg,sfcr real(r_kind) ,intent( out) :: ff10 real(r_kind),optional,intent( out) :: zz integer(i_kind):: itsfc,itsfcp integer(i_kind):: ix,iy,ixp,iyp,j real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11,dtsfc,dtsfcp,wgtmin real(r_kind):: sst00,sst01,sst10,sst11,dlat,dlon logical outside if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats_sfc,nlat_sfc,1) call grdcrd1(dlon,rlons_sfc,nlon_sfc,1) end if iy=int(dlon); ix=int(dlat) dy =dlon-iy; dx =dlat-ix dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=dx*dy ix=min(max(1,ix),nlat_sfc); iy=min(max(0,iy),nlon_sfc) ixp=min(nlat_sfc,ix+1); iyp=iy+1 if(iy==0) iy=nlon_sfc if(iyp==nlon_sfc+1) iyp=1 ! Get time interpolation factors for surface files if(obstime > hrdifsfc(1) .and. obstime <= hrdifsfc(nfldsfc))then do j=1,nfldsfc-1 if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then itsfc=j itsfcp=j+1 dtsfc=(hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j)) end if end do else if(obstime <=hrdifsfc(1))then itsfc=1 itsfcp=1 dtsfc=one else itsfc=nfldsfc itsfcp=nfldsfc dtsfc=one end if dtsfcp=one-dtsfc sst00= sst_full(ix ,iy ,itsfc)*dtsfc+sst_full(ix ,iy ,itsfcp)*dtsfcp sst01= sst_full(ix ,iyp,itsfc)*dtsfc+sst_full(ix ,iyp,itsfcp)*dtsfcp sst10= sst_full(ixp,iy ,itsfc)*dtsfc+sst_full(ixp,iy ,itsfcp)*dtsfcp sst11= sst_full(ixp,iyp,itsfc)*dtsfc+sst_full(ixp,iyp,itsfcp)*dtsfcp ! Interpolate sst to obs location tsavg=sst00*w00+sst10*w10+sst01*w01+sst11*w11 idomsfc=isli_full(ix ,iy ) wgtmin = w00 if(wgtmin < w01 )then idomsfc=isli_full(ix ,iyp) wgtmin = w01 end if if(wgtmin < w10)then idomsfc=isli_full(ixp,iy ) wgtmin = w10 end if if(wgtmin < w11)then idomsfc=isli_full(ixp,iyp) wgtmin = w11 end if if((isli_full(ix ,iy ) /= isli_full(ix ,iyp)) .or. & (isli_full(ix ,iy ) /= isli_full(ixp,iy )) .or. & (isli_full(ix ,iy ) /= isli_full(ixp,iyp)) .or. & (isli_full(ixp,iy ) /= isli_full(ix ,iyp)) .or. & (isli_full(ixp,iy ) /= isli_full(ixp,iyp)) .or. & (isli_full(ix ,iyp) /= isli_full(ixp,iyp)) ) idomsfc = idomsfc+3 ! Space-time interpolation of fields from surface wind speed ff10=(fact10_full(ix ,iy ,itsfc )*w00+ & fact10_full(ixp,iy ,itsfc )*w10+ & fact10_full(ix ,iyp,itsfc )*w01+ & fact10_full(ixp,iyp,itsfc )*w11)*dtsfc + & (fact10_full(ix ,iy ,itsfcp)*w00+ & fact10_full(ixp,iy ,itsfcp)*w10+ & fact10_full(ix ,iyp,itsfcp)*w01+ & fact10_full(ixp,iyp,itsfcp)*w11)*dtsfcp sfcr=(sfc_rough_full(ix ,iy ,itsfc )*w00+ & sfc_rough_full(ixp,iy ,itsfc )*w10+ & sfc_rough_full(ix ,iyp,itsfc )*w01+ & sfc_rough_full(ixp,iyp,itsfc )*w11)*dtsfc + & (sfc_rough_full(ix ,iy ,itsfcp)*w00+ & sfc_rough_full(ixp,iy ,itsfcp)*w10+ & sfc_rough_full(ix ,iyp,itsfcp)*w01+ & sfc_rough_full(ixp,iyp,itsfcp)*w11)*dtsfcp if(present(zz)) then if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) end if iy=int(dlon); ix=int(dlat) dy =dlon-iy; dx =dlat-ix dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx*dy1; w01=dx1*dy; w11=dx*dy ix=min(max(1,ix),nlat); iy=min(max(0,iy),nlon) ixp=min(nlat,ix+1); iyp=iy+1 if(iy==0) iy=nlon if(iyp==nlon+1) iyp=1 zz = zs_full(ix,iy) *w00 + zs_full(ixp,iy) *w10 + & zs_full(ix,iyp)*w01 + zs_full(ixp,iyp)*w11 endif return end subroutine deter_sfc2 subroutine deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,& dlon_earth_deg,expansion,obstime,isflg,idomsfc, & sfcpct,vfr,sty,vty,stp,sm,ff10,sfcr,zz,sn,ts,tsavg) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc_fov determine surface characteristics ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: determines surface characteristics within a field of view ! based on model information and the size/shape of the ! field of view. ! ! program history log: ! 2008-11-04 gayno ! 2009-12-20 gayno - modify to use relative antenna power. ! ! input argument list: ! fov_flag - is this a crosstrack or conical instrument? ! ichan - channel number - conical scanners only ! ifov - field of view number ! instr - instrument number ! sat_aziang - satellite azimuth angle ! dlat_earth_deg - latitude of fov center (degrees) ! dlon_earth_deg - longitude of fov center (degrees) ! expansion - fov expansion factor ! obstime - observation time ! ! output argument list: ! idomsfc - dominate surface type ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! sfcpct(0:3)- percentage of 4 surface types ! (0) - sea percentage ! (1) - land percentage ! (2) - sea ice percentage ! (3) - snow percentage ! vfr - vegetation fraction ! sty - dominate soil type ! vty - dominate vegetation type ! stp - top layer soil temperature ! sm - top layer soil moisture ! ff10 - wind factor ! sfcr - surface roughness lenght ! zz - terrain height ! sn - snow depth ! ts - skin temperature for each surface type ! tsavg - average skin temperature ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none ! Declare passed variables. character(len=*), intent(in ) :: fov_flag integer(i_kind) , intent(in ) :: ifov, instr integer(i_kind) , intent(in ) :: ichan real(r_kind) , intent(in ) :: dlat_earth_deg, dlon_earth_deg, sat_aziang real(r_kind) , intent(in ) :: expansion, obstime integer(i_kind) , intent( out) :: isflg, idomsfc(1) real(r_kind) , intent( out) :: sfcpct(0:3), sty, vty, vfr, stp, sm real(r_kind) , intent( out) :: ff10, sfcr, zz, sn, ts(0:3), tsavg ! Declare local variables. integer(i_kind) :: i, ii, iii, j, jj, jjj integer(i_kind) :: is(npoly), js(npoly) integer(i_kind) :: n, nearest_i, nearest_j integer(i_kind), allocatable :: max_i(:), min_i(:) integer(i_kind) :: ifull, jstart, jend integer(i_kind) :: itsfc, itsfcp integer(i_kind) :: subgrid_lengths_x, subgrid_lengths_y logical :: outside real(r_kind) :: lats_edge_fov(npoly), lons_edge_fov(npoly) real(r_kind) :: lat_mdl, lon_mdl real(r_kind) :: lat_rad, lon_rad real(r_kind) :: dtsfc, dtsfcp real(r_kind) :: x, xstart, xend, y, ystart, yend real(r_kind) :: dx_fov, dx_fov_max, dy_fov, power real(r_kind) :: del, mid, rlon1, rlon2 real(r_kind), allocatable :: y_off(:), x_off(:), powerx(:,:) type surface2 sequence real(r_kind) :: vfr real(r_kind) :: sfcr real(r_kind) :: sm real(r_kind) :: stp real(r_kind) :: ff10 real(r_kind) :: sn real(r_kind) :: ts end type surface2 type(surface2) :: sfc_mdl ! holds time interpolated surface data type surface sequence real(r_kind) :: vfr real(r_kind) :: sm real(r_kind) :: stp real(r_kind) :: ff10 real(r_kind) :: sfcr real(r_kind) :: zz real(r_kind) :: sn real(r_kind), dimension(0:3) :: ts real(r_kind), dimension(0:3) :: count real(r_kind), dimension(0:24) :: count_vty real(r_kind), dimension(0:16) :: count_sty end type surface type(surface) :: sfc_sum ! holds sums during integration across fov ! Get time interpolation factors for surface files if(obstime > hrdifsfc(1) .and. obstime <= hrdifsfc(nfldsfc))then do j=1,nfldsfc-1 if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then itsfc=j itsfcp=j+1 dtsfc=(hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j)) end if end do else if(obstime <=hrdifsfc(1))then itsfc=1 itsfcp=1 dtsfc=one else itsfc=nfldsfc itsfcp=nfldsfc dtsfc=one end if dtsfcp=one-dtsfc ! The algorithm that locates the fov breaks down if its crosses the pole. ! so, just assign surface characteristics using a nearest neighbor ! approach. if (abs(dlat_earth_deg) > 88.0_r_kind) then ! print*,'USE NEAREST NEIGHBOR NEAR POLE' if (regional) then lat_rad = dlat_earth_deg*deg2rad lon_rad = dlon_earth_deg*deg2rad call tll2xy(lon_rad,lat_rad,x,y,outside) nearest_i = nint(x) nearest_j = nint(y) else y = dlat_earth_deg*deg2rad call grdcrd1(y,rlats_sfc,nlat_sfc,1) nearest_j = nint(y) jj = nearest_j if (jj > nlat_sfc/2) jj = nlat_sfc - nearest_j + 1 x = (dlon_earth_deg/dx_gfs(jj)) + one nearest_i = nint(x) call reduce2full(nearest_i,nearest_j,ifull) nearest_i = ifull end if call time_int_sfc(nearest_i,nearest_j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) call init_sfc(sfc_sum) power = one call accum_sfc(nearest_i,nearest_j,power,sfc_mdl,sfc_sum) call calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm,stp,ff10,sfcr,zz,sn,ts,tsavg) return endif ! Determine the edge of the fov. the fov is described by a polygon with ! with npoly-1 vertices. the routine returns the lat/lon of the intersection ! of the vertices. the size of the polygon is a function of the ! expansion factor. if (fov_flag=="crosstrk")then call fov_ellipse_crosstrk(ifov,instr,sat_aziang,dlat_earth_deg,dlon_earth_deg, & lats_edge_fov,lons_edge_fov) elseif(fov_flag=="conical")then call fov_ellipse_conical(ichan,sat_aziang,dlat_earth_deg,dlon_earth_deg, & lats_edge_fov,lons_edge_fov) endif if (regional) then xstart=999999._r_kind xend=-999999._r_kind ystart=999999._r_kind yend=-999999._r_kind do n = 1, npoly lat_rad = lats_edge_fov(n)*deg2rad lon_rad = lons_edge_fov(n)*deg2rad call tll2xy(lon_rad,lat_rad,x,y,outside) ! if any part of fov is outside grid, just set is/js to ! be the near grid point at fov center. we already know the ! center grid point is inside the grid as that is checked ! in calling routine. if (outside) then ! print*,"FOV ON EDGE OF GRID" lat_rad = dlat_earth_deg*deg2rad lon_rad = dlon_earth_deg*deg2rad call tll2xy(lon_rad,lat_rad,x,y,outside) js = nint(y) is = nint(x) exit endif xstart=min(xstart,x) xend=max(xend,x) ystart=min(ystart,y) yend=max(yend,y) is(n)=nint(x) js(n)=nint(y) enddo jstart=minval(js) jend=maxval(js) allocate (max_i(jstart:jend)) allocate (min_i(jstart:jend)) max_i = -999 min_i = 999999 do j = jstart, jend do n = 1, npoly if (js(n) == j) then max_i(j) = max(max_i(j),is(n)) min_i(j) = min(min_i(j),is(n)) endif enddo enddo else ! global grid ! Locate the fov on the model grid. in the "j" direction, this is ! based on the latitudinal extent of the fov. yend = maxval(lats_edge_fov)*deg2rad call grdcrd1(yend,rlats_sfc,nlat_sfc,1) ystart = minval(lats_edge_fov)*deg2rad call grdcrd1(ystart,rlats_sfc,nlat_sfc,1) ! Note two extra rows are added for the n/s poles. jstart = nint(ystart) jstart = max(jstart,2) jend = nint(yend) jend = min(jend,nlat_sfc-1) ! Locate the extent of the fov on the model grid in the "i" direction. note, the ! algorithm works on the reduced gaussian grid. hence, the starting/ending ! "i" indices are a function of "j". allocate (max_i(jstart:jend)) allocate (min_i(jstart:jend)) do j = jstart, jend jj = j if (jj > nlat_sfc/2) jj = nlat_sfc - j + 1 x = (minval(lons_edge_fov)/dx_gfs(jj)) + one nearest_i = nint(x) min_i(j) = nearest_i x = (maxval(lons_edge_fov)/dx_gfs(jj)) + one nearest_i = nint(x) max_i(j) = nearest_i enddo end if ! is this regional or global ! In this case, the fov is located completely within one grid ! point. this is common when the fov is small compared with ! the model grid resolution. if ((jstart == jend) .and. (max_i(jstart) == min_i(jstart))) then ! print*,'ONLY ONE MODEL POINT WITHIN FOV' nearest_j = jstart nearest_i = max_i(jstart) if (.not. regional) then call reduce2full(nearest_i,nearest_j,ifull) nearest_i = ifull endif call time_int_sfc(nearest_i,nearest_j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) call init_sfc(sfc_sum) power = one call accum_sfc(nearest_i,nearest_j,power,sfc_mdl,sfc_sum) call calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm,stp,ff10,sfcr,zz,sn,ts,tsavg) deallocate(max_i,min_i) return endif ! Find the size of the fov in model grid units. if the ! fov is small compared to the model grid, there is the ! possibility that there are no model grid boxes inside the fov. ! (a grid box is "inside" the fov if its center is within ! the fov.) this situation can be avoided by ! "chopping" the model grid boxes into smaller pieces, which ! is done below. first, find the n/s size of the fov ! in grid units. if (regional) then dy_fov = abs(yend-ystart) dx_fov = abs(xend-xstart) else ! global dy_fov = abs(yend-ystart) ! n/s size of fov in grid units. ! Find the e/w size of the fov in grid units. note: the resolution ! of the model decreases towards the poles. use the max ! value of model grid spacing in this calculation to make ! the most conservative estimate as to whether or not to ! "chop" the model grid boxes. dx_fov_max = zero do j = jstart,jend jj = j if (jj > nlat_sfc/2) jj = nlat_sfc - j + 1 dx_fov_max = max(dx_fov_max, dx_gfs(jj)) enddo ! When taking the longitudinal difference, don't worry ! about greenwich or the dateline as the fov code calculates ! longitude relative to the center of the fov. rlon1 = maxval(lons_edge_fov) rlon2 = minval(lons_edge_fov) dx_fov = (rlon1-rlon2)/dx_fov_max end if ! is this regional or global? ! if the fov is small compared to the model resolution, there ! is a possibility that there will be no model points located ! within the fov. to prevent this, the model grid may be ! subdivided into smaller pieces. this subdivision is ! done separately for each direction. subgrid_lengths_x = nint(one/dx_fov) + 1 subgrid_lengths_y = nint(one/dy_fov) + 1 loop1:do ! If the fov is very small compared to the model grid, it ! is more computationally efficient to take a simple average. if (subgrid_lengths_x > 7 .or. subgrid_lengths_y > 7) then ! print*,'FOV MUCH SMALLER THAN MODEL GRID POINTS, TAKE SIMPLE AVERAGE.' call init_sfc(sfc_sum) if (regional) then do j = jstart, jend do i = min_i(j), max_i(j) call time_int_sfc(i,j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) power = one call accum_sfc(i,j,power,sfc_mdl,sfc_sum) enddo enddo else ! global do j = jstart, jend do i = min_i(j), max_i(j) ii = i call reduce2full(ii,j,ifull) call time_int_sfc(ifull,j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) power = one call accum_sfc(ifull,j,power,sfc_mdl,sfc_sum) enddo enddo endif exit loop1 endif mid = (float(subgrid_lengths_y)-one)/two + one del = one/ float(subgrid_lengths_y) allocate (y_off(subgrid_lengths_y)) do i= 1, subgrid_lengths_y y_off(i) = (float(i)-mid)*del enddo mid = (float(subgrid_lengths_x)-one)/two + one del = one / float(subgrid_lengths_x) allocate (x_off(subgrid_lengths_x)) do i= 1, subgrid_lengths_x x_off(i) = (float(i)-mid)*del enddo ! Determine the surface characteristics by integrating over the ! fov. call init_sfc(sfc_sum) if (regional) then do j = jstart, jend do i = min_i(j), max_i(j) call time_int_sfc(i,j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) do jjj = 1, subgrid_lengths_y y = float(j) + y_off(jjj) do iii = 1, subgrid_lengths_x x = float(i) + x_off(iii) call txy2ll(x,y,lon_rad,lat_rad) lat_mdl = lat_rad*rad2deg lon_mdl = lon_rad*rad2deg if (lon_mdl < zero) lon_mdl = lon_mdl + 360._r_kind if (lon_mdl >= 360._r_kind) lon_mdl = lon_mdl - 360._r_kind if (fov_flag=="crosstrk")then call inside_fov_crosstrk(instr,ifov,sat_aziang, & dlat_earth_deg,dlon_earth_deg, & lat_mdl, lon_mdl, & expansion, ichan, power ) elseif (fov_flag=="conical")then call inside_fov_conical(instr,ichan,sat_aziang, & dlat_earth_deg,dlon_earth_deg,& lat_mdl, lon_mdl, & expansion, power ) endif call accum_sfc(i,j,power,sfc_mdl,sfc_sum) enddo enddo enddo enddo else allocate(powerx(subgrid_lengths_x,subgrid_lengths_y)) do j = jstart, jend jj = j if (j > nlat_sfc/2) jj = nlat_sfc - j + 1 do i = min_i(j), max_i(j) call reduce2full(i,j,ifull) call time_int_sfc(ifull,j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) !$omp parallel do schedule(dynamic,1)private(jjj,iii,lat_mdl,lon_mdl) do jjj = 1, subgrid_lengths_y if (y_off(jjj) >= zero) then lat_mdl = (one-y_off(jjj))*rlats_sfc(j)+y_off(jjj)*rlats_sfc(j+1) else lat_mdl = (one+y_off(jjj))*rlats_sfc(j)-y_off(jjj)*rlats_sfc(j-1) endif lat_mdl = lat_mdl * rad2deg do iii = 1, subgrid_lengths_x ! Note, near greenwich, "i" index may be out of range. that is ! ok here when calculating longitude even if the value is ! greater than 360. the ellipse code works from longitude relative ! to the center of the fov. lon_mdl = (float(i)+x_off(iii) - one) * dx_gfs(jj) if (fov_flag=="crosstrk")then call inside_fov_crosstrk(instr,ifov,sat_aziang, & dlat_earth_deg,dlon_earth_deg, & lat_mdl, lon_mdl, & expansion, ichan, powerx(iii,jjj) ) elseif (fov_flag=="conical")then call inside_fov_conical(instr,ichan,sat_aziang, & dlat_earth_deg,dlon_earth_deg,& lat_mdl, lon_mdl, & expansion, powerx(iii,jjj) ) endif enddo enddo do jjj = 1, subgrid_lengths_y do iii = 1, subgrid_lengths_x call accum_sfc(ifull,j,powerx(iii,jjj),sfc_mdl,sfc_sum) enddo enddo enddo enddo deallocate(powerx) endif ! regional or global deallocate (x_off, y_off) ! If there were no model points within the fov, the model points need to be ! "chopped" into smaller pieces. if (sum(sfc_sum%count) == zero) then close(9) subgrid_lengths_x = subgrid_lengths_x + 1 subgrid_lengths_y = subgrid_lengths_y + 1 ! print*,'NO GRID POINTS INSIDE FOV, CHOP MODEL BOX INTO FINER PIECES',subgrid_lengths_x,subgrid_lengths_y else exit loop1 endif end do loop1 deallocate (max_i, min_i) call calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm,stp,ff10,sfcr,zz,sn,ts,tsavg) return end subroutine deter_sfc_fov subroutine deter_sfc_amsre_low(dlat_earth,dlon_earth,isflg,sfcpct) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc_amsre_low determine land surface type ! prgmmr: kazumori org: np2 date: 2005-10-20 ! ! abstract: determines land surface type based on surrounding land ! surface types for AMSR-E large FOV observation ! ! program history log: ! 2005-10-20 kazumori - refered from ( subroutine deter_sfc ) ! 2006-02-01 parrish - change names of sno,isli,sst ! 2008-05-28 safford - rm unused vars ! ! input argument list: ! dlat_earth - latitude ! dlon_earth - longitude ! ! output argument list: ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! sfcpct(0:3)- percentage of 4 surface types ! (0) - sea percentage ! (1) - land percentage ! (2) - sea ice percentage ! (3) - snow percentage ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth integer(i_kind) ,intent( out) :: isflg real(r_kind),dimension(0:3),intent( out) :: sfcpct integer(i_kind) jsli,it integer(i_kind):: klat1,klon1,klatp1,klonp1 real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11 real(r_kind) :: dlat,dlon logical :: outside integer(i_kind):: klat2,klon2,klatp2,klonp2 ! ! For interpolation, we usually use o points (4points for land sea decision) ! In case of lowfreq channel (Large FOV), add the check of x points(8 points) ! (klatp2,klon1),(klatp2,klonp1) ! ---#---x---x---#--- klatp2 (klatp1,klon2),(klatp1,klonp2) ! | | | | (klat1,klon2),(klat1,klonp2) ! ---x---o---o---x--- klatp1 (klat2,klon1),(klat2,klonp1) ! | | + | | ! ---x---o---o---x--- klat1 ! | | | | ! ---#---x---x---#--- klat2 ! klon1 klonp2 ! klon2 klonp1 ! ! In total, 12 points are used to make mean sst and sfc percentage. ! it=ntguessfc if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats_sfc,nlat_sfc,1) call grdcrd1(dlon,rlons_sfc,nlon_sfc,1) end if klon1=int(dlon); klat1=int(dlat) dx =dlon-klon1; dy =dlat-klat1 dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat_sfc); klon1=min(max(0,klon1),nlon_sfc) if(klon1==0) klon1=nlon_sfc klatp1=min(nlat_sfc,klat1+1); klonp1=klon1+1 if(klonp1==nlon_sfc+1) klonp1=1 klonp2 = klonp1+1 if(klonp2==nlon_sfc+1) klonp2=1 klon2=klon1-1 if(klon2==0)klon2=nlon_sfc klat2=max(1,klat1-1) klatp2=min(nlat_sfc,klatp1+1) ! Set surface type flag. Begin by assuming obs over ice-free water sfcpct = zero jsli = isli_full(klat1 ,klon1 ) if(sno_full(klat1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klon1 ) if(sno_full(klatp1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1 ,klonp1) if(sno_full(klat1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klonp1) if(sno_full(klatp1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp2,klon1) if(sno_full(klatp2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp2,klonp1) if(sno_full(klatp2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klon2) if(sno_full(klatp1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klonp2) if(sno_full(klatp1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1,klon2) if(sno_full(klat1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1,klonp2) if(sno_full(klat1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat2,klon1) if(sno_full(klat2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat2,klonp1) if(sno_full(klat2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one sfcpct=sfcpct/12.0_r_kind ! sfcpct(3)=min(sfcpct(3),sfcpct(1)) ! sfcpct(1)=max(zero,sfcpct(1)-sfcpct(3)) isflg = 0 if(sfcpct(0) > 0.99_r_kind)then isflg = 0 else if(sfcpct(1) > 0.99_r_kind)then isflg = 1 else if(sfcpct(2) > 0.99_r_kind)then isflg = 2 else if(sfcpct(3) > 0.99_r_kind)then isflg = 3 else isflg = 4 end if return end subroutine deter_sfc_amsre_low subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc_gmi determine land surface type ! prgmmr: mkim1 org: np2 date: 2017-10-04 ! ! abstract: determines land surface type based on surrounding land ! surface types for GMI large FOV observation ! ! program history log: ! 2017-10-04 mkim1 - refered from ( subroutine deter_sfc ) ! ! input argument list: ! dlat_earth - latitude ! dlon_earth - longitude ! ! output argument list: ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! sfcpct(0:3)- percentage of 4 surface types ! (0) - sea percentage ! (1) - land percentage ! (2) - sea ice percentage ! (3) - snow percentage ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth integer(i_kind) ,intent( out) :: isflg real(r_kind),dimension(0:3),intent( out) :: sfcpct integer(i_kind) jsli,it integer(i_kind):: klat1,klon1,klatp1,klonp1 real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11 real(r_kind) :: dlat,dlon logical :: outside integer(i_kind):: klat2,klon2,klatp2,klonp2 ! ! For interpolation, we usually use o points (4points for land sea decision) ! In case of lowfreq channel (Large FOV), add the check of x points(8 points) ! (klatp2,klon1),(klatp2,klonp1) ! ---#---x---x---#--- klatp2 (klatp1,klon2),(klatp1,klonp2) ! | | | | (klat1,klon2),(klat1,klonp2) ! ---x---o---o---x--- klatp1 (klat2,klon1),(klat2,klonp1) ! | | + | | ! ---x---o---o---x--- klat1 ! | | | | ! ---#---x---x---#--- klat2 ! klon1 klonp2 ! klon2 klonp1 ! ! In total, 12 points are used to make mean sst and sfc percentage. ! it=ntguessfc if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) else dlat=dlat_earth dlon=dlon_earth call grdcrd1(dlat,rlats_sfc,nlat_sfc,1) call grdcrd1(dlon,rlons_sfc,nlon_sfc,1) end if klon1=int(dlon); klat1=int(dlat) dx =dlon-klon1; dy =dlat-klat1 dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat_sfc); klon1=min(max(0,klon1),nlon_sfc) if(klon1==0) klon1=nlon_sfc klatp1=min(nlat_sfc,klat1+1); klonp1=klon1+1 if(klonp1==nlon_sfc+1) klonp1=1 klonp2 = klonp1+1 if(klonp2==nlon_sfc+1) klonp2=1 klon2=klon1-1 if(klon2==0)klon2=nlon_sfc klat2=max(1,klat1-1) klatp2=min(nlat_sfc,klatp1+1) ! Set surface type flag. Begin by assuming obs over ice-free water sfcpct = zero jsli = isli_full(klat1 ,klon1 ) if(sno_full(klat1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klon1 ) if(sno_full(klatp1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1 ,klonp1) if(sno_full(klat1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klonp1) if(sno_full(klatp1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp2,klon1) if(sno_full(klatp2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp2,klonp1) if(sno_full(klatp2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klon2) if(sno_full(klatp1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klatp1,klonp2) if(sno_full(klatp1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1,klon2) if(sno_full(klat1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat1,klonp2) if(sno_full(klat1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat2,klon1) if(sno_full(klat2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one jsli = isli_full(klat2,klonp1) if(sno_full(klat2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 sfcpct(jsli)=sfcpct(jsli)+one sfcpct=sfcpct/12.0_r_kind ! sfcpct(3)=min(sfcpct(3),sfcpct(1)) ! sfcpct(1)=max(zero,sfcpct(1)-sfcpct(3)) isflg = 0 if(sfcpct(0) > 0.99_r_kind)then isflg = 0 else if(sfcpct(1) > 0.99_r_kind)then isflg = 1 else if(sfcpct(2) > 0.99_r_kind)then isflg = 2 else if(sfcpct(3) > 0.99_r_kind)then isflg = 3 else isflg = 4 end if return end subroutine deter_sfc_gmi subroutine deter_zsfc_model(dlat,dlon,zsfc) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_zsfc_model determine model sfc elevation ! prgmmr: parrish org: np2 date: 2006-05-23 ! ! abstract: determines model sfc elevation ! ! program history log: ! 2006-05-23 parrish ! ! input argument list: ! dlat - grid relative latitude ! dlon - grid relative longitude ! ! output argument list: ! zsfc - model surface elevation (meters) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none real(r_kind),intent(in ) :: dlat,dlon real(r_kind),intent( out) :: zsfc integer(i_kind):: klat1,klon1,klatp1,klonp1 real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11 klon1=int(dlon); klat1=int(dlat) dx =dlon-klon1; dy =dlat-klat1 dx1 =one-dx; dy1 =one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if(klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if(klonp1==nlon+1) klonp1=1 ! Interpolate zsfc to obs location zsfc=w00*zs_full(klat1,klon1 ) + w10*zs_full(klatp1,klon1 ) + & w01*zs_full(klat1,klonp1) + w11*zs_full(klatp1,klonp1) return end subroutine deter_zsfc_model subroutine reduce2full(ireduce, j, ifull) !$$$ subprogram documentation block ! . . . . ! subprogram: reduce2full find "i" index on "full" gfs grid ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: for a given "i" index on the gfs "reduced" grid, find ! the corresponding "i" index on the "full" grid. ! ! program history log: ! 2008-11-04 gayno ! ! input argument list: ! ireduce - "i" index on reduced gfs grid ! j - "j" index (same for both grids) ! ! output argument list: ! ifull - "i" index on full gfs grid ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none ! Declare passed variables. integer(i_kind), intent(in ) :: ireduce, j integer(i_kind), intent( out) :: ifull ! Declare local variables. integer(i_kind) :: ii, jj, m1, m2 real(r_kind) :: r, x1 jj = j if (j > nlat_sfc/2) jj = nlat_sfc - j + 1 m2 = lpl_gfs(jj) m1 = nlon_sfc r=real(m1)/real(m2) ii = ireduce if (ii <= 0) ii = ii + lpl_gfs(jj) if (ii > lpl_gfs(jj)) ii = ii - lpl_gfs(jj) x1=(ii-1)*r ifull=mod(nint(x1),m1)+1 return end subroutine reduce2full subroutine init_sfc(sfc_sum) !$$$ subprogram documentation block ! . . . . ! subprogram: init_sfc initialize surface structure ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: initialize to zero all elements of the surface type ! data structure. ! ! program history log: ! 2008-11-04 gayno ! ! input argument list: ! none ! ! output argument list: ! sfc_sum - holds 'sums' used in the calculation of surface fields ! within a fov. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none type surface sequence real(r_kind) :: vfr real(r_kind) :: sm real(r_kind) :: stp real(r_kind) :: ff10 real(r_kind) :: sfcr real(r_kind) :: zz real(r_kind) :: sn real(r_kind), dimension(0:3) :: ts real(r_kind), dimension(0:3) :: count real(r_kind), dimension(0:24) :: count_vty real(r_kind), dimension(0:16) :: count_sty end type surface type(surface), intent(out) :: sfc_sum ! The surface characteristics of a fov are determined by: ! ! sum(sfc_field*ant_power) / sum(ant_power) ! ! This structure holds the various 'sums' and thus must ! be initialized to zero before use. sfc_sum%vfr = zero ! greenness sfc_sum%sm = zero ! soil moisture - top layer sfc_sum%stp = zero ! soil temperature - top layer sfc_sum%ff10 = zero ! wind factor sfc_sum%sfcr = zero ! roughness length sfc_sum%zz = zero ! terrain sfc_sum%sn = zero ! snow depth sfc_sum%ts = zero ! skin temperature for each surface type ! 0-water,1-land,2-ice,3-snow sfc_sum%count = zero ! sum of the antenna power for each ! surface type sfc_sum%count_vty = zero ! vegetation type sfc_sum%count_sty = zero ! soil type return end subroutine init_sfc subroutine time_int_sfc(ix,iy,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl) !$$$ subprogram documentation block ! . . . . ! subprogram: time_int_sfc time interpolate surface data ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: time interpolate surface data to the observation time ! ! program history log: ! 2008-11-04 gayno ! ! input argument list: ! ix, iy - x/y indices of grid point ! dtsfc/dtsfcp - time interpolation factors ! itsfc/itsfcp - bounding indices of data ! ! output argument list: ! sfc_mdl - holds surface information for a single model point ! valid at the observation time ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none ! Declare passed variables. integer(i_kind), intent(in ) :: ix,iy,itsfc, itsfcp real(r_kind) , intent(in ) :: dtsfc, dtsfcp type surface2 sequence real(r_kind) :: vfr real(r_kind) :: sfcr real(r_kind) :: sm real(r_kind) :: stp real(r_kind) :: ff10 real(r_kind) :: sn real(r_kind) :: ts end type surface2 type(surface2) , intent( out) :: sfc_mdl ! Note, indices are reversed (y/x). sfc_mdl%sn=sno_full(iy,ix,itsfc)*dtsfc+sno_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%vfr=veg_frac_full(iy,ix,itsfc)*dtsfc+veg_frac_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%stp=soil_temp_full(iy,ix,itsfc)*dtsfc+soil_temp_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%sm=soil_moi_full(iy,ix,itsfc)*dtsfc+soil_moi_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%ff10=fact10_full(iy,ix,itsfc)*dtsfc+fact10_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%sfcr=sfc_rough_full(iy,ix,itsfc)*dtsfc+sfc_rough_full(iy,ix,itsfcp)*dtsfcp sfc_mdl%ts=sst_full(iy,ix,itsfc)*dtsfc+sst_full(iy,ix,itsfcp)*dtsfcp return end subroutine time_int_sfc subroutine accum_sfc(i,j,power,sfc_mdl,sfc_sum) !$$$ subprogram documentation block ! . . . . ! subprogram: accum_sfc "accumulate" surface fields ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: The surface characteristics of a fov are determined by: ! sum(sfc_field*ant_power) / sum(ant_power) ! this routine determines the required "sums". ! ! program history log: ! 2008-11-04 gayno ! ! input argument list: ! i, j - i/j indices of model grid point ! power - antenna power ! sfc_sum - holds required surface data "sums" ! ! output argument list: ! sfc_mdl - holds surface information for a single model point ! valid at the observation time ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none ! Declare passed variables. integer(i_kind), intent(in ) :: i, j real(r_kind) , intent(in ) :: power type surface2 sequence real(r_kind) :: vfr ! greenness (veg fraction) real(r_kind) :: sfcr ! roughness length real(r_kind) :: sm ! soil moisture real(r_kind) :: stp ! soil temperature real(r_kind) :: ff10 ! wind factor real(r_kind) :: sn ! snow depth real(r_kind) :: ts ! skin temperature end type surface2 type(surface2) , intent(in ) :: sfc_mdl type surface sequence real(r_kind) :: vfr ! greenness (veg fraction) real(r_kind) :: sm ! soil moisture real(r_kind) :: stp ! soil temperautre real(r_kind) :: ff10 ! wind factor real(r_kind) :: sfcr ! roughness length real(r_kind) :: zz ! terrain height real(r_kind) :: sn ! snow depth real(r_kind), dimension(0:3) :: ts ! skin temp (each land type) real(r_kind), dimension(0:3) :: count ! count of each land type real(r_kind), dimension(0:24) :: count_vty ! count of each landuse type real(r_kind), dimension(0:16) :: count_sty ! count of each soil type end type surface type(surface) , intent(inout) :: sfc_sum ! Declare local parameters. real(r_kind),parameter:: minsnow=one_tenth ! Declare local variables. integer(i_kind) :: mask, sty, vty if (power == zero) return mask=isli_full(j,i) if (mask>=1.and.sfc_mdl%sn>minsnow) mask=3 if (mask==1) then ! bare (non-snow covered) land vty=nint(veg_type_full(j,i)) sfc_sum%count_vty(vty)=sfc_sum%count_vty(vty)+power sty=nint(soil_type_full(j,i)) sfc_sum%count_sty(sty)=sfc_sum%count_sty(sty)+power sfc_sum%vfr=sfc_sum%vfr + (power*sfc_mdl%vfr) sfc_sum%sm=sfc_sum%sm + (power*sfc_mdl%sm) sfc_sum%stp=sfc_sum%stp + (power*sfc_mdl%stp) elseif (mask==3) then ! snow cover land or sea ice sfc_sum%sn=sfc_sum%sn + (power*sfc_mdl%sn) endif ! wind factor, roughness length and terrain are summed ! across all surface types. sfc_sum%ff10=sfc_sum%ff10 + (power*sfc_mdl%ff10) sfc_sum%sfcr=sfc_sum%sfcr + (power*sfc_mdl%sfcr) if (regional) then sfc_sum%zz=sfc_sum%zz + (power*zs_full(j,i)) else sfc_sum%zz=sfc_sum%zz + (power*zs_full_gfs(j,i)) endif ! keep track of skin temperature for each surface type sfc_sum%ts(mask)=sfc_sum%ts(mask) + (power*sfc_mdl%ts) ! keep count of each surface type sfc_sum%count(mask) = power + sfc_sum%count(mask) return end subroutine accum_sfc subroutine calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm, & stp,ff10,sfcr,zz,sn,ts,tsavg) !$$$ subprogram documentation block ! . . . . ! subprogram: calc_sfc calculate surface fields ! prgmmr: gayno org: np2 date: 2008-11-04 ! ! abstract: The surface characteristics of a fov are determined by: ! sum(sfc_field*ant_power) / sum(ant_power) ! ! program history log: ! 2008-11-04 gayno ! ! input argument list: ! sfc_sum - holds required surface data "sums" ! ! output argument list: ! idomsfc - dominate surface type ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! sfcpct(0:3)- percentage of 4 surface types ! (0) - sea percentage ! (1) - land percentage ! (2) - sea ice percentage ! (3) - snow percentage ! vfr - vegetation fraction ! sty - dominate soil type ! vty - dominate vegetation type ! stp - top layer soil temperature ! sm - top layer soil moisture ! ff10 - wind factor ! sfcr - surface roughness lenght ! zz - terrain height ! sn - snow depth ! ts - skin temperature for each surface type ! tsavg - average skin temperature ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none ! Declare passed variables integer(i_kind), intent( out) :: isflg, idomsfc(1) real(r_kind) , intent( out) :: sm, stp, sty, vty, vfr, sfcpct(0:3) real(r_kind) , intent( out) :: ff10, sfcr, zz, sn, ts(0:3), tsavg type surface sequence real(r_kind) :: vfr real(r_kind) :: sm real(r_kind) :: stp real(r_kind) :: ff10 real(r_kind) :: sfcr real(r_kind) :: zz real(r_kind) :: sn real(r_kind), dimension(0:3) :: ts real(r_kind), dimension(0:3) :: count real(r_kind), dimension(0:24) :: count_vty real(r_kind), dimension(0:16) :: count_sty end type surface type(surface) , intent(in ) :: sfc_sum ! Declare local variables integer(i_kind) :: itmp(1), n real(r_kind) :: count_tot, count_land, count_sty_tot, & count_snow, count_vty_tot, count_ice, count_water count_tot = sum(sfc_sum%count) ! skin temperature over entire fov tsavg = sum(sfc_sum%ts(0:3))/count_tot ! landuse is predominate type count_vty_tot = sum(sfc_sum%count_vty) if (count_vty_tot==zero) then vty=zero else itmp=lbound(sfc_sum%count_vty)-1+maxloc(sfc_sum%count_vty) vty=float(itmp(1)) endif ! soil type is predominate type count_sty_tot = sum(sfc_sum%count_sty) if (count_sty_tot==zero) then sty=zero else itmp=lbound(sfc_sum%count_sty)-1+maxloc(sfc_sum%count_sty) sty=float(itmp(1)) endif ! fields for bare (non-snow covered) land count_land = sfc_sum%count(1) if (count_land>zero) then vfr = sfc_sum%vfr / count_land stp = sfc_sum%stp / count_land sm = sfc_sum%sm / count_land ts(1) = sfc_sum%ts(1) / count_land else vfr = zero stp = zero sm = one ts(1) = tsavg endif ! fields for open water count_water=sfc_sum%count(0) if(count_water>zero)then ts(0)=sfc_sum%ts(0)/count_water else ts(0)=tsavg endif ! fields for non-snow covered sea ice count_ice=sfc_sum%count(2) if(count_ice>zero)then ts(2)=sfc_sum%ts(2)/count_ice else ts(2)=tsavg endif ! fields for snow covered land and sea ice count_snow=sfc_sum%count(3) if(count_snow>zero)then sn=sfc_sum%sn/count_snow ts(3) = sfc_sum%ts(3) / count_snow else sn=zero ts(3) = tsavg endif ! percent of each surface type sfcpct=zero do n = 0, 3 sfcpct(n) = sfc_sum%count(n) / count_tot enddo idomsfc=lbound(sfcpct)+maxloc(sfcpct)-1 ! wind factor, roughness and terrain are determined ! over entire fov. ff10 = sfc_sum%ff10/count_tot sfcr = sfc_sum%sfcr/count_tot zz = sfc_sum%zz/count_tot isflg = 0 if(sfcpct(0) > 0.99_r_kind)then isflg = 0 ! open water else if(sfcpct(1) > 0.99_r_kind)then isflg = 1 ! bare land else if(sfcpct(2) > 0.99_r_kind)then isflg = 2 ! bare sea ice else if(sfcpct(3) > 0.99_r_kind)then isflg = 3 ! snow covered land and sea ice else isflg = 4 ! mixture end if return end subroutine calc_sfc end module deter_sfc_mod