!------------------------------------------------------------------ ! this is module of surface properties (reflectance, emissivity) for ! the UMD surface type classification ! ! ! 0: Water ! 1: Evergreen Needleleaf Forests ! 2: Evergreen Broadleaf Forests ! 3: Deciduous Needleleaf Forests ! 4: Deciduous Broadleaf Forests ! 5: Mixed Forests ! 6: Woodlands ! 7: Wooded Grasslands/Shrubs ! 8: Closed Bushlands or Shrublands ! 9: Open Shrublands ! 10: Grasses ! 11: Cropland ! 12: Bare ! 13: Urban and Built ! 14: Snow (not in original, based on input snow data) ! 15: Sea_Ice (not in original, based on input snow data) ! ! Note, surface albedoes range from 0 to 1.0 ! ! Author: Andrew Heidinger, NOAA ! ! 2014-10-22 - Modified by Vladimir Kondratovich, IMSG: Introduced ingest ! of full Snow IMS dataset for NH pixels ! ! Private Routines: None ! !------------------------------------------------------------------ module SURFACE_PROPERTIES use CONSTANTS use FUNDAMENTAL_CONSTANTS use PIXEL_COMMON use NWP_COMMON use SNOW_4KM implicit none private public:: SETUP_UMD_PROPS, COMPUTE_SNOW_FIELD, GET_PIXEL_SFC_TYPE, & GET_PIXEL_SFC_EMISS_FROM_SFC_TYPE, GET_SFC_PROPS, GET_SFC_PROPS_GRID integer, parameter, public:: ntype_sfc=16,nspec_sfc=6, ntype_sfc_old=15 integer, parameter, public:: nzen_emiss=10 real(kind=real4), public, save:: delta_zen_emiss real(kind=real4), dimension(nzen_emiss), public, save:: zen_emiss real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch1_sfc_alb_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch2_sfc_alb_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch3a_sfc_alb_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch3b_sfc_alb_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch3b_sfc_emiss_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch4_sfc_emiss_umd real(kind=real4), dimension(0:ntype_sfc-1), public,save:: ch5_sfc_emiss_umd real(kind=real4), dimension(0:ntype_sfc_old-1,nspec_sfc), public,save:: spec_albedo real(kind=real4), dimension(0:ntype_sfc_old-1), public, save:: alb_solzen_factor !------ integer, parameter, public:: num_lat_land_mask = 2880, num_lon_land_mask = 5760 real, parameter, public:: first_lon_land_mask = -180.0, first_lat_land_mask = 90.0, & last_lon_land_mask = 180.0, last_lat_land_mask = -90.0 real, parameter, public:: del_lon_land_mask = (last_lon_land_mask - first_lon_land_mask)/num_lon_land_mask real, parameter, public:: del_lat_land_mask = (last_lat_land_mask - first_lat_land_mask)/num_lat_land_mask integer(kind=int1),dimension(num_lon_land_mask,num_lat_land_mask), public, save:: & land_cover_umd contains !----------------------------------------------------------------- ! read into memory the radiative properties for each UMD type !----------------------------------------------------------------- subroutine SETUP_UMD_PROPS() !--- channel 1 surface albedo based on PATMOS-x ! ch1_sfc_alb_umd = (/0.0300, & !0 ch1_sfc_alb_umd = (/0.0500, & !0 0.0406, & !1 0.0444, & !2 0.0406, & !3 0.0542, & !4 0.0423, & !5 0.0423, & !6 0.0576, & !7 0.1242, & !8 0.1781, & !9 0.0788, & !10 0.0677, & !11 0.3124, & !12 0.0918, & !13 0.9763, & !14 0.9763/) !15 !--- channel 2 surface albedo based on PATMOS-x ch2_sfc_alb_umd = (/0.0163, & !0 0.1991, & !1 0.2179, & !2 0.1991, & !3 0.2943, & !4 0.2500, & !5 0.2500, & !6 0.2700, & !7 0.3072, & !8 0.2942, & !9 0.2919, & !10 0.2520, & !11 0.4019, & !12 0.2466, & !13 0.8194, & !14 0.8194/) !15 !--- channel 3a surface albedo based on PATMOS-x ch3a_sfc_alb_umd = (/0.0163, & !0 0.2020, & !1 0.2210, & !2 0.2020, & !3 0.2980, & !4 0.2500, & !5 0.2500, & !6 0.2700, & !7 0.3400, & !8 0.2580, & !9 0.2980, & !10 0.2520, & !11 0.3900, & !12 0.2650, & !13 0.1743, & !14 0.1743/) !15 !--- channel 3b surface emissivity based on SeeBor et al 2006 ch3b_sfc_emiss_umd =(/0.985, & !0 0.95, & !1 0.96, & !2 0.94, & !3 0.94, & !4 0.94, & !5 0.94, & !6 0.93, & !7 0.89, & !8 0.89, & !9 0.92, & !10 0.94, & !11 0.86, & !12 0.94, & !13 0.984, & !14 0.984/) !15 !--- channel 3b surface albedo based on SeeBor et al 2006 ch3b_sfc_alb_umd = 1.0 - ch3b_sfc_emiss_umd !--- channel 4 surface emissivity based on SeeBor et al 2006 ch4_sfc_emiss_umd = (/0.985, & !0 0.965, & !1 0.961, & !2 0.963, & !3 0.963, & !4 0.966, & !5 0.962, & !6 0.961, & !7 0.961, & !8 0.958, & !9 0.969, & !10 0.965, & !11 0.955, & !12 0.960, & !13 0.979, & !14 0.979/) !15 !--- channel 5 surface emissivity based on SeeBor et al 2006 ch5_sfc_emiss_umd = (/0.985, & !0 0.967, & !1 0.965, & !2 0.965, & !3 0.966, & !4 0.968, & !5 0.965, & !6 0.966, & !7 0.969, & !8 0.968, & !9 0.964, & !10 0.969, & !11 0.971, & !12 0.962, & !13 0.977, & !14 0.977/) !15 !--------- ADDED in for old COD algorithm - WCS3 alb_solzen_factor = reshape((/& 0.41,& ! (17) OCEAN WATER 0.40,& ! ( 1) EVERGREEN NEEDLE FOR 0.44,& ! ( 2) EVERGREEN BROAD FOR 0.32,& ! ( 3) DECIDUOUS NEEDLE FOR 0.39,& ! ( 4) DECIDUOUS BROAD FOR 0.22,& ! ( 5) MIXED FOREST 0.22,& ! ( 5) MIXED FOREST 0.15,& ! ( 8) WOODY SAVANNA 0.28,& ! ( 6) CLOSED SHRUBS 0.40, & ! ( 7) OPEN/SHRUBS 0.22,& ! (10) GRASSLAND 0.12,& ! (14) CROP MOSAIC 0.40,& ! (16) BARREN/DESERT 0.58,& ! (18) TUNDRA 0.10/),& ! (15) Antartic Snow (/ntype_sfc_old/)) spec_albedo = reshape((/& 0.0163,0.0163,0.0163,0.0163,0.0163,0.0163,& ! (17) OCEAN WATER 0.0406,0.1991,0.2020,0.2020,0.2020,0.2020,& !( 1) EVERGREEN NEEDLE FOR 0.0444,0.2179,0.2210,0.2210,0.2210,0.2210,& !( 2) EVERGREEN BROAD FOR 0.0406,0.1991,0.2020,0.2020,0.2020,0.2020,& ! ( 3) DECIDUOUS NEEDLE FOR 0.0542,0.2943,0.2980,0.2980,0.2980,0.2980,& ! ( 4) DECIDUOUS BROAD FOR 0.0423,0.2500,0.2500,0.2500,0.2500,0.2500,& ! ( 5) MIXED FOREST 0.0423,0.2500,0.2500,0.2500,0.2500,0.2500,& ! ( 5) MIXED FOREST 0.0576,0.2700,0.2700,0.2700,0.2700,0.2700,& ! ( 8) WOODY SAVANNA 0.1242,0.3072,0.3400,0.3400,0.3400,0.3400,& ! ( 6) CLOSED SHRUBS 0.1781,0.2942,0.2580,0.2580,0.2580,0.2580,& ! ( 7) OPEN/SHRUBS 0.0788,0.2919,0.2980,0.2980,0.2980,0.2980,& ! (10) GRASSLAND 0.0677,0.2520,0.2520,0.2520,0.2520,0.2520,& ! (14) CROP MOSAIC 0.3124,0.4019,0.3900,0.3900,0.3900,0.3900,& ! (16) BARREN/DESERT 0.0918,0.2466,0.2650,0.2650,0.2650,0.2650,& ! (18) TUNDRA 0.9763,0.8194,0.1743,0.1310,0.1310,0.1310/),& ! (15) ANTARCTIC SNOW (/ntype_sfc_old,nspec_sfc/),order=(/2,1/)) !------ zenith angles for emissivity tables zen_emiss = reshape((/ & 0.000,10.00,20.00,30.00,40.00,50.00,60.00,70.00,80.00,90.0/),(/nzen_emiss/)) delta_zen_emiss = zen_emiss(2) - zen_emiss(1) end subroutine SETUP_UMD_PROPS !--------------------------------------------------------------------------- ! routine to compute surface radiative properties for AVHRR channels ! input: ! lat, lon - latitude and longitude in degrees ! solzen - solar zenith angle in degrees ! output: alb1,alb2,alb3 - surface albedos as a fraction of 1 ! ems3,ems4,ems5 - surface emissivities as a fraction of 1 ! ! note - ems3 is computed as 1.0 - alb3 even for nighttime. ! ! NOTE 2 - This is only used for OLR calculation. Also, variables are all ! internal to subroutine. So current GOES numbering was used because ! of consistancy for older GSIP releases, thus ABI numbering wasn't ! used. ! !---------------------------------------------------------------------------- subroutine GET_SFC_PROPS(scene,solzen,satzen,alb1,alb2,alb3a,alb3b, & ems3b,ems4,ems5,ems6) integer(kind=int1),intent(in):: scene real (kind=ipre), intent(in):: solzen,satzen real (kind=ipre), intent(out):: alb1,alb2,alb3a,alb3b,ems3b,ems4,ems5,ems6 real (kind=real4):: zenadj,xzen, conszener integer:: izen if (solzen < 90.0) then zenadj = ( 1.0+alb_solzen_factor(scene) ) / & ( 1.0+2.0*alb_solzen_factor(scene)* & cos(solzen*3.14159265/180.0) ) else zenadj = 1.0 endif conszener=cos(satzen*dtor) !Base off of ABI EVENTUALLY FOR NOW USE GOES CHANNELS IN HERE alb1 = zenadj * ch1_sfc_alb_umd(scene) !.64 micron (GOES 1) alb2 = zenadj * ch2_sfc_alb_umd(scene) !.81 micron alb3a = zenadj * ch3a_sfc_alb_umd(scene) ! 1.6 micron alb3b = zenadj * ch3b_sfc_alb_umd(scene) ! 3.9 micron(GOES2) ems3b = zenadj * ch3b_sfc_emiss_umd(scene) izen = max(1,min(nzen_emiss-1, & int((satzen - zen_emiss(1))/delta_zen_emiss) + 1)) xzen = (satzen - zen_emiss(izen))/(zen_emiss(izen+1) - zen_emiss(izen)) ems4 = ch4_sfc_emiss_umd(scene) ! 11 micron (GOES4) ems5 = ch5_sfc_emiss_umd(scene) ! 12 micron (GOES5) if (scene == 0) then ems4 = 0.844780 + 0.328921 * conszener -0.182375*(conszener**2) ems5 = 0.775019 + 0.474005 * conszener -0.261739*(conszener**2) endif !----------------- no idea what to do for ch6 ems6 = ems5 ! 13 micron (GOES6) end subroutine GET_SFC_PROPS !----------------------------------------------------------------------- ! If no sfc emiss data base read in, base it on the surface type. ! also derive surface reflectance from surface emissivity ! also reset sfc emiss if snow covered !----------------------------------------------------------------------- subroutine GET_PIXEL_SFC_EMISS_FROM_SFC_TYPE(j1,j2) integer, intent(in):: j1, j2 integer:: i, j j_loop: do j = j1,j1+j2-1 i_loop: do i = 1,num_pix !--- check for a bad pixel if (bad_pix_mask(14,i,j) == sym%YES) then cycle endif if (space_mask(i,j) == sym%YES) then cycle endif !--- based on surface type, assign surface emissivities if seebor not used if (use_seebor_opt == sym%NO) then if (sat_info_gsip(1)%chanon(7) > 0) sfc_emiss_7(i,j) = ch3b_sfc_emiss_umd(sfc_type(i,j)) if (sat_info_gsip(1)%chanon(14) > 0) sfc_emiss_14(i,j) = ch4_sfc_emiss_umd(sfc_type(i,j)) if (sat_info_gsip(1)%chanon(15) > 0) sfc_emiss_15(i,j) = ch5_sfc_emiss_umd(sfc_type(i,j)) if (sat_info_gsip(1)%chanon(16) > 0) sfc_emiss_16(i,j) = ch5_sfc_emiss_umd(sfc_type(i,j)) !If SEEBOR not used endif !--- for ocean, use this parameterization from Nick Nalli - overwrite seebor if (sfc_type(i,j) == 0) then if (sat_info_gsip(1)%chanon(7) > 0) sfc_emiss_7(i,j) = emiss_water_sfc(7) if (sat_info_gsip(1)%chanon(14) > 0) sfc_emiss_14(i,j) = emiss_water_sfc(14) if (sat_info_gsip(1)%chanon(15) > 0) sfc_emiss_15(i,j) = emiss_water_sfc(15) if (sat_info_gsip(1)%chanon(16) > 0) sfc_emiss_16(i,j) = emiss_water_sfc(16) endif !--- if (snow_mask if (snow_mask(i,j) /= sym%NO_SNOW) then if ((sfc_type(i,j) /= sym%EVERGREEN_NEEDLE_SFC) .and. & (sfc_type(i,j) /= sym%EVERGREEN_BROAD_SFC)) then if (sat_info_gsip(1)%chanon(7) > 0) sfc_emiss_7(i,j) = ch3b_sfc_emiss_umd(14) if (sat_info_gsip(1)%chanon(14) > 0) sfc_emiss_14(i,j) = ch4_sfc_emiss_umd(14) if (sat_info_gsip(1)%chanon(15) > 0) sfc_emiss_15(i,j) = ch5_sfc_emiss_umd(14) if (sat_info_gsip(1)%chanon(16) > 0) sfc_emiss_16(i,j) = ch5_sfc_emiss_umd(14) endif endif end do i_loop end do j_loop end subroutine GET_PIXEL_SFC_EMISS_FROM_SFC_TYPE !--------------------------------------------------------------------------- ! routine to compute surface radiative properties for AVHRR channels. Used in ! the calculation of longwave fluxes at grid level ! ! input: !---------------------------------------------------------------------------- subroutine GET_SFC_PROPS_GRID(scene,solzen,satzen,alb1,alb3b, & ems3b,ems4) integer(kind=int1),intent(in):: scene real (kind=ipre), intent(in):: solzen,satzen real (kind=ipre), intent(out):: alb1,alb3b,ems3b,ems4 real (kind=real4):: zenadj,xzen, conszener integer:: izen if (solzen < 90.0) then zenadj = ( 1.0+alb_solzen_factor(scene) ) / & ( 1.0+2.0*alb_solzen_factor(scene)* & cos(solzen*3.14159265/180.0) ) else zenadj = 1.0 endif conszener=cos(satzen*dtor) !Base off of ABI EVENTUALLY FOR NOW USE GOES CHANNELS IN HERE alb1 = zenadj * ch1_sfc_alb_umd(scene) !.64 micron alb3b = zenadj * ch3b_sfc_alb_umd(scene) ! 3.9 micron ems3b = zenadj * ch3b_sfc_emiss_umd(scene) izen = max(1,min(nzen_emiss-1, & int((satzen - zen_emiss(1))/delta_zen_emiss) + 1)) xzen = (satzen - zen_emiss(izen))/(zen_emiss(izen+1) - zen_emiss(izen)) ems4 = ch4_sfc_emiss_umd(scene) ! 11 micron (GOES4) if (scene == 0) then ems4 = 0.844780 + 0.328921 * conszener -0.182375*(conszener**2) endif end subroutine GET_SFC_PROPS_GRID !------------------------------------------------------------------------------- !--- populate the snow_input array based on all available sources of snow data !------------------------------------------------------------------------------- subroutine COMPUTE_SNOW_FIELD(j1,j2) integer, intent(in):: j1,j2 integer(kind=int4):: i, j, inwp,jnwp, snow_input, xp, xm, yp, ym real(kind=real4) :: x, y, snowsum, weight, latv, lonv, dist integer(kind=int4), parameter :: nneib=4 real(kind=real4), dimension(nneib) :: snowval real, parameter:: missing_gfs = 9.999E+20 snow_mask = 0 do j = j1,j1+j2-1 do i = 1, num_pix snow_input = 0 !--- check for space scans if (space_mask(i,j) == sym%YES) then cycle end if inwp = i_nwp(i,j) jnwp = j_nwp(i,j) if (sfc_type(i,j) == 14) then snow_input = 1 end if if (use_sst_anal .eq. sym%YES) then if ((sfc_type(i,j) == 0) .and. (sst_anal_cice(i,j) > 10) .and. (sst_anal_cice(i,j) /= 122)) then snow_input = 1 end if end if if (nwp_flag > 0 .and. inwp > 0 .and. jnwp > 0) then if (weasd_nwp(inwp, jnwp) > 1.0 .and. weasd_nwp(inwp, jnwp) /= missing_gfs) then snow_input = 1 endif end if !------------------------------------------------------- ! use some basic tests to modify snow_input ! remember results for snow_frac !------------------------------------------------------- snow_mask(i,j) = snow_input !-- can't be snow if warm if ((scinfo(sc_ind)%chanon(14) == sym%YES) .and. (bt14(i,j) > 277.0)) then snow_mask(i,j) = 0 end if !--- some day-specific tests if (solzen(i,j) < 85.0) then ! day check !--- conditions where snow is not possible if ((scinfo(sc_ind)%chanon(2) == sym%YES) .and. (ref2(i,j) < 10.0)) then snow_mask(i,j) = 0 end if if ((scinfo(sc_ind)%chanon(7) == sym%YES).and.(ref7(i,j) > 5.0)) then snow_mask(i,j) = 0 end if !---- conditions where I will add snow above that in ancillary data !---- note, atmospherically corrected values not available to this routine ! --- WCS NOTE: ONLY VALID IF CHANNELS ARE READ IN. Also no clear composite test. if (scinfo(sc_ind)%chanon(14) == sym%YES) then if (scinfo(sc_ind)%chanon(2) == sym%YES) then if (scinfo(sc_ind)%chanon(7) == sym%YES) then if (sfc_type(i,j) /= 0 .and. inwp > 0 .and. jnwp > 0) then if (abs(bt14(i,j) - tmpsfc_nwp(inwp,jnwp)) < 10.0) then if (bt14(i,j) < 277.0 ) then if (ref2(i,j) > 20.0 ) then if ((ref7(i,j) < 3.0)) then snow_mask(i,j) = 1 end if end if end if end if end if end if end if end if end if !day check !---------------------------------------------------------------------- !--- set snow pixels over desert to be not desert !---------------------------------------------------------------------- if (snow_mask(i,j) == sym%YES) then desert_mask(i,j) = sym%NO end if !--- evaluation of the snow_mask ended. Now start evaluation of the snow_frac snow_frac(i,j)=real(snow_mask(i,j)) !--- this is it in southern hemisphere (no Snow IMS data there). !--- Use Snow IMS information in the Northern Hemisphere if ((lat(i,j) >= 0.).and.use_ims_command) then snow_frac(i,j)=0. call ll2xy_ims4(lat(i,j),lon(i,j),x,y) !--- Establish square compising the pixel center xm=floor(x) xp=xm+1 ym=floor(y) yp=ym+1 !--- Discard squares sticking out of the map if ((xm.lt.1).or.(xp.gt.x_ims)) cycle if ((ym.lt.1).or.(yp.gt.y_ims)) cycle !--- Analyze snow in verices according to Snow IMS data snowval=0. if (snow_ims(yp,xm).gt.2) snowval(1)=1. if (snow_ims(ym,xm).gt.2) snowval(2)=1. if (snow_ims(yp,xp).gt.2) snowval(3)=1. if (snow_ims(ym,xp).gt.2) snowval(4)=1. snowsum=sum(snowval) !--- If everywhere is snow or no snow, assign the same to the pixel: if (snowsum.eq.0.) then snow_frac(i,j)=0. else if (snowsum.eq.4.) then snow_frac(i,j)=1. else !--- This is boundary case - interpolate. Go from vertex to vertex: snowsum=0. weight=0 call xy2ll_ims4(xm,yp,latv,lonv) dist=sqrt((latv-lat(i,j))**2+(lonv-lon(i,j))**2) if (dist.eq.0.) then snow_frac(i,j)=snowval(1) go to 1 else snowsum=snowsum+snowval(1)/dist weight=weight+1./dist end if call xy2ll_ims4(xm,ym,latv,lonv) dist=sqrt((latv-lat(i,j))**2+(lonv-lon(i,j))**2) if (dist.eq.0.) then snow_frac(i,j)=snowval(2) go to 1 else snowsum=snowsum+snowval(2)/dist weight=weight+1./dist end if call xy2ll_ims4(xp,yp,latv,lonv) dist=sqrt((latv-lat(i,j))**2+(lonv-lon(i,j))**2) if (dist.eq.0.) then snow_frac(i,j)=snowval(3) go to 1 else snowsum=snowsum+snowval(3)/dist weight=weight+1./dist end if call xy2ll_ims4(xp,ym,latv,lonv) dist=sqrt((latv-lat(i,j))**2+(lonv-lon(i,j))**2) if (dist.eq.0.) then snow_frac(i,j)=snowval(4) go to 1 else snowsum=snowsum+snowval(4)/dist weight=weight+1./dist end if !--- At this end, we get truly snow fraction (for the boundary) snow_frac(i,j)=snowsum/weight end if 1 continue !--- Here we end analysis of the Snow IMS data (yesterday snow). end if pix_snow_frac(i,j)=snow_frac(i,j) snow_mask(i,j)=min(1,max(0,nint(snow_frac(i,j)))) end do ! i_loop end do ! j_loop end subroutine COMPUTE_SNOW_FIELD !----------------------------------------------------------------------- ! this routine interpolate the UMD surface type to each AVHRR pixel ! using nearest neighor sampling !----------------------------------------------------------------------- subroutine GET_PIXEL_SFC_TYPE(j1,j2) integer, intent(in):: j1, j2 integer:: ilon_land_mask, ilat_land_mask integer:: i, j i_loop: do j = j1,j2 j_loop: do i = 1,num_pix !------ compute pixel surface type ilon_land_mask = max(1,min(num_lon_land_mask, & int((lon(i,j) - first_lon_land_mask)/del_lon_land_mask) + 1)) ilat_land_mask = max(1,min(num_lat_land_mask, & int((lat(i,j) - first_lat_land_mask)/del_lat_land_mask) + 1)) !--- compute desert flags from sfc type desert_mask(i,j)=sym%NO_DESERT if ((sfc_type(i,j)==12)) desert_mask(i,j)=sym%BRIGHT_DESERT if ((sfc_type(i,j) <= 12).and.(sfc_type(i,j) >= 7)) desert_mask(i,j)=sym%NIR_DESERT end do j_loop end do i_loop end subroutine GET_PIXEL_SFC_TYPE end module SURFACE_PROPERTIES