PROGRAM MODELSELECT C**************************************************************************** C Program to select model data at user input Class I Area sites, and write C the results to an ASCII file, along with other site-specific inputs C read in on csv file needed to perform regional haze impact estimates. Also C write model results at the selected sites to netCDF for visualization/QA C*************************************************************************** USE M3UTILIO IMPLICIT NONE c17 ifort -fPIC -check -traceback -extend_source -zero -O3 -mp1 -o model_select.exe model_select.F -L/home/local-rhel7/apps/netcdf-4.4.1/intel-17.0/lib -lnetcdf -lnetcdff -L/home/local-rhel7/apps/ioapi-3.2/intel-17.0/lib -lioapi -module /home/local-rhel7/apps/ioapi-3.2/intel-17.0/Linux2_x86_64ifort C........... Parameters CHARACTER *16 PNAME PARAMETER (PNAME='MODELSELECT') INTEGER, PARAMETER :: MXTCOL = 40 ! max columns in sitefile CHARACTER*1, PARAMETER:: SEP = '/' ! separator for calendar date C........... Local variables CHARACTER *16 INFILE ! input netCDF model data file CHARACTER *16 SITEFILE ! input (csv) monitor site file CHARACTER *16 OUTFILE ! output netCDF model data file CHARACTER *16 ASCIIFILE ! output ASCII model data file CHARACTER *16 DEFAULT CHARACTER *256 MESG CHARACTER(LEN=1000) LINE ! input line buffer INTEGER NCOLS ! number of grid columns INTEGER NROWS ! number of grid rows INTEGER NLAYS ! number of vertical layers INTEGER NVARS ! number of species INTEGER NSPC ! number of chemical species for ASCII output INTEGER NSITES ! number of Class I areas INTEGER NSRC ! number of source tags for ASCII output INTEGER NMON ! number of months of RH factors (ASCII output) INTEGER NVMAX ! ASCII output species maximum INTEGER TSTEP ! time step of output file INTEGER NSECS ! time step (s) of output file INTEGER JDATE INTEGER JTIME INTEGER YEAR, IYEAR! date vars INTEGER IDAY ! day var INTEGER IMON ! month var CHARACTER*20 DATE ! calendar date INTEGER NSTEP INTEGER STDATE, STTIME ! start date and time (daily avg file) INTEGER EDATE, ETIME ! end date and time (daily avg file) INTEGER STATUS ! allocation-status CHARACTER*30 SEGMENT( MXTCOL ) ! Input line fields REAL LAT, LON ! site lat, signed lon, unsigned lon for LL2LAM REAL DX, DY, XORIG, YORIG ! x and y grid spacing (m), origin x and y (m) REAL ANG1, ANG2, ANG3 ! projection angles REAL XCEN, YCEN ! projection center REAL X, Y ! site x and y (m) (output of LL2LAM) INTEGER C ! species loop counter INTEGER R ! species loop counter INTEGER L INTEGER S INTEGER V INTEGER T ! time step loop counter INTEGER ISITE ! monitor site counter LOGICAL FIRSTIME ! firstime flag for site loop CHARACTER*16 STR ! scratch string variable INTEGER GDTYP INTEGER SITEUNIT INTEGER ASCIIUNIT CHARACTER*32 FMTHDR CHARACTER*64 FMTDATA LOGICAL, ALLOCATABLE:: MASK(:,:) REAL, ALLOCATABLE:: A(:,:,:,:) REAL, ALLOCATABLE:: RAY(:,:) REAL, ALLOCATABLE:: FSRH(:,:,:) REAL, ALLOCATABLE:: FLRH(:,:,:) REAL, ALLOCATABLE::FSSRH(:,:,:) REAL, ALLOCATABLE:: FRH(:,:,:) REAL, ALLOCATABLE:: AIN(:,:,:,:) REAL, ALLOCATABLE:: OUT(:,:,:,:) CHARACTER *16, ALLOCATABLE:: SCNAME(:) ! source tag array CHARACTER *5, ALLOCATABLE:: SPNAME(:) ! chemical species array CHARACTER*5, ALLOCATABLE::SITENAME(:) INTEGER, ALLOCATABLE:: ROWST(:) ! row # of site INTEGER, ALLOCATABLE:: COLST(:) ! col # of site C ======================================================================= C...... beginning of the code DEFAULT = 'INFILE' CALL GETSTR('Enter name for input NetCDF data file ', $ DEFAULT, INFILE) DEFAULT = 'OUTFILE' CALL GETSTR('Enter name for output averaged NetCDF data file ', $ DEFAULT, OUTFILE) DEFAULT = 'SITEFILE' CALL GETSTR('Enter name for monitor site data file ', $ DEFAULT, SITEFILE) SITEUNIT = GETEFILE(SITEFILE, .FALSE., .TRUE., PNAME) IF (SITEUNIT .LT. 0) then CALL M3EXIT (PNAME, 0, 0, & 'SITEFILE file opening error', 2 ) END IF ! if site input failed STDATE= 0 STDATE = GETNUM( 0, 9999999,STDATE,'Enter start date') STTIME = 0 STTIME = GETNUM( 0, 9999999,STTIME,'Enter start time') EDATE = 0 EDATE = GETNUM( 0, 9999999, EDATE, 'Enter end date') ETIME = 0 ETIME = GETNUM( 0, 9999999, ETIME, 'Enter end time') NSITES = 0 NSITES = GETNUM( 0, 999, NSITES, 'Enter # of Cl I Areas in grid') NSRC = 0 NSRC = GETNUM( 0, 9999999, NSRC, 'Enter # of source tags') NSPC = 0 NSPC = GETNUM( 0, 9999999, NSPC, 'Enter # of chemical species') NMON = 0 NMON = GETNUM( 0, 9999999, NMON, 'Enter # months for RH factors') IF ( .NOT. OPEN3( INFILE, FSREAD3, PNAME ) ) THEN MESG = 'Could not open input file ' // INFILE CALL M3EXIT( PNAME, 0, 0, MESG, 2 ) END IF IF ( .NOT. DESC3( INFILE ) ) THEN MESG = 'Could not DESC3 file ' // INFILE CALL M3EXIT( PNAME, 0, 0, MESG, 2 ) END IF NROWS = NROWS3D NCOLS = NCOLS3D NLAYS = NLAYS3D NVARS = NVARS3D GDTYP = GDTYP3D TSTEP = TSTEP3D DX = XCELL3D DY = YCELL3D XORIG = XORIG3D YORIG = YORIG3D C set up lambert grid ANG1 = P_ALP3D ! lower latitude of projection cone ANG2 = P_BET3D ! upper latitude of projection cone ANG3 = P_GAM3D ! central meridian of projection XCEN = XCENT3D ! origin lon YCEN = YCENT3D ! origin lat IF (.NOT. SETLAM (ANG1,ANG2,ANG3,XCEN,YCEN) ) then CALL M3EXIT (PNAME, 0, 0, & 'Lambert projection setup error', 2 ) end if ! if SETLAM failed SDATE3D = STDATE STIME3D = STTIME NSTEP = MXREC3D NSECS = TIME2SEC(TSTEP) NSTEP = SECSDIFF(STDATE, STTIME, EDATE, ETIME)/NSECS + 1 C.....We are borrowing description from INFILE IF ( .NOT. OPEN3( OUTFILE, FSNEW3, PNAME ) ) THEN MESG = 'Could not open input file ' // OUTFILE CALL M3EXIT( PNAME, 0, 0, MESG, 2 ) END IF C...............Allocate buffers ALLOCATE( $ MASK ( NCOLS, NROWS ), RAY ( NCOLS, NROWS ), & A ( NCOLS, NROWS, NLAYS, NVARS ), & AIN ( NCOLS, NROWS, NLAYS, NVARS ), $ OUT ( NCOLS, NROWS, NVARS, NSTEP ), $ FSRH ( NCOLS, NROWS, NMON ), FLRH ( NCOLS, NROWS, NMON ), $ FSSRH( NCOLS, NROWS, NMON ), FRH ( NCOLS, NROWS, NMON ), $ SCNAME( NSRC ), SPNAME( NSPC ), SITENAME( NSITES ), $ COLST( NSITES ), ROWST( NSITES ), $ STAT = STATUS ) IF ( STATUS .NE. 0 ) THEN WRITE( MESG, '( A, I10)' ) & 'Memory allocation failed: STAT=', STATUS CALL M3EXIT( PNAME, 0, 0, MESG, 2 ) END IF JDATE = STDATE JTIME = STTIME C Initialize MASK and other variables DO IMON = 1, NMON DO R = 1, NROWS DO C = 1, NCOLS MASK (C, R) = .FALSE. RAY (C, R) = 0. FSRH (C, R, IMON) = 0. FLRH (C, R, IMON) = 0. FSSRH(C, R, IMON) = 0. FRH (C, R, IMON) = 0. END DO END DO END DO C Read IMPROVE site lat-lon (csv file) and convert to LCC coordinates; set site C mask and site-specific Rayleigh (annual) extinction; populate monthly RH adjustment factor arrays DO ISITE = 1, NSITES FIRSTIME = .TRUE. COLST( ISITE ) = 0 ROWST( ISITE ) = 0 DO IMON = 1, NMON READ(SITEUNIT,1000) LINE CALL PARSLINE( LINE, MXTCOL, SEGMENT ) ! break lines into segments IF (FIRSTIME) THEN SITENAME( ISITE ) = SEGMENT( 1 ) LAT = STR2REAL( SEGMENT( 2 ) ) LON = STR2REAL( SEGMENT( 3 ) ) IF (.NOT. LL2LAM (LON, LAT, X, Y) ) then WRITE(*,*) ' LAMBERT CONVERSION ERROR ' STOP END IF C = INT( (X-XORIG)/DX ) + 1 R = INT( (Y-YORIG)/DY ) + 1 IF ( (C .LT. 1 .OR. C .GT. NCOLS) .OR. $ (R .LT. 1 .OR. R .GT. NROWS) ) THEN WRITE ( MESG, '(A, A5, A)' ), $ 'LL2LAM failed to find site named '//SITENAME(ISITE)// $ ' in domain' CALL M3WARN( PNAME,JDATE, JTIME, MESG ) ELSE ! set mask variable for this location MASK (C, R) = .TRUE. RAY (C, R) = STR2REAL( SEGMENT(4) ) COLST( ISITE ) = C ROWST( ISITE ) = R FIRSTIME = .FALSE. WRITE(*,*) COLST(ISITE), ROWST(ISITE), SITENAME(ISITE) ENDIF ENDIF IF ( MASK(C, R) ) THEN ! if site found set RH adjustment factor arrays FSRH (C, R, IMON) = STR2REAL( SEGMENT(5) ) FLRH (C, R, IMON) = STR2REAL( SEGMENT(6) ) FSSRH(C, R, IMON) = STR2REAL( SEGMENT(7) ) FRH (C, R, IMON) = STR2REAL( SEGMENT(8) ) ENDIF END DO ! end IMON loop for this site END DO ! end read on sites DEFAULT = 'ASCIIFILE' CALL GETSTR('ENTER NAME FOR OUTPUT ASCII DATA FILE ', $ DEFAULT, ASCIIFILE) ASCIIUNIT = GETEFILE(ASCIIFILE, .FALSE., .TRUE., PNAME) IF (ASCIIUNIT .LT. 0) THEN CALL M3EXIT (PNAME, 0, 0, & 'ASCIIFILE FILE OPENING ERROR', 2 ) END IF ! if output failed WRITE(FMTHDR, '(A,I2.2,A)') '(1X,A6,A5,8A6,A8,',NSPC,'A14)' WRITE(FMTDATA,'(A,I2.2,A)') $ '(1X,A6,A5,F6.0,3I6,F6.2,F6.1,F6.2,F6.1,A10,',NSPC,'G14.6)' C Extract source tags from tagged tracer names DO S = 1, NSRC ! loop on tagged sources WRITE (STR, '(A1,I3.3)') 'S',S ! # 1 = all non-tagged sources SCNAME(S) = STR END DO C Extract chemical species names from tagged tracer names DO S = 1, NSPC ! loop on chemical species V = (S - 1) * NSRC + 1 STR = VNAME3D(V) IF (NSRC .GT. 1) THEN SPNAME(S) = STR(1:3) ! for CAMx sa file species names ELSE SPNAME(S) = STR ! for CAMX avrg file species names ENDIF END DO C Write ASCII file header WRITE (ASCIIUNIT, FMTHDR) 'Tracer', ' site', $ ' ssray', ' mon ', ' day ', ' year ', ' fsrh ', ' flrh ', $ 'fssrh ', ' frh ', ' date ' , (SPNAME(S),S=1,NSPC) C Time step loop over desired time interval DO T = 1, NSTEP C Read and accumulate IF ( .NOT.READ3( INFILE, ALLVAR3, ALLAYS3, $ JDATE, JTIME, AIN)) THEN MESG = 'Could not read all variables from ' // INFILE CALL M3EXIT( PNAME, JDATE, JTIME, MESG, 2 ) END IF C Apply the site mask to select model output at IMPROVE sites DO S = 1, NVARS DO L = 1, NLAYS DO R = 1, NROWS DO C = 1, NCOLS IF ( MASK( C, R ) ) THEN A( C,R,L,S ) = AIN( C,R,L,S ) IF (L .EQ. 1) OUT( C,R,S,T ) = A( C,R,L,S ) ! for ascii ELSE A(C,R,L,S) = BADVAL3 IF (L .EQ. 1) OUT(C, R, S, T) = 0. END IF END DO ! if C > NCOLS END DO ! if R > NROWS END DO ! if L > NLAYS END DO ! if S > NVARS C Write output IF ( .NOT. WRITE3( OUTFILE, ALLVAR3, JDATE, JTIME, A ) ) THEN MESG = 'Could not write all vars to file ' // OUTFILE CALL M3EXIT( PNAME, 0, 0, MESG, 2 ) END IF CALL NEXTIME(JDATE, JTIME, TSTEP) END DO ! if T > NSTEP C Write output for all species each day at each Class I area within the domain C for each tagged source (need time step loop to be innermost) DO S = 1, NSRC ! loop on source tags NVMAX = (NSPC - 1) * NSRC + S DO ISITE = 1, NSITES C = COLST( ISITE ) R = ROWST( ISITE ) IF ( MASK (C, R) ) THEN JDATE = STDATE ! initialize date/time JTIME = STTIME ! " " DO T = 1, NSTEP YEAR = INT(JDATE/1000.) IYEAR = YEAR - 2000 ! last two digits of yr CALL DAYMON(JDATE, IMON, IDAY ) WRITE (DATE,'(I4.2,A1,I2.2,A1,I2.2)') $ IMON, SEP, IDAY, SEP, IYEAR ! MM/DD/YY format WRITE (ASCIIUNIT,FMTDATA) SCNAME(S), SITENAME(ISITE), $ RAY(C, R), IMON, IDAY, YEAR, $ FSRH(C, R, IMON), FLRH(C, R, IMON), $ FSSRH(C, R, IMON), FRH(C, R, IMON), $ DATE, (OUT(C,R,V,T), V = S, NVMAX, NSRC) CALL NEXTIME(JDATE, JTIME, TSTEP) END DO ! if T > NSTEP ENDIF ! skip if monitor site not in domain END DO ! if ISITE > NSITES END DO ! if S > max source tag CLOSE (asciiunit) CALL M3EXIT( PNAME, 0, 0, 'Normal completion', 0 ) 1000 FORMAT( 40A ) STOP END C ========================================================================== C ========================================================================== C ========================================================================== SUBROUTINE PARSLINE( LINE, N, SEGMENT ) C*********************************************************************** C subroutine body starts at line C C DESCRIPTION: C This subroutine separates a "list-formatted" line of strings in which C the segments may or may not have quotes. Although fortran requires C the quotes for true list-formatting, this subroutine can be used when C the quotes are only present to enclose a character (such as space, comma, C or semi-colon) that would otherwise be a delimiter. If an "!" is C encountered, everything after it is treated as a comment. C C PRECONDITIONS REQUIRED: C C SUBROUTINES AND FUNCTIONS CALLED: C C REVISION HISTORY: C Created by M. Houyoux 3/99 C C****************************************************************************/ C C Project Title: Sparse Matrix Operator Kernel Emissions (SMOKE) Modeling C System C File: @(#)$Id: parsline.f,v 1.10 2003/05/02 13:46:31 cas Exp $ C C COPYRIGHT (C) 2000, MCNC--North Carolina Supercomputing Center C All Rights Reserved C C See file COPYRIGHT for conditions of use. C C Environmental Programs Group C MCNC--North Carolina Supercomputing Center C P.O. Box 12889 C Research Triangle Park, NC 27709-2889 C C env_progs@mcnc.org C C Pathname: $Source: /afs/isis/depts/cep/emc/apps/archive/edss_tools/edss_tools/src/lib/parsline.f,v $ C Last updated: $Date: 2003/05/02 13:46:31 $ C C*************************************************************************** IMPLICIT NONE C........... EXTERNAL FUNCTIONS CHARACTER*2 CRLF INTEGER FINDC EXTERNAL CRLF, FINDC C........... SUBROUTINE ARGUMENTS CHARACTER(*), INTENT (IN) :: LINE ! character string to parse INTEGER , INTENT (IN) :: N ! maximum array length CHARACTER(*), INTENT(OUT) :: SEGMENT( N ) ! parsed string C........... Local parameters INTEGER , PARAMETER :: NDELIM = 4 CHARACTER*1, PARAMETER :: DELIMLST( NDELIM ) = & (/ ',', ' ', ';', ' ' /) C........... Array of 1-char strings for processing CHARACTER*1 ARRSTR( 5120 ) ! 256 * 20 C........... Arrays for sorting non-delimiters on a per-machine basis INTEGER NDINDX ( NDELIM ) CHARACTER*1, SAVE :: DELIMSRT( NDELIM ) C........... Other local variables INTEGER I, J, L, L1, L2 ! counters and indices INTEGER IXP ! index to non-delimeters INTEGER :: NCNT ! count of fields LOGICAL :: ALPHA ! true when within alpha-numeric LOGICAL :: DELIM ! true when within or past delimiter LOGICAL, SAVE:: FIRSTIME = .TRUE.! true first time routine is called LOGICAL :: PREVDELIM = .TRUE. ! true when last char was a delim LOGICAL :: NUMBER ! true when within number in string LOGICAL :: QUOTED ! true when within quotes in string LOGICAL :: THISNMBR ! true when current iteration is numbr CHARACTER*1 CBUF ! temporary buffer CHARACTER*1 :: DOUBLEQ = '"' CHARACTER*1 :: SINGLEQ = "'" CHARACTER*1 :: PERIOD = '.' CHARACTER*1 QUOTVAL ! value of starting quote CHARACTER*300 MESG ! message buffer CHARACTER*16 :: PROGNAME = 'PARSLINE' ! program name C*********************************************************************** C begin body of subroutine PARSLINE C......... The first time the routine is called, sort the list of delimiters IF( FIRSTIME ) THEN DO I = 1, NDELIM NDINDX( I ) = I END DO CALL SORTIC( NDELIM, NDINDX, DELIMLST ) DO I = 1, NDELIM J = NDINDX( I ) DELIMSRT( I ) = DELIMLST( J ) END DO FIRSTIME = .FALSE. END IF L2 = LEN_TRIM( LINE ) C......... Check for comments, and use to set the end of the line L = INDEX( LINE( 1:L2 ), '!' ) IF( L .LE. 0 ) THEN L = L2 ELSE L = L - 1 END IF C......... Skip blank lines IF( L .EQ. 0 ) RETURN C......... Initialize count, flags, and segments (npte, initializing in C the variable definitions is insufficient) NCNT = 0 SEGMENT = ' ' ! array ALPHA = .FALSE. DELIM = .TRUE. NUMBER = .FALSE. QUOTED = .FALSE. C......... Process LINE 1-character at a time DO I = 1, L CBUF = LINE( I:I ) C............. Look for character in delimiters IXP = FINDC( CBUF, NDELIM, DELIMSRT ) C............. Evaluate the current character for number or not THISNMBR = ( CBUF .GE. '0' .AND. CBUF .LE. '9' ) C............. Waiting for next field... IF( DELIM ) THEN NUMBER = THISNMBR ALPHA = ( .NOT. NUMBER .AND. IXP .LE. 0 ) IF( CBUF .EQ. SINGLEQ ) THEN QUOTED = .TRUE. DELIM = .FALSE. QUOTVAL = SINGLEQ PREVDELIM = .FALSE. L1 = I + 1 NCNT = NCNT + 1 ELSE IF( CBUF .EQ. DOUBLEQ ) THEN QUOTED = .TRUE. DELIM = .FALSE. QUOTVAL = DOUBLEQ PREVDELIM = .FALSE. L1 = I + 1 NCNT = NCNT + 1 ELSE IF( ALPHA ) THEN DELIM = .FALSE. PREVDELIM = .FALSE. L1 = I NCNT = NCNT + 1 ELSE IF( NUMBER ) THEN DELIM = .FALSE. PREVDELIM = .FALSE. L1 = I NCNT = NCNT + 1 C............... If another delimeter, then another field, but last C field was blank UNLESS delim is a space ELSE IF( CBUF .NE. DELIMLST( 2 ) ) THEN IF( PREVDELIM ) THEN NCNT = NCNT + 1 ELSE PREVDELIM = .TRUE. END IF END IF ! Else its a space delimiter C............. In a quoted field, skip everything unless it is an end quote ELSE IF( QUOTED ) THEN IF( CBUF .EQ. QUOTVAL ) THEN QUOTED = .FALSE. DELIM = .TRUE. PREVDELIM = .FALSE. L2 = I - 1 CALL STORE_SEGMENT END IF C............. If start of field was a number, but adjacent character is not C a delimiter, then turn field into an alpha ELSE IF( NUMBER .AND. .NOT. THISNMBR .AND. IXP .LE. 0 ) THEN ALPHA = .TRUE. NUMBER = .FALSE. C............. If start of field was a number or alpha, and this is a C delimiter, then end of number has been reached ELSE IF( IXP .GT. 0 ) THEN ALPHA = .FALSE. NUMBER = .FALSE. DELIM = .TRUE. PREVDELIM = .TRUE. L2 = I - 1 CALL STORE_SEGMENT END IF END DO C......... Store final segment IF( CBUF .EQ. QUOTVAL ) L = L - 1 L2 = L IF( IXP .LE. 0 ) CALL STORE_SEGMENT RETURN C****************** FORMAT STATEMENTS ************************************ C........... Internal buffering formats............ 94xxx 94010 FORMAT( 10( A, :, I8, :, 1X ) ) C****************** INTERNAL SUBPROGRAMS ***************************** CONTAINS C............. This subprogram stores the segment from the input string SUBROUTINE STORE_SEGMENT IF( NCNT .LE. N ) THEN SEGMENT( NCNT ) = ADJUSTL( LINE( L1:L2 ) ) ELSE MESG = 'ERROR: Overflow prevented while '// & 'parsing line ' // PROGNAME CALL M3MSG2( MESG ) MESG = 'First 200 characters of line contents are:' CALL M3MSG2( MESG ) MESG = LINE( 1:200 ) CALL M3MSG2( MESG ) MESG = 'Formatting problem.' CALL M3EXIT( PROGNAME, 0, 0, MESG, 2 ) END IF END SUBROUTINE STORE_SEGMENT END SUBROUTINE PARSLINE