!121017: added output of elv. write(lun_output) lat, lon, elv !121017 added output of elv. ! added more stations !130204: for 5km run use domainXDim=270 and domainYDim=420 ! for 1deg rrun use domainXDim=360 and domainYDim=180 !130205: removed output of 'elv', those input parameter can ! read from the input files. ! Since the inut parameters are not output after model run, ! to keep consistency, all modifications of the input ! variables should be done in 'prep_input.m'!!! NOT in this program. !130529: add control file 'modis_params.f90' to control the 5km or 1deg ! selection. Add run for subregion function controled by ! 'domain' variable. For 5km output 'elv'. !131220: v3.2.1; Add 14 QC ground statins from Karl in get_aPixel() !140128: v3.2.2; Add 31+20 TAO and PIRATA Ocean stations !160422: added capability to deal with user-input aer model. !160501: replace the modis_params.F90 with a namelist file. PROGRAM main USE Module_SRB ,ONLY: SRB_userAer_Form1, Modflx IMPLICIT NONE integer :: i,j,k, il, ix,iy, iosts integer :: iargc integer :: argCount integer :: lun_flist = 20, & lun_input = 30, & lun_output = 40, & lun_config = 50 character :: listFile*128, kstr*32, pstr*128, configFile*128 character :: resultPath*128, & inputFile*128, & outputFile*128 integer :: year, month, day, doy, hrmin integer :: doy2mon !function declaration real :: gmt character :: yyyy*4, mm*4, dd*2, ddd*3, hhmm*4, satName*5 integer(4):: swathXDim, swathYDim logical :: swathGood real :: missInput real :: lat,lon, tLon !latitude, longitude real :: mu0 !consine of solar zenith angle real :: elv, sfc !elevation, surface type real :: alb(4) !albedo, [bsVIS,bsNIR,wsVIS,wsNIR] real :: pwtr, uo3, prs !p.water, o3, surface pressure. real :: aod !aerosol optical thickness real :: wfrc,wcod,wcre !wat. fraction, wat. optical depth, wat. drop size real :: ifrc,icod,icre !ice. fraction, ice. optical depth, ice. drop size real :: ctp,wctp,ictp !cloud top pressure. real :: tfrc !total cloud fraction. real :: albdir(7) real :: albdif(7) real :: flxall(5,3) ! all sky real :: flxclr(5,3) ! clear sky real, dimension(:,:), ALLOCATABLE ::& ssda, ssua, ssdf, vsda, vsua, vsdf,& ssdc, ssuc, stda, stua, stuc logical :: actPix logical :: isActPixel !function declarition. logical :: prnt=.FALSE. character :: header*256 integer, parameter :: NWVL=7 real, parameter :: SOLCON=1367.6!1355.04 real, parameter :: MU0MIN=0.1 real, parameter :: MU0MAX=1.0 type(SRB_UserAer_Form1) :: userAer !type(SRB_UserAer_Form2) :: userAer integer :: domainXDim, domainYDim real :: gridR !--- Define the namelist character(5) :: VERFlag = 'GRID' !='SWATH' or 'GRID' logical :: stationOnly = .FALSE. integer :: userInputAer = 0 !=0 built-in aer; =1 user aer form1, =2 user aer form2 real :: domain(4) = [ -180,180,-90,90 ] !=[w-lon, e-lon, s-lat, n-lat] real :: MISSREAL = -9999. !Missing value namelist /MODIS_SRB_Params/ VERFlag, & stationOnly, & userInputAer, & domain, & MISSREAL !-------------------- Excutable codes start here ------------------- argCount = iargc() if (argCount/=2) then print*, '--- usage: srbexe configFile inputFileList' stop endif call getarg(1, configFile) call getarg(2, listFile) !-- open the namelist file open(lun_config, file=configFile) read(lun_config, nml=MODIS_SRB_Params) if (VERFlag == 'SWATH') then domainXDim = 270 !5km swath domainYDim = 420 !5km swath else domainXDim = 360 !1-deg domain domainYDim = 180 !1-deg domain endif if ( stationOnly ) then gridR = 100 ! R=100km else gridR=huge(0.0) endif allocate( ssda(domainXDim,domainYDim), & ssua(domainXDim,domainYDim), & ssdf(domainXDim,domainYDim), & vsda(domainXDim,domainYDim), & vsua(domainXDim,domainYDim), & vsdf(domainXDim,domainYDim), & ssdc(domainXDim,domainYDim), & ssuc(domainXDim,domainYDim), & stda(domainXDim,domainYDim), & stua(domainXDim,domainYDim), & stuc(domainXDim,domainYDim) ) !-- open the file list and read in the file pathes ! open(lun_flist, file=listFile, form='formatted', & access='sequential', status='old') do il = 1,1 read(lun_flist, '(A)') kstr read(lun_flist, '(A)') pstr select case (trim(adjustl(kstr))) case ('resultPath:'); resultPath = trim(adjustl(pstr)) case DEFAULT print*, '--- main(): Invalid input path!', trim(kstr) stop end select enddo !do il = 1,1 DO !loop over files in the listfile !-- one swath a time read(lun_flist, '(A)',iostat=iosts) inputFile if (iosts /=0) EXIT !reach the end of the list file. i = index(inputFile, '.input', back=.true.) j = index(inputFile, 'modis.', back=.true.) yyyy = inputFile(i-12:i-9 ) ddd = inputFile(i-8 :i-6 ) hhmm = inputFile(i-4 :i-1 ) satName = inputFile(j+6 :i-13) if (trim(satName)=='aqua.') satName='aqua' read(yyyy, '(I4)') year read(ddd, '(I3)') doy read(hhmm, '(I4)') hrmin gmt = int(hrmin/100) + mod(hrmin,100)/60. month = doy2mon(doy,year) !-- Open input file for processing open(lun_input, file=inputFile, form='unformatted', & access='stream', status='old') read(lun_input) swathXDim, swathYDim read(lun_input) missInput !-- Open output file for writting outputFile = trim(resultPath)//& 'modis.'//trim(satName)//'.'//yyyy//ddd//'.'//hhmm//'.dat' open(lun_output, file=outputFile, form='unformatted', & access='stream', status='replace') write(lun_output) swathXDim, swathYDim swathGood=.FALSE. do iy =1,swathYDim do ix =1,swathXDim call get_aPixel( lun_input, &!swathYDim, swathXDim, & lat,lon, mu0, elv, sfc, alb, & pwtr, uo3, prs, aod, & wfrc,wcod,wcre, ifrc,icod,icre, ctp, tfrc ) !--- If station only, check if the pixel is in the neighborhood. if (stationOnly) then actPix = isActPixel(lat,lon,gridR) else tLon = lon; if (tLon>180.) tLon=tLon-360. !lon=-180..180 if (tLondomain(2) .or.& latdomain(4) ) then actPix = .FALSE. else actPix = .TRUE. endif endif !--- If within region of interest, get further information. if (actPix) then !--- if (userInputAer==1) then call getUserAerModel( userAer, lat, lon, sfc ) userAer%Tauaer = aod elseif (userInputAer==2) then else !userAer%AerID = '' !userAer%Nzsra = 0 endif !--- albdir(1:4)=alb(1) albdir(5:7)=alb(2) albdif(1:4)=alb(3) albdif(5:7)=alb(4) !--- Sanity check if ( & ( lat<-90. .or. lat>90. ) .or.& ( lon<-180. .or. lon>360. ) .or.& ( mu0<0. .or. mu0>1. ) .or.& ( any(alb<0.) .or. any(alb>1.) ) .or.& ( aod<0. ) .or.& ( tfrc<0. .or. tfrc>1. ) .or.& ( tfrc>0. .and. wcod<0. .and. icod<0. ) .or.& ( tfrc>0. .and. wcre<0. .and. icre<0. ) .or.& ( tfrc>0. .and. ctp<0. ) .or.& ( tfrc>0. .and. elv<-500. ) ) then actPix = .FALSE. endif endif !-- call SRB model ! * lat, lon, mu0, alb, aod, cld.parms, ctp, connot be missing or invalid. ! prs, pwtr, uo3 could be missing. ! elv connot be missing when cloud are precent. ! * Flxall(Ip,Iv) : Ip=1 to Nfpar, Iv=1 to Nintv, all-sky flux; ! Ip=1, top of atmosphere down-flux ! Ip=2, top of atmosphere up-flux ! Ip=3, surface down-flux ! Ip=4, surface up-flux ! Ip=5, diffuse surface down-flux ! Iv=1, shortwave ! Iv=2, visible (PAR) ! Iv=3, near infrared ! Flxclr(Ip,Iv) : as Flxall but clear-sky flux ! if ( .not.actPix ) then !CYCLE !no processing and no output. flxall(:,:) = MISSREAL flxclr(:,:) = MISSREAL else if (userInputAer==1 .or. userInputAer==2) then CALL Modflx( Prnt, Header, & NWVL, SOLCON, MU0MIN, MU0MAX, missInput,& lat, lon, month, gmt, doy, & wfrc, ifrc, pwtr, uo3, mu0, elv, prs, userAer, & wcod, icod, wcre, icre, ctp, albdir, albdif, & flxall, flxclr ) else CALL Modflx( Prnt, Header, & NWVL, SOLCON, MU0MIN, MU0MAX, missInput,& lat, lon, month, gmt, doy, & wfrc, ifrc, pwtr, uo3, mu0, elv, prs, aod, & wcod, icod, wcre, icre, ctp, albdir, albdif, & flxall, flxclr ) endif if (missInput/=MISSREAL) then where(flxall==missInput) flxall=MISSREAL where(flxclr==missInput) flxclr=MISSREAL endif swathGood=.TRUE. endif if (trim(verFlag)=='SWATH') then write(lun_output) lat, lon, elv elseif (trim(verFlag)=='GRID') then write(lun_output) lat, lon !130205 removed output of 'elv', those input parameter can read from the input files. else print*, '--- Invalid verFlag value!', verFlag STOP endif write(lun_output) flxall(:,:) write(lun_output) flxclr(:,:) ! ssda(ix,iy)=flxall(3,1) ! ssua(ix,iy)=flxall(4,1) ! ssdf(ix,iy)=flxall(5,1) ! vsda(ix,iy)=flxall(3,2) ! vsua(ix,iy)=flxall(4,2) ! vsdf(ix,iy)=flxall(5,2) ! ssdc(ix,iy)=flxclr(3,1) ! ssuc(ix,iy)=flxclr(4,1) ! stda(ix,iy)=flxall(1,1) ! stua(ix,iy)=flxall(2,1) ! stuc(ix,iy)=flxclr(2,1) !write(11,'(5F10.2)') ssda(ix,iy), ssua(ix,iy), ssdf(ix,iy), ssdc(ix,iy), stua(ix,iy) enddo !do iy =1,swathYDim enddo !do ix =1,swathXDim !-- close the current input and output file. close(lun_input) close(lun_output) if (swathGood) print*, '--- Done: '//yyyy//ddd//'.'//hhmm//'_'//satName ENDDO !loop over flist close(lun_flist) deallocate( ssda, ssua, ssdf, vsda, vsua, vsdf,& ssdc, ssuc, stda, stua, stuc ) END PROGRAM !*********************************************************************** SUBROUTINE get_aPixel( lun_input, &!numY, numX, & lat,lon, mu0, elv, sfc, alb, & pwtr, uo3, prs, aod, & wfrc,wcod,wcre, ifrc,icod,icre, ctp, tfrc ) !*********************************************************************** IMPLICIT NONE integer, intent(in) :: lun_input !integer, intent(in) :: numY,numX,iy,ix !pixel (row,col) in the swath. real(4), intent(out) :: lat,lon,mu0 !latitude, longitude, consine of solar zenith angle. real(4), intent(out) :: elv,sfc !elevation real(4), intent(out) :: alb(4) !albedo, [bsVIS,bsNIR,wsVIS,wsNIR] real(4), intent(out) :: pwtr, uo3, prs !p.water, O3, sfc.pressure. real(4), intent(out) :: aod !aerosol optical thickness real(4), intent(out) :: wfrc,wcod,wcre !wat. fraction, wat. optical depth, wat. drop size real(4), intent(out) :: ifrc,icod,icre !ice. fraction, ice. optical depth, ice. drop size real(4), intent(out) :: ctp !cloud top pressure real(4), intent(out) :: tfrc !total cloud fraction. !% sfcTyp: !% 1 evergreen needleleaf forest % 10 grasslands !% 2 evergreen broadleaf forest % 11 permanent wetlands !% 3 deciduous needleleaf forest % 12 croplands !% 4 deciduous broadleaf forest % 13 urban and built-up !% 5 mixed forests % 14 cropland/natural vegetation mosaic !% 6 closed shrubland % 15 snow and ice !% 7 open shrublands % 16 barren or sparsely vegetated !% 8 woody savannas % 0 water !% 9 savannas % -2 unclassified !% cfrc values could be NaN,0 or number(=0..1), "0" is special (means clear). !% ccre values could be NaN or number, "0" is not special. !% ccod values could be NaN or number, "0" is not special. !% cctp values could be NaN or number, "0" is not special. ! cctp in hPa ! cfrc =0..1 ! ccre in micron ! pwatr in cm ! uo3 in atm-m ! elv in meters ! prs in hPa read(lun_input) lat, lon,& ctp, wfrc,wcre,wcod, ifrc,icre,icod, tfrc,& !cld frc=0..1. aod, mu0, pwtr, uo3, elv, prs,& sfc, alb(1:4) END SUBROUTINE !*********************************************************************** FUNCTION isActPixel(lat,lon,gridR) result(actPix) !*********************************************************************** IMPLICIT NONE logical :: actPix real, intent(in) :: lat,lon,gridR integer :: i,j,k real, parameter :: PI=3.1415927 !PI=2.0*asin(1.0) real, parameter :: D2R=PI/180. real, parameter :: earthR=6378. !integer :: actStn(34) = 19+8+[(i,i=1,34)] !integer :: actStn(6) = 19+8+34+[1,2,3,4,5,6] !SURFRAD stations !integer :: actStn(19+8+34) = [(i,i=1,19+8+34)] !integer :: actStn(50) = 19+8+34+6+[(i,i=1,50)] !New 120719 stations !integer :: actStn(34+50) = [19+8+[(i,i=1,34)], 19+8+34+6+[(i,i=1,50)]] !integer :: actStn(19+8+34+6+50) = [(i,i=1,19+8+34+6+50)] !integer :: actStn(19+8+34+6+50+14) = [(i,i=1,19+8+34+6+50+14)] !All statoins. integer :: actStn(19+8+34+6+50+14+31+20) = [(i,i=1,19+8+34+6+50+14+31+20)] !All statoins+Ocean sites.. integer, parameter :: nStns = 19+8+34+6+50+14+31+20 real :: stnLatLon(2,nStns) !lat,lon integer :: it, istn real :: lonP,latP, stnLat, stnLon, latS, lonS, & dlat, dlon, dy, dx, dist !print*, stnLatLon !stop latP = lat *D2R lonP = lon *D2R IF ( lonP>PI ) lonP = lonP - 2*PI !=> -180~180 actPix = .FALSE. do it = 1, size(actStn) istn = actStn(it) stnLat = stnLatLon(1,istn) stnLon = stnLatLon(2,istn) if (stnLon>180.) stnLon=stnLon-360. !if ((stnLat<60.).AND.(stnLat>-60.)) CYCLE !Polar staions only. latS = stnLat*D2R !; %in radians lonS = stnLon*D2R !; %in radians dlat = abs(latS-latP) dlon = abs(lonS-lonP) !%-- big circle distance. dy = sqrt((cos(latP)*sin(dlon))**2 + & (cos(latS)*sin(latP) - & sin(latS)*cos(latP)*cos(dlon))**2) dx = sin(latS)*sin(latP) + & cos(latS)*cos(latP)*cos(dlon) dist = earthR*atan2(dy,dx); !dy = dlat*earthR !dx = dlon*earthR*cos(latP) !dist = sqrt(dx*dx + dy*dy) if (distMlsra .OR. nWL>Mnwl) then STOP '--- main(): number of layer or wavelength exceeds the limit!' endif userAer%AerID = 'testAerModel' userAer%Nzsra = nLay+1 userAer%Za(1:nLay+1) = [ 35., 12., 2., 0. ] userAer%Sclht(1:nLay) = [ 99., 8., 8. ] userAer%Bsra(1:nWL,1:nLay) = reshape([ 2.99E-4, 5.16E-3, 1.07E-1, & 2.54E-4, 3.90E-3, 8.19E-2, & 2.14E-4, 2.98E-3, 6.38E-2, & 1.80E-4, 2.33E-3, 5.08E-2, & 1.19E-4, 1.40E-3, 3.22E-2, & 4.31E-5, 4.80E-4, 1.36E-2, & 1.64E-5, 1.82E-4, 6.84E-3 ], & [nWL,nLay], order=[2,1] ) userAer%Osra(1:nWL,1:nLay) = reshape([ 1.000, 0.901, 0.917, & 1.000, 0.901, 0.918, & 1.000, 0.894, 0.912, & 1.000, 0.885, 0.905, & 0.999, 0.851, 0.875, & 0.997, 0.777, 0.805, & 0.408, 0.504, 0.525 ], & [nWL,nLay], order=[2,1] ) userAer%Gsra(1:nWL,1:nLay) = reshape([ 0.717, 0.699, 0.715, & 0.723, 0.686, 0.703, & 0.717, 0.674, 0.690, & 0.707, 0.661, 0.678, & 0.683, 0.641, 0.655, & 0.620, 0.633, 0.636, & 0.485, 0.758, 0.729 ], & [nWL,nLay], order=[2,1] ) END SUBROUTINE !*********************************************************************** FUNCTION doy2mon(doy,year) result (mon) !*********************************************************************** IMPLICIT NONE integer, intent(in) :: doy,year integer :: mon !1 2 3 4 5 6 7 8 9 10 11 12 integer :: nday(12)=[31,28,31,30,31,30,31,31,30,31,30,31] integer :: i,j,k,td logical :: isLeap if (isleap(year)) then if (doy<1 .or. doy>366) then print*, '---doy2mon(): Invalid value of doy!' stop endif nday(2)=29 else if (doy<1 .or. doy>365) then print*, '---doy2mon(): Invalid value of doy!' stop endif endif td=doy do i=1,12 td = td-nday(i) if (td<=0) exit enddo mon=i END FUNCTION !*********************************************************************** FUNCTION dayOfyear(year,month,day) result (jday) !*********************************************************************** IMPLICIT NONE integer, intent(in) :: year,month,day integer :: jday integer,parameter :: nday(13) = [1, 32, 60, 91, 121, 152, & 182, 213,244,274,305,335,366] integer,parameter :: nday_leap(13) = [1, 32, 61, 92, 122, 153, & 183, 214,245,275,306,336,367] logical :: isLeap if (isleap(year)) then jday = nday_leap(month)+day-1 else jday = nday(month)+day-1 endif END FUNCTION !*********************************************************************** function isLeap(yr) result(isALeap) !*********************************************************************** IMPLICIT NONE logical :: isALeap integer, intent(in) :: yr isAleap = .FALSE. if (mod(yr, 4) == 0 .AND. & mod(yr, 100) /= 0 .OR. & mod(yr, 400) == 0) then isALeap = .TRUE. endif end function