c*********************************************************************** c Djangoh radiative correction results container * C DELTA=d Sigma (Born+RC) / d Sigma (Born) - 1D0 * c Last update: 14.04.2009 (minor changes) * c * c Entries: CALL DJ_RES(X,Y,DELTA,ERR) * c Entries: CALL DJ_RES_ADV(X,Y,DELTA,ERR_STA,ERR_SYS) * c Input: X,Y * c Output: DELTA,ERR(ERR_STA,ERR_SYS) * c * c Settings common block: * c COMMON/DJ_RES_SETTS/NLOG,LOGUNIT,DJ_RES_ERROR * c >NLOG: number of calls to be logged. Default:0. * c >LOGUNIT: output unit to print logs/errors. Default:42 * c >DJ_RES_UNKNOWN: 'unknown' return value. Default:1D6 * c >DJ_RES_ERRUNKNOWN: 'unknown' error value. Default:1D6 * c Do not change DJ_RES_UNKNOWN after initialization (first call)! * c * c Uses djangoh_r.dat values created by djangoh_u_r.f user routine. * c. File format * c. #(int) IVER: file version (reserved, currently 2) * c. #(int) I2ERRS: 1(0) the stat. & sys. errors are separated * c. #(int) S: S-Mp2-Ml2 * c. #(int) NPTSX,NPTSY,NMODEX,NMODEY: number of points and grid mode * c. #(float) XPTMIN,XPTMAX,YPTMIN,YPTMAX: min/max point values * c. #(int) DELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY): results array * c. #(int) NCTRLSUM: currently NPTSX*NPTSY * c Initialization is performed at first call or using DJ_RES_INIT. * c * c To check initialization: * c 1. Use IDJRES_INITOK() function (returns 1 or 0) * c 2. Or check (NSTATUS.EQ.1) variable in common/DJRES_DATA/ * c NSTATUS: 0(not inited), 1(ok), 2(error) * c * c See webpage for details: http://cern.ch/~makarenko/dj_res/ * c * c Functions/subroutines list: * c DJ_RES: Get relative correction and error * c DJ_RES_ADV: Get relative correction end stat./sys. errors. * c DJ_RES_INIT: Loads data file. Invoked automatically at first call. * c IDJRES_INITOK: Check initialization (returns 1 or 0). * c*********************************************************************** C. Function settings block BLOCK DATA DJ_RES_SETTS INTEGER NLOG,LOGUNIT REAL*8 DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN COMMON/DJRES_SETS/NLOG,LOGUNIT,DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN DATA LOGUNIT,NLOG/42,0/ DATA DJ_RES_UNKNOWN/1.0D6/ DATA DJ_RES_ERRUNKNOWN/1.0D6/ END C. Data block C. RPTSX,RPTSY: data points arrays C. NPTSX,NPTSY,NPTS_MAX: data table size C. XPTMIN,XPTMAX,YPTMIN,YPTMAX: data table margins C. NMODEX,NMODEY: data grid modes (0: linear, 1: logarythmic) C. NSTATUS: initialization flag (0:no data read, 1:ready, 2:error) C. NCTRLSUM: table control sum C. GSP: S-Mp2-Ml2 BLOCK DATA DJ_RES_DATA INTEGER NPTS_MAX,NPTS_MAX2,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,GSP,RESDELTA & ,ERRSTA,ERRSYS COMMON/DJRES_DATA/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,GSP & ,RESDELTA(NPTS_MAX,NPTS_MAX),ERRSTA(NPTS_MAX,NPTS_MAX) & ,ERRSYS(NPTS_MAX,NPTS_MAX),NSTATUS,NCTRLSUM C. Default init DATA NPTSX,NPTSY/0,0/ DATA RPTSX/NPTS_MAX*0D0/ DATA RPTSY/NPTS_MAX*0D0/ DATA NSTATUS,NCTRLSUM/0,0/ DATA XPTMIN,XPTMAX,YPTMIN,YPTMAX/0D0,0D0,0D0,0D0/ END C. ******************************************************************** C. DJ_RES: Get relative correction and error C. DELTA=d Sigma (Born+RC) / d Sigma (Born) - 1D0 C. ******************************************************************** SUBROUTINE DJ_RES(X,Y,DELTA,ERR) IMPLICIT NONE REAL*8 X,Y,DELTA,ERR,ERR_STA,ERR_SYS C. Function settings INTEGER NLOG,LOGUNIT REAL*8 DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN COMMON/DJRES_SETS/NLOG,LOGUNIT,DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN DELTA=DJ_RES_UNKNOWN ERR_STA=0D0 ERR_SYS=DJ_RES_ERRUNKNOWN CALL DJ_RES_ADV(X,Y,DELTA,ERR_STA,ERR_SYS) IF((ERR_STA.LT.DJ_RES_ERRUNKNOWN).AND & .(ERR_SYS.LT.DJ_RES_ERRUNKNOWN))THEN ERR=DSQRT(ERR_STA*ERR_STA+ERR_SYS*ERR_SYS) ELSE ERR=DJ_RES_ERRUNKNOWN ENDIF RETURN END C. ******************************************************************** C. DJ_RES_ADD: Get relative correction and error C. DELTA=d Sigma (Born+RC) / d Sigma (Born) - 1D0 C. ******************************************************************** SUBROUTINE DJ_RES_ADV(X,Y,DELTA,ERR_STA,ERR_SYS) IMPLICIT NONE REAL*8 X,Y,DELTA,ERR_STA,ERR_SYS C. Function settings INTEGER NLOG,LOGUNIT REAL*8 DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN COMMON/DJRES_SETS/NLOG,LOGUNIT,DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN C. Data block INTEGER NPTS_MAX,NPTS_MAX2,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,GSP,RESDELTA & ,ERRSTA,ERRSYS COMMON/DJRES_DATA/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,GSP & ,RESDELTA(NPTS_MAX,NPTS_MAX),ERRSTA(NPTS_MAX,NPTS_MAX) & ,ERRSYS(NPTS_MAX,NPTS_MAX),NSTATUS,NCTRLSUM SAVE/DJRES_DATA/ C. Local variables INTEGER NCALL,I1,I2 REAL*8 TEMP1,TEMP2 INTEGER NXLEFT,NYLEFT DATA NCALL/0/ C. For averaging REAL*8 VALLL,VALRL,VALLR,VALRR REAL*8 ER1LL,ER1RL,ER1LR,ER1RR,ER2LL,ER2RL,ER2LR,ER2RR REAL*8 WGAPX,WGAPY,WDXL,WDXR,WDYL,WDYR REAL*8 VALAL,VALAR,VALAA,ER1AL,ER1AR,ER1AA,ER2AL,ER2AR,ER2AA C. ----------------------------------------------- NCALL=NCALL+1 IF(NCALL.LT.NLOG)WRITE(LOGUNIT,*)'DJ_RES CALL N,X,Y:',NCALL,X,Y C. Set erroneous result by default DELTA=DJ_RES_UNKNOWN ERR_STA=0D0 ERR_SYS=DJ_RES_ERRUNKNOWN C. Initialize/Load data file. IF (NSTATUS.EQ.0) CALL DJ_RES_INIT() C. Check for internal error. IF (NSTATUS.GT.1) THEN IF (NCALL.LE.NLOG) WRITE(LOGUNIT,*)'DJ_RES ERROR' RETURN ENDIF C. Check parameters. IF((X.LT.XPTMIN).OR.(X.GT.XPTMAX).OR & .(Y.LT.YPTMIN).OR.(Y.GT.YPTMAX))THEN IF(NCALL.LE.NLOG)WRITE(LOGUNIT,*)'Invalid parameters:',X,Y RETURN ENDIF C. Find left X,Y-points in RPTSY CALL IDJR_FIND_LEFTPT(X,Y,NXLEFT,NYLEFT) C. The right points (NXLEFT+1) and (NYLEFT+1) must be valid. IF((NXLEFT.EQ.NPTSX).AND.(NXLEFT.GT.0))NXLEFT=NXLEFT-1 IF((NYLEFT.EQ.NPTSY).AND.(NYLEFT.GT.0))NYLEFT=NYLEFT-1 C. Check IF((NXLEFT.LT.1).OR.(NXLEFT.GT.NPTSX).OR & .(NYLEFT.LT.1).OR.(NYLEFT.GT.NPTSY)) THEN IF(NCALL.LE.NLOG)WRITE(LOGUNIT,*)'Error: cant find point:',X,Y RETURN ENDIF C. Calculate values for the four nearest points VALLL=RESDELTA(NXLEFT,NYLEFT) ER1LL=ERRSTA(NXLEFT,NYLEFT) ER2LL=ERRSYS(NXLEFT,NYLEFT) VALRL=RESDELTA(NXLEFT+1,NYLEFT) ER1RL=ERRSTA(NXLEFT+1,NYLEFT) ER2RL=ERRSYS(NXLEFT+1,NYLEFT) VALLR=RESDELTA(NXLEFT,NYLEFT+1) ER1LR=ERRSTA(NXLEFT,NYLEFT+1) ER2LR=ERRSYS(NXLEFT,NYLEFT+1) VALRR=RESDELTA(NXLEFT+1,NYLEFT+1) ER1RR=ERRSTA(NXLEFT+1,NYLEFT+1) ER2RR=ERRSYS(NXLEFT+1,NYLEFT+1) C. Check point for 'unknown' IF((VALLL.EQ.DJ_RES_UNKNOWN).AND.(VALRL.EQ.DJ_RES_UNKNOWN).AND & .(VALLR.EQ.DJ_RES_UNKNOWN).AND.(VALRR.EQ.DJ_RES_UNKNOWN))THEN IF(NCALL.LE.NLOG)WRITE(LOGUNIT,*)'DJ_RES: NO DATA FOR {X,Y}' RETURN ENDIF C. Fill single 'unknown' edge by neighbours + error IF(VALLL.EQ.DJ_RES_UNKNOWN)THEN IF(VALLR.NE.DJ_RES_UNKNOWN)THEN VALLL=VALLR ELSEIF(VALRL.NE.DJ_RES_UNKNOWN)THEN VALLL=VALRL ELSEIF(VALRR.NE.DJ_RES_UNKNOWN)THEN VALLL=VALRR ENDIF ER1LL=0D0 ER2LL=1D0 ENDIF IF(VALLR.EQ.DJ_RES_UNKNOWN)THEN IF(VALLL.NE.DJ_RES_UNKNOWN)THEN VALLR=VALLL ELSEIF(VALRR.NE.DJ_RES_UNKNOWN)THEN VALLR=VALRR ELSEIF(VALRL.NE.DJ_RES_UNKNOWN)THEN VALLR=VALRL ENDIF ER1LR=0D0 ER2LR=1D0 ENDIF IF(VALRL.EQ.DJ_RES_UNKNOWN)THEN IF(VALLL.NE.DJ_RES_UNKNOWN)THEN VALRL=VALLL ELSEIF(VALRR.NE.DJ_RES_UNKNOWN)THEN VALRL=VALRR ELSEIF(VALLR.NE.DJ_RES_UNKNOWN)THEN VALRL=VALLR ENDIF ER1RL=0D0 ER2RL=1D0 ENDIF IF(VALRR.EQ.DJ_RES_UNKNOWN)THEN IF(VALLR.NE.DJ_RES_UNKNOWN)THEN VALRR=VALLR ELSEIF(VALRL.NE.DJ_RES_UNKNOWN)THEN VALRR=VALRL ELSEIF(VALLL.NE.DJ_RES_UNKNOWN)THEN VALRR=VALLL ENDIF ER1RR=0D0 ER2RR=1D0 ENDIF C. Average over 4 edge points WGAPX=RPTSX(NXLEFT+1)-RPTSX(NXLEFT) WGAPY=RPTSY(NYLEFT+1)-RPTSY(NYLEFT) WDXL=X-RPTSX(NXLEFT) WDXR=RPTSX(NXLEFT+1)-X WDYL=Y-RPTSY(NYLEFT) WDYR=RPTSY(NYLEFT+1)-Y IF(WGAPX.GT.1D-14)THEN VALAL=(VALLL*WDXR+VALRL*WDXL)/WGAPX VALAR=(VALLR*WDXR+VALRR*WDXL)/WGAPX ER1AL=(ER1LL*WDXR+ER1RL*WDXL)/WGAPX ER1AR=(ER1LR*WDXR+ER1RR*WDXL)/WGAPX ER2AL=(ER2LL*WDXR+ER2RL*WDXL)/WGAPX ER2AR=(ER2LR*WDXR+ER2RR*WDXL)/WGAPX ELSE WRITE(LOGUNIT,*)'DJ_RES Warning. WGAPX too small:',WGAPX VALAL=(VALLL+VALRL)/2D0 VALAR=(VALLR+VALRR)/2D0 ER1AL=(ER1LL+ER1RL)/2D0 ER1AR=(ER1LR+ER1RR)/2D0 ER2AL=(ER2LL+ER2RL)/2D0 ER2AR=(ER2LR+ER2RR)/2D0 ENDIF IF(WGAPY.GT.1D-14)THEN VALAA=(VALAL*WDYR+VALAR*WDYL)/WGAPY ER1AA=(ER1AL*WDYR+ER1AR*WDYL)/WGAPY ER2AA=(ER2AL*WDYR+ER2AR*WDYL)/WGAPY ELSE WRITE(LOGUNIT,*)'DJ_RES Warning. WGAPY too small:',WGAPY VALAA=(VALAL+VALAR)/2 ER1AA=(ER1AL+ER1AR)/2 ER2AA=(ER2AL+ER2AR)/2 ENDIF DELTA=VALAA ERR_STA=ER1AA ERR_SYS=ER2AA IF(NCALL.LE.NLOG)WRITE(LOGUNIT,*)'DJ_RES=',DELTA,'+/-',ERR_STA & ,'+/-',ERR_SYS RETURN END C. ******************************************************************** C. IDJR_FIND_LEFTPT: Returns the nearest left point number C. Note: right point may not exist (if X==XPTMAX or Y==YPTMAX) C. ******************************************************************** SUBROUTINE IDJR_FIND_LEFTPT(X,Y,NXLEFT,NYLEFT) IMPLICIT NONE REAL*8 X,Y INTEGER NXLEFT,NYLEFT,I1,I2 C. block common/DJRES_DATA/ INTEGER NPTS_MAX,NPTS_MAX2,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,GSP,RESDELTA & ,ERRSTA,ERRSYS COMMON/DJRES_DATA/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,GSP & ,RESDELTA(NPTS_MAX,NPTS_MAX),ERRSTA(NPTS_MAX,NPTS_MAX) & ,ERRSYS(NPTS_MAX,NPTS_MAX),NSTATUS,NCTRLSUM C. Out-of-range checks IF((X.LT.XPTMIN).OR.(X.GT.XPTMAX).OR & .(Y.LT.YPTMIN).OR.(Y.GT.YPTMAX).OR & .(NPTSX.LE.0).OR.(NPTSY.LE.0)) THEN NXLEFT=-1 NYLEFT=-1 RETURN ENDIF C. Find left X-point in RPTSX() set -> NXLEFT I1=1 I2=NPTSX 120 NXLEFT=I1+(I2-I1)/2 IF(RPTSX(NXLEFT).LT.X)THEN I1=NXLEFT IF(I2.GT.(I1+1))GOTO 120 C. Additional check for X==XPTMAX case if(RPTSX(I2).EQ.X)NXLEFT=I2 ELSEIF(RPTSX(NXLEFT).GT.X)THEN I2=NXLEFT IF(I2.GT.(I1+1))GOTO 120 NXLEFT=I1 ENDIF C. Find left Y-point in RPTSY set -> NYLEFT I1=1 I2=NPTSY 121 NYLEFT=I1+(I2-I1)/2 IF(RPTSY(NYLEFT).LT.Y)THEN I1=NYLEFT IF(I2.GT.(I1+1))GOTO 121 C. Additional check for Y==YPTMAX case if(RPTSY(I2).EQ.Y)NYLEFT=I2 ELSEIF(RPTSY(NYLEFT).GT.Y)THEN I2=NYLEFT IF(I2.GT.(I1+1))GOTO 121 NYLEFT=I1 ENDIF RETURN END c. **************************************************************** c. DJ_RES_INIT: load data file c. c. File format: One number per line c. #(int) IVER: file version (reserved, currently 2) c. #(int) I2ERRS: 1(0) the stat. & sys. errors are separated c. #(int) S: S-Mp2-Ml2 c. #(int) NPTSX,NPTSY,NMODEX,NMODEY: number of points and grid mode c. #(float) XPTMIN,XPTMAX,YPTMIN,YPTMAX: min/max point values c. #(int) RESDELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY): results array c. #(int) NCTRLSUM: currently NPTSX*NPTSY SUBROUTINE DJ_RES_INIT() IMPLICIT NONE C. common/DJRES_SETS/ INTEGER NLOG,LOGUNIT REAL*8 DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN COMMON/DJRES_SETS/NLOG,LOGUNIT,DJ_RES_UNKNOWN,DJ_RES_ERRUNKNOWN C. common/DJRES_DATA/ INTEGER NPTS_MAX,NPTS_MAX2,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,GSP,RESDELTA & ,ERRSTA,ERRSYS COMMON/DJRES_DATA/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,GSP & ,RESDELTA(NPTS_MAX,NPTS_MAX),ERRSTA(NPTS_MAX,NPTS_MAX) & ,ERRSYS(NPTS_MAX,NPTS_MAX),NSTATUS,NCTRLSUM SAVE/DJRES_DATA/ C. 'unknown' value in dat-file (check for unknown as (VAL.GE.VAL_NO_DATA)) C. To be replaced with DJ_RES_UNKNOWN in memory. REAL*8 VAL_NO_DATA DATA VAL_NO_DATA/1D5/ C. local variables INTEGER I1,I2,ITEMP,IVER,NDATFILE,I2ERRS REAL*8 TEMP1,TEMP2,TEMP3 PARAMETER(NDATFILE=32) C. ------------------------------------------------ NSTATUS=2 I2ERRS=0 IVER=0 OPEN(NDATFILE,FILE='djangoh_r.dat',ERR=107,STATUS='OLD') C. TODO: Implement comment lines in dat-file READ(NDATFILE,*)IVER IF (IVER.NE.2) THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Unexpected data file format',IVER GOTO 106 ENDIF IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Data file loading' READ(NDATFILE,*)I2ERRS IF((I2ERRS.GT.1).OR.(I2ERRS.LT.0))THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Unexpected I2ERRS value: ',I2ERRS GOTO 106 ENDIF C. (real*8) GSP READ(NDATFILE,*)GSP IF(GSP.LE.0D0)THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Invalid GSP value: ',GSP GOTO 106 ENDIF C. (int) NPTSX,NPTSY,NMODEX,NMODEY READ(NDATFILE,*) NPTSX READ(NDATFILE,*) NPTSY READ(NDATFILE,*) NMODEX READ(NDATFILE,*) NMODEY IF(NLOG.GT.0)WRITE(LOGUNIT,*)'NPTSX,NPTSY:',NPTSX,NPTSY IF(NLOG.GT.0)WRITE(LOGUNIT,*)'NMODEX,NMODEY:',NMODEX,NMODEY IF((NPTSX.LE.0).OR.(NPTSX.GT.NPTS_MAX).OR & .(NPTSY.LT.0).OR.(NPTSY.GT.NPTS_MAX)) THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*) & 'Invalid data file: unexpected points number',NPTSX,NPTSY GOTO 106 ENDIF IF((NMODEX.LT.0).OR.(NMODEX.GT.1). & OR.(NMODEY.LT.0).OR.(NMODEY.GT.1)) THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*) & 'Invalid data file: unexpected grid mode',NMODEX,NMODEY GOTO 106 ENDIF C. (real*8) XPTMIN,XPTMAX,YPTMIN,YPTMAX READ(NDATFILE,*) XPTMIN READ(NDATFILE,*) XPTMAX READ(NDATFILE,*) YPTMIN READ(NDATFILE,*) YPTMAX IF(NLOG.GT.0)WRITE(LOGUNIT,*)'X min,max:',XPTMIN,XPTMAX IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Y min,max:',YPTMIN,YPTMAX IF((XPTMIN.GT.XPTMAX).OR.(YPTMIN.GT.YPTMAX)) THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*) & 'Invalid data file: unexpected grid margins', & XPTMIN,XPTMAX,YPTMIN,YPTMAX GOTO 106 ENDIF IF(((XPTMIN.LE.0D0).OR.(XPTMAX.LE.0D0)).AND.(NMODEX.EQ.1))THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*) & 'Invalid data file: NEGATIVE X MARGINS + LOG-SCALE', & XPTMIN,XPTMAX,NMODEX GOTO 106 ENDIF IF(((YPTMIN.LE.0D0).OR.(YPTMAX.LE.0D0)).AND.(NMODEY.EQ.1))THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*) & 'Invalid data file: NEGATIVE Y MARGINS + LOG-SCALE', & YPTMIN,YPTMAX,NMODEY GOTO 106 ENDIF DO 102 I1=1,NPTSX DO 102 I2=1,NPTSY READ(NDATFILE,*)TEMP1 READ(NDATFILE,*)TEMP2 IF(I2ERRS.EQ.1)THEN READ(NDATFILE,*)TEMP3 ELSE TEMP3=0D0 ENDIF IF(TEMP1.GE.VAL_NO_DATA)THEN TEMP1=DJ_RES_UNKNOWN TEMP2=0D0 TEMP3=DJ_RES_ERRUNKNOWN ENDIF RESDELTA(I1,I2)=TEMP1 ERRSTA(I1,I2)=TEMP2 ERRSYS(I1,I2)=TEMP3 102 CONTINUE READ(NDATFILE,*)NCTRLSUM IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Control sum:',NCTRLSUM IF (NCTRLSUM.NE.(NPTSX*NPTSY)) THEN IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Invalid control sum: ',NCTRLSUM GOTO 106 ENDIF IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Data file succesfully loaded' CLOSE(NDATFILE) GOTO 108 106 IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Error reading data file' CLOSE(NDATFILE) RETURN 107 IF(NLOG.GT.0)WRITE(LOGUNIT,*)'Error: no data file found' RETURN 108 CONTINUE C. Initialize points IF(NMODEX.EQ.1)THEN TEMP1=DLOG(XPTMIN) TEMP2=DLOG(XPTMAX) ENDIF DO 111 I1=1,NPTSX IF(NMODEX.EQ.0)THEN IF (I1.EQ.NPTSX) THEN RPTSX(I1)=XPTMAX ELSE RPTSX(I1)=XPTMIN+(XPTMAX-XPTMIN)*(I1-1)/(NPTSX-1) ENDIF ELSE IF (I1.EQ.1) THEN RPTSX(I1)=XPTMIN ELSEIF (I1.EQ.NPTSX) THEN RPTSX(I1)=XPTMAX ELSE RPTSX(I1)=DEXP(TEMP1+(TEMP2-TEMP1)*(I1-1)/(NPTSX-1)) ENDIF ENDIF 111 CONTINUE IF(NMODEY.EQ.1)THEN TEMP1=DLOG(YPTMIN) TEMP2=DLOG(YPTMAX) ENDIF DO 112 I1=1,NPTSY IF(NMODEY.EQ.0)THEN IF (I1.EQ.NPTSY) THEN RPTSY(I1)=YPTMAX ELSE RPTSY(I1)=YPTMIN+(YPTMAX-YPTMIN)*(I1-1)/(NPTSY-1) ENDIF ELSE IF (I1.EQ.1) THEN RPTSY(I1)=YPTMIN ELSEIF (I1.EQ.NPTSY) THEN RPTSY(I1)=YPTMAX ELSE RPTSY(I1)=DEXP(TEMP1+(TEMP2-TEMP1)*(I1-1)/(NPTSY-1)) ENDIF ENDIF 112 CONTINUE NSTATUS=1 IF(NLOG.GT.0)WRITE(LOGUNIT,*)'DJ_RES: Init OK' RETURN END FUNCTION IDJRES_INITOK() IMPLICIT NONE INTEGER IDJRES_INITOK INTEGER NPTS_MAX,NPTS_MAX2,NPTSX,NPTSY,NMODEX,NMODEY,NCTRLSUM & ,NSTATUS PARAMETER(NPTS_MAX=1000) PARAMETER(NPTS_MAX2=NPTS_MAX*NPTS_MAX) REAL*8 RPTSX,RPTSY,XPTMIN,XPTMAX,YPTMIN,YPTMAX,GSP,RESDELTA & ,ERRSTA,ERRSYS COMMON/DJRES_DATA/RPTSX(NPTS_MAX),RPTSY(NPTS_MAX) & ,XPTMIN,XPTMAX,YPTMIN,YPTMAX,NPTSX,NPTSY,NMODEX,NMODEY,GSP & ,RESDELTA(NPTS_MAX,NPTS_MAX),ERRSTA(NPTS_MAX,NPTS_MAX) & ,ERRSYS(NPTS_MAX,NPTS_MAX),NSTATUS,NCTRLSUM IDJRES_INITOK=0 IF(NSTATUS.EQ.1)IDJRES_INITOK=1 RETURN END