SUBROUTINE SETUPCOLDEN( COLDEN, NDEN ) C ------------------------------------------------------------------- C PURPOSE: This subroutine sets the column densities up for C generating g-factor profiles. C INPUT: INTEGER NDEN, the number elements in COLDEN C OUTPUT: REAL COLDEN(NDEN), column densities C CALLED BY: GFACDRIVER C HISTORY: March 1995 by S. J. Bertin, CPI C Adapted from IDL version. C ------------------------------------------------------------------- INTEGER NDEN REAL COLDEN(NDEN) INTEGER I DO 10 I=1,NDEN COLDEN(I) = 10.0**(10.0+(I/(NDEN/10.0))) 10 CONTINUE END SUBROUTINE LOADCOEFF( FILENAME ) C ------------------------------------------------------------------- C PURPOSE: This subroutine reads a g-factor database file and C sets the variables in the GFACTORS common block. C INPUTS: CHARACTER*(*) FILENAME - name of g-factor database file C CALLED BY: GFACDRIVER C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA C Adapted from IDL version. C ------------------------------------------------------------------- CHARACTER*(*) FILENAME INCLUDE 'params_gfac.inc' INCLUDE 'gfactors.inc' INTEGER I, IF10P7, ISZA, IFEAT CHARACTER*80 JUNK IF (FILENAME.NE.GFACFILENAME) THEN GFACFILENAME='' OPEN(UNIT=COEFFUNIT, FILE=FILENAME, STATUS='OLD', ERR=9010) DO 10 I=1,3 READ(COEFFUNIT,6000,ERR=9000) JUNK 10 CONTINUE READ(COEFFUNIT,5000,ERR=9000) NFEAT DO 20 I=1,NFEAT+1 READ(COEFFUNIT,6000,ERR=9000) JUNK 20 CONTINUE READ(COEFFUNIT,5000,ERR=9000) NF10P7 READ(COEFFUNIT,*,ERR=9000) (F10P7(I),I=1,NF10P7) READ(COEFFUNIT,6000,ERR=9000) JUNK READ(COEFFUNIT,5000,ERR=9000) NSZA READ(COEFFUNIT,*,ERR=9000) (SZA(I),I=1,NSZA) READ(COEFFUNIT,6000,ERR=9000) JUNK READ(COEFFUNIT,5000,ERR=9000) NCOEFF DO 100 IF10P7=1,NF10P7 DO 110 ISZA=1,NSZA DO 30 I=1,3 READ(COEFFUNIT,6000,ERR=9000) JUNK 30 CONTINUE DO 120 IFEAT=1,NFEAT READ(COEFFUNIT,6000,ERR=9000) JUNK READ(COEFFUNIT,*,ERR=9000) + (COEFF(I,IFEAT,IF10P7,ISZA),I=1,NCOEFF) 120 CONTINUE 110 CONTINUE 100 CONTINUE GFACFILENAME=FILENAME ENDIF CLOSE(UNIT=COEFFUNIT) GOTO 9999 5000 FORMAT(6X,I6) 6000 FORMAT(A) C ERROR HANDLING 9000 CLOSE(UNIT=COEFFUNIT) 9010 PRINT*, 'ERROR OPENING OR READING FROM ',FILENAME 9999 CONTINUE END SUBROUTINE CHECKRANGES( FEATURE ) C ------------------------------------------------------------------- C PURPOSE: This subroutine checks the feature to verify that it is C within the valid range. C INPUT: INTEGER FEATURE - Feature to check. C OUTPUT: Error message printed and FEATURE set to one if out of C range. C CALLED BY: GFACDRIVER C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from IDL version. C ------------------------------------------------------------------- INTEGER FEATURE INCLUDE 'params_gfac.inc' INCLUDE 'gfactors.inc' IF(FEATURE.LT.1 .OR. FEATURE.GT.NFEAT) THEN PRINT*,'ERROR: FEATURE ',FEATURE,' IS INVALID.' PRINT*,' MUST BE BETWEEN 1 AND ',NFEAT,' INCLUSIVE.' PRINT*,' SETTING FEATURE TO 1' FEATURE = 1 ENDIF END INTEGER FUNCTION FIND_NEAREST( POINT, GRID, NPOINTS ) C ------------------------------------------------------------------- C PURPOSE: This function finds the index in GRID whose value C is closest to POINT C INPUT: REAL POINT - value to find. C REAL GRID(NPOINTS) - array of possible values. C INTEGER NPOINTS - number of points in grid. C OUTPUT: INTEGER - the index of the nearest point. C CALLED BY: SMARTINTERPOL C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C ------------------------------------------------------------------- INTEGER NPOINTS REAL POINT, GRID(NPOINTS) INTEGER NEAREST NEAREST = 0 10 NEAREST = NEAREST+1 IF((POINT.GT.GRID(NEAREST)) .AND. + (NEAREST.LT.(NPOINTS-1))) GOTO 10 IF((GRID(NEAREST+1)-POINT).LT.(POINT-GRID(NEAREST))) THEN NEAREST = NEAREST+1 ENDIF FIND_NEAREST = NEAREST END SUBROUTINE SMARTINTERPOL( XPOINT, YPOINT, XGRID, YGRID, ZDATA, + ZPOINT, NX, NY, RZX, RZY ) C ------------------------------------------------------------------- C PURPOSE: This subroutine performs a bilinear C interpolation/extrapolation. C INPUT: REAL XPOINT - the coordinate on the x grid where Z is C to be determined. C REAL YPOINT - the coordinate on the y grid where Z is C to be determined. C REAL XGRID(NX) - the x grid. C REAL YGRID(NY) - the y grid. C REAL ZDATA(RZX, RZY) - the two dimensional data array C INTEGER NX - Number of points in XGRID. C INTEGER NY - Number of points in YGRID. C INTEGER RZX - First dimension of ZDATA. C INTEGER RZY - Second dimension of ZDATA. C OUTPUT: REAL ZPOINT - the value of Z at the specified point C COMMENTS: This procedure will extrapolate if the specified point is outside C the X or Y grids. The user should check his data before C calling smartinterpol if extrapolation is not desired. C CALLED BY: COEFFINTERP C ROUTINES CALLED: C FIND_NEAREST C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA C ------------------------------------------------------------------- INTEGER NX, NY, RZX, RZY REAL XPOINT, YPOINT, XGRID(NX), YGRID(NY), ZDATA(RZX,RZY), ZPOINT INTEGER XI, YI, FIND_NEAREST REAL M1, M2 XI = FIND_NEAREST(XPOINT, XGRID, NX) YI = FIND_NEAREST(YPOINT, YGRID, NY) IF(XI.EQ.NX .OR. XPOINT.LT.XGRID(XI)) XI = XI-1 IF(YI.EQ.NY .OR. YPOINT.LT.YGRID(YI)) YI = YI-1 M1 = (XPOINT-XGRID(XI))/(XGRID(XI+1)-XGRID(XI)) M2 = (YPOINT-YGRID(YI))/(YGRID(YI+1)-YGRID(YI)) ZPOINT = (1.0-M1) * (1.0-M2) * ZDATA(XI , YI ) + + M1 * (1.0-M2) * ZDATA(XI+1, YI ) + + M1 * M2 * ZDATA(XI+1, YI+1) + + (1.0-M1) * M2 * ZDATA(XI , YI+1) END SUBROUTINE COEFFINTERP( ASZA, AF10P7, FEATURE, COEFFOUTPUT ) C ------------------------------------------------------------------- C PURPOSE: This subroutine interpolates on the SZA and F10P7 grids C to get coefficients at the points specified in the C SZA and F10P7 arrays. C INPUT: REAL ASZA - solar zenith angles at which the g-factor C profiles is desired C REAL AF10P7 - corresponding 10.7 cm flux value C INTEGER FEATURE - feature to use for profile. C OUTPUT: REAL COEFFOUTPUT(MNCOEFF) - the coefficients at the C selected SZA and F10P7. C CALLED BY: GFACDRIVER C ROUTINES CALLED: C SMARTINTERPOL C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from IDL Version. C ------------------------------------------------------------------- INCLUDE 'params_gfac.inc' INCLUDE 'gfactors.inc' INTEGER FEATURE REAL ASZA, AF10P7, COEFFOUTPUT(MNCOEFF) INTEGER ICOEFF, ISZA, IF10P7 REAL TCOEFF(MNF10P7,MNSZA) DO 10 ICOEFF=1,NCOEFF DO 20 IF10P7=1,NF10P7 DO 30 ISZA=1,NSZA TCOEFF(IF10P7,ISZA) = + COEFF(ICOEFF,FEATURE,IF10P7,ISZA) 30 CONTINUE 20 CONTINUE CALL SMARTINTERPOL( AF10P7, ASZA, F10P7(1:NF10P7), + SZA(1:NSZA), + TCOEFF, COEFFOUTPUT(ICOEFF), + NF10P7, NSZA, MNF10P7, MNSZA) 10 CONTINUE END SUBROUTINE BANKS( COLDEN, NDEN, COEFF, GFACTOR ) C ------------------------------------------------------------------- C PURPOSE: This subroutine generates gfactor profiles using the C coefficients with the following function: C C gfactor = A * ( 2.0 - ( R + 1 ) * exp ( -C * ColDen ) ) * C exp ( - D * ColDen ) C C A = Coeff(1) C R = Coeff(2) C C = Coeff(3) C D = Coeff(4) C C INPUTS: REAL COLDEN(NDEN) - column density data. C INTEGER NDEN - number of column densities. C REAL COEFF(4) - array of Coefficients for the above C function C OUTPUTS: REAL GFACTOR(NDEN) - g-factors corresponding to given C column densities. C CALLED BY: BANKS_GAUSSIAN C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from IDL version. C ------------------------------------------------------------------- INTEGER NDEN REAL COLDEN(NDEN), COEFF(4), GFACTOR(NDEN) INTEGER I DO 10 I=1,NDEN GFACTOR(I) = COEFF(1) * ( 2.0 - ( COEFF(2) + 1.0 ) * + EXP( -COEFF(3) * COLDEN(I) ) ) * + EXP( -COEFF(4) * COLDEN(I) ) 10 CONTINUE END SUBROUTINE GAUSSIAN( COLDEN, NDEN, COEFF, GFACTOR ) C ------------------------------------------------------------------- C PURPOSE: This procedure generates gfactor using the Coefficients C provided with the following function: C C gfactor = E * exp ( - ((alog(ColDen) - alog(F))/alog(G))^2 ) C C E = Coeff(1) C F = Coeff(2) C G = Coeff(3) C C INPUTS: REAL COLDEN(NDEN) - column density data C REAL NDEN - number of column densities C REAL COEFF(3) - array of coefficients for the above C function C OUTPUTS: REAL GFACTOR(NDEN) - g-factors corresponding to given C column densities. C CALLED BY: BANKS_GAUSSIAN C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from IDL version. C ------------------------------------------------------------------- INTEGER NDEN REAL COLDEN(NDEN), COEFF(3), GFACTOR(NDEN) INTEGER I DO 10 I=1,NDEN GFACTOR(I) = COEFF(1) * EXP( -( (LOG(COLDEN(I)) - + LOG(COEFF(2)) ) / LOG(COEFF(3)) )**2 ) 10 CONTINUE END SUBROUTINE BANKS_GAUSSIAN( COLDEN, NDEN, COEFF, GFACTOR ) C ------------------------------------------------------------------- C PURPOSE: This procedure generates g-factor profiles using the C coefficients with the following function: C C gfactor = A * ( 2.0 - ( R + 1 ) * exp ( -C * ColDen ) ) * C exp ( - D * ColDen ) + C E * exp ( -( (alog(ColDen) - alog(F)) / alog(G) )^2 ) + C T * exp ( -( (alog(ColDen) - alog(U)) / alog(V) )^2 ) + C A = Coeff(1) C R = Coeff(2) C C = Coeff(3) C D = Coeff(4) C E = Coeff(5) C F = Coeff(6) C G = Coeff(7) C T = Coeff(8) C U = Coeff(9) C V = Coeff(10) C C INPUTS: REAL COLDEN(NDEN) - column density data C REAL NDEN - number of column densities C REAL COEFF(3) - array of coefficients for the above C function C OUTPUTS: REAL GFACTOR(NDEN) - g-factors corresponding to given C column densities. C CALLED BY: GFACDRIVER C ROUTINES CALLED: C BANKS - function for generating Banks portion of fit C GAUSSIAN - function for generating Gaussian portion of fit C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from IDL version. C ------------------------------------------------------------------- INCLUDE 'params_gfac.inc' INTEGER NDEN REAL COLDEN(NDEN), COEFF(MNCOEFF), GFACTOR(NDEN) INTEGER I REAL TGFAC(MNDEN) REAL BANKSCOEFF(4), GAUSS1COEFF(3), GAUSS2COEFF(3) DO 3 I=1,4 BANKSCOEFF(I) = COEFF(I) 3 CONTINUE DO 4 I=1,3 GAUSS1COEFF(I) = COEFF(I+4) 4 CONTINUE DO 5 I=1,3 GAUSS2COEFF(I) = COEFF(I+7) 5 CONTINUE CALL BANKS(COLDEN, NDEN, BANKSCOEFF, TGFAC(1:NDEN)) DO 10 I=1,NDEN GFACTOR(I) = TGFAC(I) 10 CONTINUE CALL GAUSSIAN(COLDEN, NDEN, GAUSS1COEFF, TGFAC(1:NDEN)) DO 20 I=1,NDEN GFACTOR(I) = GFACTOR(I) + TGFAC(I) 20 CONTINUE CALL GAUSSIAN(COLDEN, NDEN, GAUSS2COEFF, TGFAC(1:NDEN)) DO 30 I=1,NDEN GFACTOR(I) = GFACTOR(I) + TGFAC(I) 30 CONTINUE END SUBROUTINE GFACDRIVER( FILENAME, SZA, F10P7, GFACTOR, COLDEN, + FEATURE, NDEN ) C ------------------------------------------------------------------- C PURPOSE: This procedure loads a gfactor coefficient database, C interpolates to get coefficients at the specified C SZA and F10P7, and generates gfactor profiles. C INPUT: CHARACTER*(*) FILENAME - the name of the coefficient C database file C REAL SZA - solar zenith angle at which the g-factor C profile is desired C REAL F10P7 - corresponding 10.7 cm flux value C REAL COLDEN(NDEN) - this variable may be used to C specify the column densities to use. C INTEGER FEATURE - the feature used for the profile. C INTEGER NDEN - the nember of column densities to use. C OUTPUT: REAL COLDEN(NDEN) - this variable returns the generated C column densities from the routine if COLDEN(1) is C -1.0 on input. C REAL GFACTOR(NDEN) - gfactors for each case specified. C ROUTINES CALLED: C SETUPCOLDEN, LOADCOEFF, CHECKRANGES, COEFFINTERP, C BANKS_GAUSSIAN C HISTORY: March 1995, S. J. Bertin, CPI, Fairfax, VA. C Adapted from the IDL version. C ------------------------------------------------------------------- CHARACTER*(*) FILENAME INTEGER FEATURE, NDEN REAL SZA, F10P7, GFACTOR(NDEN), COLDEN(NDEN) INCLUDE 'params_gfac.inc' REAL NEWCOEFF(MNCOEFF) IF(NDEN.GT.MNDEN) THEN PRINT*, 'ERROR: NDEN TOO LARGE, SETTING TO ', MNDEN NDEN = MNDEN ENDIF IF(COLDEN(1).EQ.-1) THEN CALL SETUPCOLDEN(COLDEN(1:NDEN), NDEN) ENDIF CALL LOADCOEFF( FILENAME ) CALL CHECKRANGES( FEATURE ) CALL COEFFINTERP( SZA, F10P7, FEATURE, NEWCOEFF ) CALL BANKS_GAUSSIAN( COLDEN(1:NDEN), NDEN, NEWCOEFF, + GFACTOR(1:NDEN) ) END