MODULE visibmod !==================================================================== ! ! Name: visibmod ! ! Function: ! Evaluate Visibility ! ! Description: ! In clear areas, aerosol Visibility is evaluated based on visible ! aerosol optical depth and planetary boundary layer depth. ! In cloudy areas, fog visibility is evaluated based on cloud optical ! depth and fog probability and fog depth. ! Cloud optical depth and effective cloud radius will be read from ! MODIS_MOD06 files. ! Aerosol optical depth will be read from MODIS_MOD04 files and then ! interpolated to pixel resolution. ! ! Reference: ! ! Inputs: ! Input passed though GEOCAT structure: satellite, nwp, out2 ! Cloud optical depth and effective cloud radius will be read from ! MODIS_MOD06 files. ! Aerosol optical depth will be read from MODIS_MOD04 files. ! ! Outputs: ! ! out2(Algo_Num)%visib = surface visibility ! ! Diagnostic outputs: ! out2(Algo_Num)%vis_type = surface visibility type ! 0 - visibility is from aerosol optical depth ! 1 - visibility is from fog visibility estimate ! ! out2(Algo_Num)%Fguess_fog_visib = first guess visibility from fog ! out2(Algo_Num)%fog_visib = bias corrected visibility from fog ! out2(Algo_Num)%Fguess_aod_visib = first guess visibility from aod ! out2(Algo_Num)%aod_visib = bias corrected visibility from aod ! out2(Algo_Num)%pblcld = 1 from cloud with planetary boundary layer ! out2(Algo_Num)%zpbl_vis - Planetary boundary layer height ! ! Public Routines: ! visib_main ! ! Public Variables: ! ! Dependencies: ! Geocat objects - visible aersol optical depth, ! nwp planetary boundary layer depth, ! fog probability, fog depth ! cloud mask ! for pblcld flag need cloud top height,surface height ! cloud optical thickness ! effective cloud drop radius ! ! Restrictions: ! No paths greater than 256 chars ! at present image not over a pole or over date line ! ! History: ! 7/14/2010 - Allen Lenzen - Version 1 - fog part of visibility to be added ! august 2010 - added fog visibility ! september 2010 -added visible aerosol optical depth from GOES East aerosol ! /smoke product (GASP) ( now removed) ! added readgaspaod subroutine (now removed) ! when aod_vis is available in out2, algorithm will use that ! instead of using readgaspaod routine (Now removed) ! April 2011 - add code to read cloud optical thickness (COT) from MODIS_MOD06 ! files ! - add code to read cloud effective radius (reff) from MODIS_MOD06 ! files ! - add code to read aerosol optical depth (AOD) from MODIS_MOD04 ! files ! May 2011 - add interpolation for aerosol optical depth from MOD04 resolution ! to image resolution ! - added updated look up tables (LUT) for aerosol visibility based ! on AOD value and ! fog visibility based on COT value ! - change fog probability to 80 percent for fog present ! - aod interpolation changed to handle Aqua data properly !==================================================================== !--- module use statements use ALGORITHM_MODULE_USAGE use ATMOSPHERE implicit none private !--- module-level subroutines and functions public:: visib_main ! main entry subroutine !--- module-level variables contains SUBROUTINE visib_main(Algo_num) !==================================================================== ! ! Name: visib ! ! Function: ! Main visibility algoritm ! ! Description: ! ! In clear areas, aerosol Visibility is evaluated based on visible ! aerosol optical depth and planetary boundary layer depth. ! In cloudy areas, fog visibility is evaluated based on cloud optical ! depth and fog probability and fog depth. ! ! Reference: ! ! Calling Sequence: CALL vis(ialgo) ! ! Inputs: ! Input passed though GEOCAT structure: satellite, nwp, out2 ! ! Outputs: ! out2(Algo_Num)%visib = surface visibility ! Diagnostic outputs: ! out2(Algo_Num)%vis_type = surface visibility type ! 0 - visibility is from aerosol optical depth ! 1 - visibility is from fog visibility estimate ! ! out2(Algo_Num)%Fguess_fog_visib = first guess visibility from fog ! out2(Algo_Num)%fog_visib = bias correcte visibility from fog ! out2(Algo_Num)%Fguess_aod_visib = first guess visibility from aod ! out2(Algo_Num)%aod_visib = bias corrected visibility from aod ! out2(Algo_Num)%pblcld = 1 from cloud with planetary boundary layer ! out2(Algo_Num)%zpbl_vis - Planetary boundary layer height ! ! ! Dependencies: ! Input passed though GEOCAT structure: satellite, nwp,out2 ! Geocat objects - visible aersol optical depth, ! nwp planetary boundary layer depth, ! fog probability, fog depth ! cloud mask ! for pblcld flag need cloud top height,surface height ! bias corrections read in with subroutine readvisbiascor ! ! Restrictions: ! ! History: ! !==================================================================== use ATMOSPHERE !--- arguments INTEGER(kind=int4), intent(in) :: Algo_num ! algorithm index !--- local variables INTEGER(kind=int4) :: Stat ! status INTEGER(kind=int4) :: iline ! line index INTEGER(kind=int4) :: ielem ! element index INTEGER(kind=int4) :: NWP_Lon_Cell INTEGER(kind=int4) :: NWP_Lat_Cell integer(kind=int4) idim(2) REAL(kind=real4), dimension(:,:), POINTER :: Visibility REAL(kind=real4), dimension(:,:), pointer :: aod_visible REAL(kind=real4), dimension(:,:), POINTER :: Fguess_aod_vis REAL(kind=real4), dimension(:,:), POINTER :: Fguess_fog_vis REAL(kind=real4), dimension(:,:), POINTER :: fog_visib REAL(kind=real4), dimension(:,:), POINTER :: aerosol_visib REAL(kind=real4), dimension(:,:), POINTER :: zpblpnt REAL(kind=real8), dimension(:,:), POINTER,save :: COT REAL(kind=real8), dimension(:,:), POINTER,save :: reff REAL(kind=real8), dimension(:,:), POINTER,save :: CWP REAL(kind=real8), dimension(:,:), POINTER,save :: AOD04 REAL(kind=real8), dimension(:,:), POINTER,save :: MOD06_Lon REAL(kind=real8), dimension(:,:), POINTER,save :: MOD06_Lat REAL(kind=real8), dimension(:,:), POINTER,save :: MOD04_Lon REAL(kind=real8), dimension(:,:), POINTER,save :: MOD04_Lat character(Len=256) :: Err_Message integer (kind=int1), dimension(:,:), POINTER :: Visib_type REAL(kind=real4) visib,visib_fog,aodvis integer(kind=int4) :: aod_cat,visib_cat integer(kind=int1), dimension(:,:), POINTER :: pblcldpnt REAL(kind=real4) :: Pblh logical notclear ! cloudy, probably cloudy, probably clear logical clear logical modis logical abi ! cloud mask 0-cloudy, 1-probably cloudy, 2-probably clear, 3-clear ! specifly number of months and catagories INTEGER(kind=int4),parameter :: nmonbias=12,ncatbias=4 ! aod and fog bias corrections and slope REAL (kind=real4),dimension(nmonbias,ncatbias),save :: aod_bias REAL (kind=real4),dimension(nmonbias,ncatbias),save :: aod_slope REAL (kind=real4),dimension(nmonbias,ncatbias),save :: fog_bias REAL (kind=real4),dimension(nmonbias,ncatbias),save :: fog_slope integer(kind=int4) modissource6,modissource4 LOGICAL first ! paths to input data character *200 dirbias integer(kind=int4) imon integer (kind=int4)nc,nr,nc6,nr6 integer (kind=int4) nc4,nr4 ! so only read input once save first,modis,abi data first/.true./ real(kind=real8), allocatable :: lwp(:,:),fog_depth(:,:) real(kind=real8) :: h2o_density real(kind=real4) :: lwc real(kind=real4) :: maxvisib data maxvisib/60./ ! data maxvisib/30./ integer(kind=int4) ifog ! algorithm for fog integer(kind=int4) :: Lines_Add ! this needs to be changed if fog is not the fourth algorithm ifog=4 allocate (lwp(sat%nx,sat%ny),fog_depth(sat%nx,sat%ny)) lwp=missing_value_real8 ! set up bias correction table on first call if(first)then ! print *,'scinfo',scinfo(sc_ind)%name first=.false. dirbias=TRIM(out2(Algo_Num)%ancil_subdir) if(scinfo(sc_ind)%series.eq.'modis')then modis=.true. abi=.false. elseif(scinfo(sc_ind)%series.eq.'goesr')then modis=.false. abi=.true. else print *,'not valid series',scinfo(sc_ind)%series stop endif write(6,*)'VISIBILITY read bias, slope files (INFORMATION) maxvisib',maxvisib call readvisbiascor (dirbias,nmonbias,ncatbias,aod_bias,aod_slope, & fog_bias,fog_slope,modis,abi) endif imon=timestr%month !------------------------------ ! Alias surface visibility Visibility => out2(Algo_Num)%visib Visib_type => out2(Algo_Num)%vis_type ! alias zpbl, pblcld, zpbl zpblpnt => out2(Algo_Num)%zpbl_vis pblcldpnt => out2(Algo_Num)%pblcld ! intialize to missing zpblpnt=missing_value_real4 Visibility=missing_value_real4 Visib_type=missing_value_int1 pblcldpnt=0 ! alias First guess visibilities Fguess_aod_vis => out2(Algo_Num)%fguess_aerosol_visib Fguess_fog_vis => out2(Algo_Num)%fguess_fog_visib ! alias visibilities fog_visib => out2(Algo_Num)%fog_visib aerosol_visib => out2(Algo_Num)%aerosol_visib aerosol_visib=missing_value_real4 Fguess_fog_vis=missing_value_real4 Fguess_aod_vis=missing_value_real4 fog_visib=missing_value_real4 ! ! Alias aerosol optical depth - visible Aod_visible => out2(Algo_num)%aod_visible aod_visible=missing_value_real4 !--- preprocessing !--- main pixel loop ! this is for Version 1 will use aod_vis when available in out2 in future nc=sat%nx nr=sat%ny ! in future for ABI if(abi)then ! AOD, ((COT,reff) or LWP ) from sat or out2 structure ! fog_probability, fog_depth from sat or out2 structure endif ! MODIS NOW if(modis)then ! read in COT from MODIS_MOD06 if(scinfo(sc_ind)%name.eq.'Terra')then modissource6=SOURCE_MODIS_MOD06 modissource4=SOURCE_MODIS_MOD04 else modissource6=SOURCE_MODIS_MYD06 modissource4=SOURCE_MODIS_MYD04 endif CALL Get_Atmosphere_Data(modissource6, & 'Cloud_Optical_Thickness', & COT, & Native_Lon=MOD06_Lon, & Native_Lat=MOD06_Lat, & Stat=Stat) IF( Stat /= STATUS_OK ) THEN WRITE (Err_Message, *) & "COT: Error from Get_Atmosphere_Data: ", & Stat CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE) RETURN END IF ! read in reff from MODIS_MOD06 CALL Get_Atmosphere_Data(modissource6, & 'Cloud_Effective_Radius', & reff, & Native_Lon=MOD06_Lon, & Native_Lat=MOD06_Lat, & Stat=Stat) IF( Stat /= STATUS_OK ) THEN WRITE (Err_Message, *) & "COT: Error from Get_Atmosphere_Data: ", & Stat CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE) RETURN END IF ! read in cloud water path from MODIS_MOD06 CALL Get_Atmosphere_Data(modissource6, & 'Cloud_Water_Path', & CWP, & Native_Lon=MOD06_Lon, & Native_Lat=MOD06_Lat, & Stat=Stat) idim=shape(MOD06_Lon) nc6=idim(1) nr6=idim(2) IF( Stat /= STATUS_OK ) THEN WRITE (Err_Message, *) & "COT: Error from Get_Atmosphere_Data: ", & Stat CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE) RETURN END IF ! read in AOD from MODIS_MOD04 !-------------------------------------------------- ! AOD at native res !-------------------------------------------------- Lines_Add=2 ! this could probably be 1 but 2 is safer and very slight effect ! on speed or memory size CALL Get_Atmosphere_Data(modissource4, & 'Optical_Depth_Land_And_Ocean', & AOD04, & Pad_Lines=Lines_Add, & Native_Lon=MOD04_Lon, & Native_Lat=MOD04_Lat, & Stat=Stat) idim=shape(MOD04_lon) nc4=idim(1) nr4=idim(2) call aod04to06(nc6,nr6,MOD06_Lon,MOD06_Lat,& nc4,nr4,MOD04_Lon,MOD04_Lat,nc,nr,aod_visible,AOD04,sat%ystart,sat%ny_all) IF( Stat /= STATUS_OK ) THEN WRITE (Err_Message, *) & "AOD small: Error from Get_Atmosphere_Data: ", & Stat CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE) RETURN END IF endif ! MODIS h2o_density=998.2071 ! kg/m3 at 1 atm/20 degrees C lwc=0.3 ! coefficient for fogdepth from lwp fog_depth=missing_value_real4 line_loop: do iline=1, sat%ny element_loop: do ielem=1, sat%nx !--- check for space pixel if (sat%space_mask(ielem,iline) == sym%SPACE) then cycle endif ! calculate lwp (liquid water path) if(COT(ielem,iline)>0 .and. reff(ielem,iline)>0)then lwp(ielem,iline)=(5./9.)*h2o_density*COT(ielem,iline)* & reff(ielem,iline)*1.e-6*1.e3 fog_depth(ielem,iline)=lwp(ielem,iline)/lwc ! units meters endif !--- define aliases for this loop iteration !--- perform operations on the current pixel NWP_Lon_Cell = sat%x_nwp(ielem,iline) NWP_Lat_Cell = sat%y_nwp(ielem,iline) Pblh = nwp%dat(NWP_Lon_Cell,NWP_Lat_Cell)%zpbl !--- write output to the out2 struct ! logical for notclear and clear notclear=sat%cldmask(ielem,iline).ne.3 clear=sat%cldmask(ielem,iline).eq.3 zpblpnt(ielem,iline)=nwp%dat(NWP_Lon_Cell,NWP_Lat_Cell)%zpbl ! set up flag for pblcld for Version 1 diagnostics if(sat%cldz(ielem,iline) > 0.0 )then if(sat%cldz(ielem,iline)-nwp%dat(NWP_Lon_Cell,NWP_lat_Cell)%zsfc & < zpblpnt(ielem,iline) ) then pblcldpnt(ielem,iline)=1 endif endif if (Aod_visible(ielem,iline) > 0.0 .and. clear) then ! use WMO value now which is 3.0 Visibility(ielem,iline)=3./(Aod_visible(ielem,iline)/(Pblh/1000.)) !bias correction of visibility ! now visib in killometers visib=Visibility(ielem,iline) Fguess_aod_vis(ielem,iline)=visib aodvis=Aod_visible(ielem,iline) aod_cat=0 if(aodvis .ge. 30.)then aod_cat=1 elseif(aodvis .ge. 10.)then aod_cat=2 elseif(aodvis .ge. 2.)then aod_cat=3 elseif(aodvis .gt. 0.)then aod_cat=4 endif if(aod_cat .ne. 0 )then if(modis)then visib=.25*visib+.75*(aod_bias(imon,aod_cat) + & aod_slope(imon,aod_cat)*Aod_visible(ielem,iline)) elseif(abi)then visib=.20*visib+.80*(aod_bias(imon,aod_cat) + & aod_slope(imon,aod_cat)*Aod_visible(ielem,iline)) endif endif ! vibibility in km ! Visibility(ielem,iline)=min(30.,visib) Visibility(ielem,iline)=min(maxvisib,visib) ! new 6/9/2011 aerosol_visib(ielem,iline)=Visibility(ielem,iline) Visib_type(ielem,iline)=0 endif zpblpnt(ielem,iline)=nwp%dat(NWP_Lon_Cell,NWP_Lat_Cell)%zpbl ! set up flag for pblcld for Version 1 diagnostics if(sat%cldz(ielem,iline) > 0.0 )then if(sat%cldz(ielem,iline)-nwp%dat(NWP_Lon_Cell,NWP_lat_Cell)%zsfc & < zpblpnt(ielem,iline) ) then pblcldpnt(ielem,iline)=1 endif endif if (.not.clear) then ! now 80 percent threshold if(out2(ifog)%MVFR_fog_prob(ielem,iline) > 80. .and. & fog_depth(ielem,iline) .ne. -999.)then if(sat%solzen(ielem,iline)<90.and.COT(ielem,iline)>0.0)then visib_fog=3.0/(COT(ielem,iline)/(fog_depth(ielem,iline)/1000.)) Fguess_fog_vis(ielem,iline)=visib_fog visib_cat=0 if(visib_fog .ge. 30.)then visib_cat=1 elseif(visib_fog .ge. 10.)then visib_cat=2 elseif(visib_fog .ge. 2.)then visib_cat=3 elseif(visib_fog .gt. 0.)then visib_cat=4 endif if(visib_cat .ne. 0 )then if(visib_cat .eq. 4)then if(modis.or.abi)then visib_fog=.4*visib_fog+.6*(fog_bias(imon,visib_cat) + & fog_slope(imon,visib_cat)*COT(ielem,iline)) endif if(Visibility(ielem,iline).ne.missing_value_real4)then write(6,*)'should not be in fog',ielem,iline endif ! Visibility(ielem,iline)=min(30.,visib_fog) Visibility(ielem,iline)=min(maxvisib,visib_fog) ! 6/9/2011 ! fog_visib(ielem,iline)=min(30.,visib_fog) fog_visib(ielem,iline)=min(maxvisib,visib_fog) ! 6/9/2011 Visib_type(ielem,iline)=1 else ! don't expect other categories write(6,*)'unexp' call flush(6) write(6,*)'unexpected fog visibility category at ',ielem, & iline,' cat ',visib_cat,' visib_fog ',visib_fog,' COT ', & COT(ielem,iline),' fog_depth ',fog_depth(ielem,iline) ! leave as a missing point for now endif endif endif endif endif end do element_loop end do line_loop !--- postprocessing IF(Sat%Iseg == (Sat%NSeg - 1)) THEN DEALLOCATE(COT) DEALLOCATE(CWP) DEALLOCATE(MOD04_Lon) DEALLOCATE(MOD04_Lat) DEALLOCATE(MOD06_Lon) DEALLOCATE(MOD06_Lat) END IF deallocate(lwp) deallocate(fog_depth) write(6,*)'bottom of visib_main INFORMATIONAL' call flush(6) !--- postprocessing end SUBROUTINE visib_main ! end main subroutine ! !==================================================================== ! ! Name: Readvisbiascor ! ! Function: ! Reads bias corrections for visibility ! ! Description: ! Reads the coefficients into a properly-shaped 2D array ! ! ! Reference: ! ! Inputs: ! directory to find bias correction coefficients ! dimension of bias ! ! Outputs: ! ! Bias correction coefficients (Bias and slope) ! ! Public Routines: ! ! Public Variables: ! ! Dependencies: ! ! ! ! Restrictions: ! No paths greater than 256 chars ! ! History: SUBROUTINE readvisbiascor(dirbias,nmonbias,ncatbias,aod_bias, & aod_slope,fog_bias,fog_slope,modis,abi) implicit none logical modis,abi integer(kind=int4) :: nmonbias,ncatbias,nmonbiasi,ncatbiasi real (kind=real4),dimension(nmonbias,ncatbias) :: aod_bias,aod_slope real (kind=real4),dimension(nmonbias,ncatbias) :: fog_bias,fog_slope character *200 filebias,dirbias*(*) ! New for AOD bias and slope ! one for ABI and one for MODIS ! for ABI if(abi)then filebias=trim(dirbias)//'/ABI.AOD.LUT.WMO.harmonic.mean.visibility.dat' endif ! for MODIS if(modis)then filebias=trim(dirbias)//'/MODIS.AOD.LUT.WMO.harmonic.mean.visibility.dat' endif open(10,file=filebias,form='unformatted') ! read in dimension of bias corrrections read(10)nmonbiasi,ncatbiasi ! check dimensions if(nmonbiasi.ne.nmonbias.or.ncatbiasi.ne.ncatbiasi)then print *,'error dimensions wrong readbiascor' stop endif ! read bias and slope read(10)aod_bias,aod_slope close(10) if(abi)then filebias=trim(dirbias)//'/ABI.COT.LUT.WMO.harmonic.mean.visibility.80prob.dat' elseif(modis)then filebias=trim(dirbias)//'/MODIS.COT.LUT.WMO.harmonic.mean.visibility.80prob.dat' endif !filebias=trim(dirbias)//'/GOES_ASOS_bias_correction.dat' open(10,file=filebias,form='unformatted') ! read in dimension of bias corrrections read(10)nmonbiasi,ncatbiasi ! check dimensions if(nmonbiasi.ne.nmonbias.or.ncatbiasi.ne.ncatbiasi)then print *,'error dimensions wrong readbiascor' stop endif ! read bias and slope read(10)fog_bias,fog_slope close(10) return end SUBROUTINE readvisbiascor ! interpolate from MOD04 resolution to MOD06 resolution( pixel resolution) ! restrictions ! not over a pole or over dateline subroutine aod04to06(nc6,nr6,lon06,lat06,nc4,nr4,lon04,lat04,nco,nro, & aod06,AOD04,ystart,nr6tot) implicit none integer(kind=int4) nco,nro,ystart integer(kind=int4)nc4,nr4,nr6tot real(kind=real8) lon04(nc4,nr4),lat04(nc4,nr4),AOD04(nc4,nr4) ! need for sat lat lon !real(kind=real4) lon06(nc6,nr6),lat06(nc6,nr6) real(kind=real8) lon06(nc6,nr6),lat06(nc6,nr6) real(kind=real4) aod06(nco,nro) real(kind=real8) lon04b(nc4,nr4),lat04b(nc4,nr4) real(kind=real8) dist2,dist2c,drad,pi,cosb real(kind=real8) lat12,lon12,lat34,lon34,aod12,aod34 real(kind=real8) delta(3),alp0,bet0 real(kind=real8) latbar integer(kind=int4) ndiff,i1,i2,i3,i4,j1,j2,j3,j4,ipass ! patthern 1 4 ! 2 3 ! but either 1 or 4 could be north most or 2 and 3 south most ! but 4 is east of 1 and 3 is east or 2 and 1 is north of 2 and 4 is north of 3 ! ! set up search pattern for search around nearest box to find box in parameter (ndiff=3) integer(kind=int4) nsearch(ndiff) data nsearch/1,8,16/ integer(kind=int4) idiff(16,ndiff),jdiff(16,ndiff) data idiff/16*0,-1,0,0,1,-1,-1,1,1,8*0,-2,0,0,0,-2,-2,2,2,-1,-1,1,1,-2,2,2,2/ data jdiff/16*0,0,-1,1,0,-1,1,1,-1,8*0,0,-2,2,2,-1,1,1,-1,-2,2,2,-2,-2,2,2,-2/ data delta/.1,.01,.001/ integer(kind=int4) maxn,i,j,jorder,ii,jj,iii,jjj integer(kind=int4)iout,nc6,nr6,iorder integer(kind=int4), allocatable :: num(:,:),ilist(:,:,:),jlist(:,:,:) integer(kind=int4) nis,nie,njs,nje integer(kind=int4) maxout,nnum,ib,jb,ibb,jbb real(kind=real8) alp,bet,latat(0:100,0:100),lonat(0:100,0:100),alpm,betm real(kind=real8)llat,llon,lonws,lones parameter (maxout=4) ! maximum boxs can search away from center box real(kind=real8) dll,a,b,lat,lon real(kind=real8) lonww(nc4-1,nr4-1),lonee(nc4-1,nr4-1),latnn(nc4-1,nr4-1) real(kind=real8) latss(nc4-1,nr4-1) real(kind=real8) lats(nc4-1,nr4-1),latn(nc4-1,nr4-1),lone(nc4-1,nr4-1) real(kind=real8) lonw(nc4-1,nr4-1) integer(kind=int4) is,ie,js,je,l,n integer jbstart integer entry save entry data entry/0/ entry=entry+1 pi=4.*atan(1.) drad=pi/180. maxn=1000 dll=0.5 if(lon04(2,1)>lon04(1,1))then iorder=1 else iorder=-1 endif if(lat04(1,2)>lat04(1,1))then jorder=1 else jorder=-1 endif ! set up offset to find pointers to points 1,2,3,4 if(iorder.eq.1)then if(jorder.eq.-1)then i1=0 j1=0 i2=0 j2=1 i3=1 j3=1 i4=1 j4=0 else i1=0 j1=1 i2=0 j2=0 i3=1 j3=0 i4=1 j4=1 endif else if(jorder.eq.-1)then i1=1 j1=0 i2=1 j2=1 i3=0 j3=1 i4=0 j4=0 else i1=1 j1=1 i2=1 j2=0 i3=0 j3=0 i4=0 j4=1 endif endif if(entry.eq.1)then print *,'VISIBILITY aod interpolation (INFORMATION) iorder', & iorder,jorder endif ! final area with non missing lat, lon on MOD04 grid is=1000000 ie=-100000 js=1000000 je=-1000000 do j=1,nr4 do i=1,nc4 if(lat04(i,j) .le. -998. .or. lon04(i,j) .le. -998. )then ! print *,'missing',i,j else is=min(i,is) js=min(j,js) ie=max(i,ie) je=max(j,je) endif end do end do ! check if a row has zeros do j=js,je do i=is+1,ie-1 if (lat04(i,j) .eq. 0.0 .and. lon04(i,j) .eq. 0.0)then if(lat04(i+1,j) .eq. 0.0 .and. lon04(i+1,j) .eq. 0.0)then print *,'cant have row of zeros' stop 'zeroes' endif endif end do end do nis=100000 njs=100000 nie=-100000 nje=-100000 latnn=1.e32 latss=1.e32 lonee=1.e32 lonww=1.e32 lats=1.e32 latn=1.e32 lone=1.e32 lonw=1.e32 do j=js,je-1 do i=is,ie-1 ! need a fix if cross +-180 longitude ! define central point of quadrilateral lon04b(i,j)=.25*(lon04(i,j)+lon04(i+1,j)+lon04(i,j+1)+lon04(i+1,j+1)) lat04b(i,j)=.25*(lat04(i,j)+lat04(i+1,j)+lat04(i,j+1)+lat04(i+1,j+1)) ii=lon04b(i,j)/dll jj=lat04b(i,j)/dll nis=min(ii,nis) njs=min(jj,njs) nje=max(jj,nje) nie=max(ii,nie) ib=i jb=j ! set up out of bounds limits for ib,jb quadrilateral area lonww(ib,jb)= & min(lon04(i+i1,j+j1),lon04(i+i2,j+j2),lon04(i+i3,j+j3),lon04(i+i4,j+j4)) lonee(ib,jb)= & max(lon04(i+i1,j+j1),lon04(i+i2,j+j2),lon04(i+i3,j+j3),lon04(i+i4,j+j4)) latnn(ib,jb)= & max(lat04(i+i1,j+j1),lat04(i+i2,j+j2),lat04(i+i3,j+j3),lat04(i+i4,j+j4)) latss(ib,jb)= & min(lat04(i+i1,j+j1),lat04(i+i2,j+j2),lat04(i+i3,j+j3),lat04(i+i4,j+j4)) ! set up bounds of inside quasi-rectangular area latn(ib,jb)=max(lat04(i+i2,j+j2),lat04(i+i3,j+j3)) lats(ib,jb)=min(lat04(i+i1,j+j1),lat04(i+i4,j+j4)) lonw(ib,jb)=min(lon04(i+i3,j+j3),lon04(i+i4,j+j4)) lone(ib,jb)=max(lon04(i+j1,j+j1),lon04(i+i2,j+j2)) end do end do allocate (num(nis:nie,njs:nje)) num=0 allocate (ilist(nis:nie,njs:nje,maxn),jlist(nis:nie,njs:nje,maxn)) do j=js,je-1 do i=is,ie-1 ! set up hash table i.e. MOD 04 lists to speed up search if(lat04(i,j).gt.-998. .and. lon04(i,j).gt.-998.)then ii=lon04b(i,j)/dll jj=lat04b(i,j)/dll if(num(ii,jj)5.and.i5.and.j+ystart-1nie.or.jjjnje)go to 1004 if(num(iii,jjj).ne.0)then do l=1,num(iii,jjj) jb=jlist(iii,jjj,l) ib=ilist(iii,jjj,l) ! determine if outside ib,jb area if(lon06(i,j)lonee(ib,jb))go to 1021 ! out if(lat06(i,j)>latnn(ib,jb))goto 1021 ! out if(lat06(i,j)=latn(ib,jb).and. & lon06(i,j)<=lonw(ib,jb).and.lon06(i,j)>=lone(ib,jb))then go to 1020 ! inside rectangle endif ! check if with big box ! check if west of west line i1,i2 a=lon04(ib+i2,jb+j2) b=(lon04(ib+i1,jb+j1)-a)/(lat04(ib+i1,jb+j1)-lat04(ib+i2,jb+j2)) lon=a+b*(lat06(i,j)-lat04(ib+i2,jb+j2)) lonws=lon if(lon06(i,j)lon)then go to 1021 ! outside endif ! check if south of north line i1,i4 a=lat04(ib+i1,jb+j1) b=(lat04(ib+i4,jb+j4)-a)/(lon04(ib+i4,jb+j4)-lon04(ib+i1,jb+j1)) lat=a+b*(lon06(i,j)-lon04(ib+i1,jb+j1)) if(lat06(i,j)>lat)then go to 1021 ! outside endif go to 1020 ! inside big pixel 1021 continue end do ! l endif 1004 continue end do ! n end do ! iout go to 999 1020 continue ! box is ib,jb ! since in box refine search to get alp,bet for 06 point dist2=1.e10 alp0=0.0 bet0=0.0 do ipass=1,2 do jj=0,10 do ii=0,10 alp=delta(ipass)*float(ii)+alp0 bet=delta(ipass)*float(jj)+bet0 lat12=lat04(ib+i1,jb+j1)*bet+lat04(ib+i2,jb+j2)*(1.-bet) lat34=lat04(ib+i4,jb+j4)*bet+lat04(ib+i3,jb+j3)*(1.-bet) lon12=lon04(ib+i1,jb+j1)*bet+lon04(ib+i2,jb+j2)*(1.-bet) lon34=lon04(ib+i4,jb+j4)*bet+lon04(ib+i3,jb+j3)*(1.-bet) latat(ii,jj)=alp*lat34+(1.-alp)*lat12 lonat(ii,jj)=alp*lon34+(1.-alp)*lon12 latbar=.5*(lat06(i,j)+latat(ii,jj))*drad cosb=cos(latbar) dist2c=(lat06(i,j)-latat(ii,jj))**2.+ & ((lon06(i,j)-lonat(ii,jj))*cosb)**2. if(dist2c=0.0.and.AOD04(ib+i2,jb+j2)>=0.0.and. & AOD04(ib+i3,jb+j3)>=0.0.and.AOD04(ib+i4,jb+j4)>=0.0)then aod12=AOD04(ib+i1,jb+j1)*bet+AOD04(ib+i2,jb+j2)*(1.-bet) aod34=AOD04(ib+i4,jb+j4)*bet+AOD04(ib+i3,jb+j3)*(1.-bet) aod06(i,j)=alp*aod34+(1.-alp)*aod12 endif 999 continue endif end do ! i end do ! j deallocate (num,ilist,jlist) return end subroutine aod04to06 end MODULE visibmod