text:
	subroutine makefile(ch)
c
C	THIS IS MAKEFILE20
C
C	A PROGRAM TO SEARCH FOR LANGMUIR WAVES, MOSTLY UPSTREAM
C		 AND TO SAVE AMPLITUDE AND AMPLITUDE OF HARMONIC
C
C	ORIGINALLY, THE PROGRAM DID NOT COMPUTE ANGLE TO B.  ON 13 SEP 2007
C	THIS WAS ADDED AND A NEW SET OF RESULTS WRITTEN TO FOR080.DAT
C	note, 9 June 2008  I find only two criteria, bandwidth lt .02 freq
C		and freq between .75 fp and 1.25 fp
C
	COMMON /TDSF_BACKGROUND/ TDSFBCK(1025,2,4),NBCK(2,4)
C
	integer*4 ch,ok,okt,OK1,OK2,SCETI4(2),N1DATA(2048),N2DATA(2048)
	INTEGER*4 W_CHANNEL_CLOSE,W_EVENT,RET_SIZE,W_MESSAGES_OFF
	INTEGER*4 W_ITEM_I4,W_ITEM_R4,W_ITEM_R8
	INTEGER*4 TDS_CHANNEL,ISRC,SUNCLOCK,ERTDAY,ERT(2)
	REAL*8 SCET8,XGSE,YGSE,ZGSE,RGSE
	REAL DBSPEC1(1025),V1DATA(2048)
	REAL DBSPEC2(1025),V2DATA(2048)
	REAL ZCROSS(1024),Z1INT(1024),Z2INT(1024)
	REAL FREQ(1025),FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	REAL F2COUNT(2),F2AVR(2),F2STD(2),F2BW(2)
	REAL FPK(5),BWPK(5),AMPK(5)
	REAL VMATRIX(2,2),EVAL(2)
	character*32 ITEM
	character*4 event
	DATA TWOPI /6.2831853/
	DATA YXRATIO /11.8/
C	ABOVE BASED ON LEFFX = 41.1, LEFFY = 3.79
	DATA PXTOT,PYTOT /0.,0./
	DATA RE /6.378E3/
	DATA ANGTOB,TDIST /0.,0./
	DATA NFLLIMP,NFUPLIMP /40,1023/
	data event /'TDSF'/
C	data event /'FILL'/
C
	YXSQ = YXRATIO**2
	IFOUND = 0
C
C
C	CALCULATE FREQUENCIES FOR SPECTRUM
C
	SPS = 120000.
	DO N = 1,1025
	  FREQ(N) = (N-1)*SPS/2048.    ! FREQ(N) GOES WITH DBSPECT(N)
C	  DBSPECT(2) IS LOWEST NON-ZERO FREQUENCY, DBSPECT(1) IS DC
	ENDDO
C
C	OPEN(UNIT=90,FILE='MAKEFILE20.RESULTS',STATUS='OLD',ACCESS='APPEND')
 100	ok = w_event(ch,event)
C
C	CHECK FOR END OF RECORD
C
	if (ok.ne.1) then
		if(ok.eq.82) then
		   ok = w_channel_close(ch)
		   CLOSE(UNIT=90)
		   return
	        endif
		write(6,*) 'cannot open ',event, ', ok=', ok
	endif
C
C	END OF END OF RECORD CHECK
C
	OKT = W_MESSAGES_OFF(ch)
C
C
	ITEM = 'CHANNEL'
	ok = W_ITEM_I4(ch, item, TDS_CHANNEL, 1, ret_size)
C
	IPROS = 4
	IF(TDS_CHANNEL.EQ.1) THEN
C	  ITEM = 'DATA'
C	  ok1 = W_ITEM_I4(ch, item, N1DATA, 2048, ret_size)
	  CALL TDS_PHYS(CH,IPROS,N1DATA,V1DATA,DBSPEC1)
	  ITEM = 'EVENT_NUMBER'
	  ok = W_ITEM_I4(ch, item, NO_EVT1, 1, ret_size)
		MAXCH1 = 0
		DO N = 1,2048
		  MAXCH1 = MAX0(MAXCH1,IABS(N1DATA(N)-128))
		ENDDO
	  GO TO 100
	ELSEIF(TDS_CHANNEL.EQ.2) THEN
C	  ITEM = 'DATA'
C	  ok2 = W_ITEM_I4(ch, item, N2DATA, 2048, ret_size)
	  CALL TDS_PHYS(CH,IPROS,N2DATA,V2DATA,DBSPEC2)
	  ITEM = 'EVENT_NUMBER'
	  ok = W_ITEM_I4(ch, item, NO_EVT, 1, ret_size)
	  IF(NO_EVT.NE.NO_EVT1) GO TO 100
		MAXCH2 = 0
		DO N = 1,2048
		  MAXCH2 = MAX0(MAXCH2,IABS(N2DATA(N)-128))
		ENDDO
	ELSE
	  GO TO 100
	ENDIF

C
CDIAG	PRINT*,'******START EVENT',NO_EVT1,NO_EVT,' *****'
	PRINT*,'******START EVENT',NO_EVT1,NO_EVT,' *****'
C
	IF(TDS_CHANNEL.EQ.1) GO TO 100		! GO GET CHANNEL 2
	IF(TDS_CHANNEL.NE.2) GO TO 100
C
C	CALCULATE VMATRIX, FOR ANGLE TO B
C
        VMATRIX(1,1) = 0.
        VMATRIX(1,2) = 0.
        VMATRIX(2,2) = 0.
	XEFFLEN = 41.1
	YEFFLEN = 3.79
        DO N = 1,2048
          VMATRIX(1,1) = VMATRIX(1,1) + (V1DATA(N)/XEFFLEN)**2
          VMATRIX(1,2) = VMATRIX(1,2) + V1DATA(N)*V2DATA(N)/XEFFLEN/
     1	YEFFLEN
          VMATRIX(2,2) = VMATRIX(2,2) + (V2DATA(N)/YEFFLEN)**2
        ENDDO
C
	XRMS = 1000.*SQRT(VMATRIX(1,1)/2048.)    ! in mV/m
	YRMS = 1000.*SQRT(VMATRIX(2,2)/2048.)
C
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, XGSE, 1, ret_size)
	XRE = XGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, YGSE, 1, ret_size)
	YRE = YGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, ZGSE, 1, ret_size)
	ZRE = ZGSE/RE
	RRE = SQRT(XRE**2+YRE**2+ZRE**2)
C
C	IF(RGSE.LT.6.4D05) THEN			! LIMIT 100 RE
C	    ADVANCE BY ONE HOUR AND GET ANOTHER EVENT
C	    ITEM = 'EVENT_SCET_R8'
C	    OK = W_ITEM_R8(ch, item, SCET8, 1, ret_size)
	    ITEM = 'EVENT_SCET'
	    ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
C	    PRINT*,RGSE/6.38E3,'INSIDE 100 RE',SCETI4
C	    SCET8 = SCET8 + 1.D00/24.D00
C	    CALL W_CHANNEL_POSITION(CH,SCET8)
C	    GO TO 100
C	  ENDIF
C
	ITEM = 'SOURCE'
	ok = W_ITEM_I4(ch, item, ISRC, 1, ret_size)
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(ch, item, NO_EVT, 1, ret_size)
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
	ITEM = 'FAST_RX_SPEED'
	ok = W_ITEM_I4(ch, item, ISPS, 1, ret_size)
	IF(ISPS.NE.0) GO TO 100
	ITEM = 'WIND_MFI_BX(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BX, 1, ret_size)
	ITEM = 'WIND_MFI_BY(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BY, 1, ret_size)
	ITEM = 'WIND_MFI_BZ(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BZ, 1, ret_size)
	ITEM = 'WIND_MFI_BPHI(GSE)_R4'
	OK = W_ITEM_R4(CH, ITEM, AZMAG,1,RET_SIZE)
	ITEM = 'WIND_MFI_BTHETA(GSE)_R4'
	OK = W_ITEM_R4(CH, ITEM, THETAB,1,RET_SIZE)
	ITEM = 'SUN_ANGLE'
	ok = W_ITEM_I4(ch, item, SUNCLOCK, 1, ret_size)
C	ITEM = 'WIND_3DP_E_TEMP_R4'
C	ok = W_ITEM_R4(ch, item, TEMPE, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_TEMP_R4'
C	ok = W_ITEM_R4(ch, item, TEMPI, 1, ret_size)
	ITEM = 'WIND_3DP_ION_DENSITY_R4'
	ok = W_ITEM_R4(ch, item, DENS, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_VX(GSE)_R4'
C	ok = W_ITEM_R4(ch, item, VX, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_VY(GSE)_R4'
C	ok = W_ITEM_R4(ch, item, VY, 1, ret_size)
C	ITEM = 'WIND_MFI_BMAG_R4'
C	ok = W_ITEM_R4(ch, item, BMAG, 1, ret_size)
	ITEM = 'WIND_SPIN_RATE_R4'
	ok = W_ITEM_R4(ch, item, SPINRATE, 1, ret_size)
	if(ok.ne.1) then
	  spinrate = 2.
	endif
	ITEM = 'EVENT_TM_SCET_I4'
	ok = W_ITEM_I4(ch, item, ERT, 2, ret_size)
	ERTDAY = MOD(ERT(1),100)
C
	EANGLE =  -360.*(SUNCLOCK-14.)/4096. - 45. ! ANGLE SUN TO +EX AT END
	IF(EANGLE.LT.-180.) EANGLE = EANGLE + 360.
	IF(EANGLE.GT.180.)  EANGLE = EANGLE - 360.
	DANG = SPINRATE*360./SPS/TWOPI
	THANT = EANGLE + 3072.*DANG  ! ANGLE SUN TO +EX AT START 16nov99
	THANTM = THANT + 1024.*DANG  ! ANGLE SUN TO +EX AT MIDDLE 16nov99
C
C	EVALUATE MINIMUM VARIANCE STUFF USING FACT THAT TRACE
C		AND DETERMINANT ARE UNCHANGED BY ROTATIONS.
C		THIS IS ENOUGH FOR A 2 X 2 MATRIX
C
	DET = VMATRIX(1,1)*VMATRIX(2,2) - VMATRIX(1,2)**2
	TRACE = VMATRIX(1,1) + VMATRIX(2,2)
	DISC = SQRT(TRACE**2 - 4.*DET)
	EVAL(1) = .5*(TRACE + DISC)
	EVAL(2) = .5*(TRACE - DISC)
	AXRATIO = SQRT(EVAL(2)/EVAL(1))
C	AX MAJOR AND AX MINOR ARE RATIOS OF Y COMPONENT TO X COMPONENT OF
C		SLOPES OF MAJOR AND MINOR EIGENVECTORS
C		AND GIVE DIRECTION OF E.
c
	AXMAJOR = (EVAL(1)-VMATRIX(1,1))/VMATRIX(1,2)
	AXMINOR = (EVAL(2)-VMATRIX(1,1))/VMATRIX(1,2)
C	SLOPE = YXRATIO*AXMAJOR
	SLOPE = AXMAJOR
	print*,'vmatrix',vmatrix(1,1),vmatrix(1,2),vmatrix(2,2)
	GSEANG = THANTM - ATAND(SLOPE)
	ANGTOB = GSEANG - AZMAG
c	print*,'eval',eval
c	print*,'thant',thant
c	print*,'azmag,bx,by',azmag,bx,by
c	print*,'axmajor,axminor,slope',axmajor,axminor,slope
c	print*,'angtob',angtob
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
c        print*,'angtob',angtob
c        if(no_evt.eq.12115717) stop
C
C	ITEM = 'CAL_PA'
C	ok = W_ITEM_I4(ch, item, ICAL, 1, ret_size)
C	IF(ICAL.NE.0) GO TO 100
C	HAVE TO EXCLUDE CALS BY TIMING.
c		to exclude cals
	  if(sceti4(2).ge.010000.and.sceti4(2).lt.011000) THEN
		PRINT*,'CALIBRATE'
		GO TO 100
	  endif
C
	CALL DIFFFK(XRE,YRE,ZRE,BX,BY,BZ,DDIFF,XT,YT,ZT,XI,YI,ZI,IN)
C
	TDIST = (XRE-XT)**2 + (YRE-YT)**2 + (ZRE-ZT)**2
	TDIST = SQRT(TDIST)
C
C
C	  FIND MAXIMUM IN SPECTRUM
C
	  SMAX = -1000.
	  DO N = 35,1022			! 35 IS ABOUT 2 kHz
	    IF(DBSPEC1(N).GT.SMAX) THEN
	      SMAX = DBSPEC1(N)
	      NSPMAX = N
	    ENDIF
	  ENDDO
	PRINT*,'X SPECTRUM MAX:N,F,SMAX',NSPMAX,(NSPMAX-1)*SPS/2048.,SMAX
	  FREQMAX1 = (NSPMAX-1)*SPS/2048.
C
	  SMAX = -1000.
	  DO N = 35,1022			! 35 IS ABOUT 2 kHz
	    IF(DBSPEC2(N).GT.SMAX) THEN
	      SMAX = DBSPEC2(N)
	      NSPMAX = N
	    ENDIF
	  ENDDO
	PRINT*,'Y SPECTRUM MAX:N,F,SMAX',NSPMAX,(NSPMAX-1)*SPS/2048.,SMAX
	  FREQMAX2 = (NSPMAX-1)*SPS/2048.
C
C	  AVERAGE OVER ADJACENT POINTS, THIS IS TEMPORARY FOR COMPARISON
	  FTOT = 0.
	  STOT = 0.
	  DO NN = -1,1
	    FNN = (NSPMAX-1+NN)*SPS/2048.
	    SNN = 10.**(.1*DBSPEC1(NSPMAX+NN))
	    FTOT = FTOT+FNN*SNN
	    STOT = STOT+SNN
	  ENDDO
	  FREQMAX = FTOT/STOT
	  PRINT*,'INTERPOLATED FREQ',FREQMAX
C	  WRITE(76,*) 'INTERPOLATED FREQ',FREQMAX
C
C	CRITERIA
C
	print*,'start criteria, diff=',ddiff
c	IF(DDIFF.GT.0.) GO TO 100
C	if(no_evt.ne.74375) go to 100
	call FANDBW(CH,FREQ,DBSPEC1,DBSPEC2,NPEAK,FAVR,FBW,SKEW)

C
	XMAX = 0.
	YMAX = 0.
	XMIN = 0.
	YMIN = 0.
	DO N = 1,2048
	  XMAX = AMAX1(XMAX,V1DATA(N))
	  YMAX = AMAX1(YMAX,V2DATA(N))
	  XMIN = AMIN1(XMIN,V1DATA(N))
	  YMIN = AMIN1(YMIN,V2DATA(N))
	ENDDO
C	print*,'xmin,max,y',xmin,xmax,ymin,ymax
	XMAX = AMAX1(XMAX,-XMIN)
	YMAX = AMAX1(YMAX,-YMIN)
C	CHANGE TO mV/m
C	 BASED ON LEFFX = 41.1, LEFFY = 3.79
	XMAX = 1000.*XMAX/41.1
	YMAX = 1000.*YMAX/3.79
	  IFOUND = 0
C	  GO TO 100
C
	FP = 0.
	IF (DENS.GT.0.) FP = AMAX1(8.E3,9.E3*SQRT(DENS))     ! PLASMA FREQ IN KHZ
	IF(FREQMAX.LT..75*FP.OR.FREQMAX.GT.1.25*FP) THEN
	  IF(FP.NE.0.) THEN
     	   PRINT*,'FREQ NOT CLOSE TO PLASMA FREQ, FP,RATIO-',FP,FREQMAX/FP
	  ELSE
     	   PRINT*,'FREQ NOT CLOSE TO PLASMA FREQ, FP,F ',FP,FREQMAX
	  ENDIF
	  GO TO 100
	ENDIF
C
C	ADD TO BACKGROUND
C
	IBCKG = 0
	IF(IBCKG.NE.0) THEN
	  XMXLIM = .1
	  IF(XMAX.LT.XMXLIM.AND.TDS_CHANNEL.EQ.1) THEN
	    DO N = 1,1025
	      TDSFBCK(N,TDS_CHANNEL,ISPS+1) =
     1      TDSFBCK(N,TDS_CHANNEL,ISPS+1) + DBSPEC1(N)
	    ENDDO
	    NBCK(TDS_CHANNEL,ISPS+1) = NBCK(TDS_CHANNEL,ISPS+1) +1
	  ENDIF
C
	  XMXLIM = .2
	  IF(YMAX.LT.XMXLIM.AND.TDS_CHANNEL.EQ.2) THEN
	    DO N = 1,1025
	      TDSFBCK(N,TDS_CHANNEL,ISPS+1) =
     1      TDSFBCK(N,TDS_CHANNEL,ISPS+1) + DBSPEC1(N)
	    ENDDO
	    NBCK(TDS_CHANNEL,ISPS+1) = NBCK(TDS_CHANNEL,ISPS+1) +1
	  ENDIF
	ENDIF		! END OF ADD TO BACKGROUND
C
C	CALCULATE BANDWIDTH, ETC.
	NFLLIM = NSPMAX/2
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFUPLIM = (3*NSPMAX)/2
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	PRINT*,'DO FREQS',FREQ(NFLLIM),FREQ(NFUPLIM)
	CONST1 = (1000./41.1)**2
	CONST2 = (1000./3.79)**2
	DO NC = 1,2
	  FCOUNT(NC) = 1.E-20
	  FAVR(NC) = 0.
	  FSTD(NC) = 0.
	  F3MOM(NC) = 0.
	  DO N = NFLLIM,NFUPLIM
	    IF(NC.EQ.1) THEN
		VLT = CONST1*10.**(.1*DBSPEC1(N))
	    ELSE
		VLT = CONST2*10.**(.1*DBSPEC2(N))
	    ENDIF
	    FCOUNT(NC) = FCOUNT(NC) + VLT
	    FAVR(NC) = FAVR(NC) + VLT*FREQ(N)
	    FSTD(NC) = FSTD(NC) + VLT*FREQ(N)**2
	    F3MOM(NC) = F3MOM(NC) + VLT*FREQ(N)**3
	  ENDDO
	  FAVR(NC) = FAVR(NC)/FCOUNT(NC)
	  FSTD(NC) = FSTD(NC)/FCOUNT(NC) - FAVR(NC)**2
	  FBW(NC) = SQRT(AMAX1(FSTD(NC),0.))
	  SKEW(NC) = F3MOM(NC)/FCOUNT(NC) - 3.*FAVR(NC)*FSTD(NC)
     1		 - FAVR(NC)**3
	print*,' bw and .02*f',fbw(nc),.02*favr(nc)
	  IF(FBW(NC).GT..02*FAVR(NC)) GO TO 100
	ENDDO
C
C	DO SAME FOR HARMONIC
C
C	CALCULATE BANDWIDTH, ETC.
	NBW = 3
	NFLLIM = 2*NSPMAX - 1 - NBW
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFLLIM = MIN0(NFLLIM,NFUPLIMP)
	NFUPLIM = 2*NSPMAX - 1 + NBW
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	PRINT*,'CHECK NFUPLIM, ETC.',NFLLIM,NFUPLIM
	IF(NFUPLIM.GT.0)
     1		PRINT*,'DO HARM FREQS',FREQ(NFLLIM),FREQ(NFUPLIM)
	DO NC = 1,2
	  F2COUNT(NC) = 1.E-20
	  F2AVR(NC) = 0.
	  F2STD(NC) = 0.
	  DO N = NFLLIM,NFUPLIM
	    IF(NC.EQ.1) THEN
		PWR = CONST1*10.**(.1*DBSPEC1(N))
	    ELSE
		PWR = CONST2*10.**(.1*DBSPEC2(N))
	    ENDIF
	    F2COUNT(NC) = F2COUNT(NC) + PWR
	    F2AVR(NC) = F2AVR(NC) + PWR*FREQ(N)
	    F2STD(NC) = F2STD(NC) + PWR*FREQ(N)**2
	  ENDDO
	  F2AVR(NC) = F2AVR(NC)/F2COUNT(NC)
	  F2STD(NC) = F2STD(NC)/F2COUNT(NC) - F2AVR(NC)**2
	  F2BW(NC) = SQRT(AMAX1(F2STD(NC),0.))
	ENDDO
C
c*************
c	CALL FANDBW(SPS,TDS_CHANNEL,DBSPEC1,DBSPEC2,NPEAK,FPK,BWPK,AMPK)
c	print*,'fandbw',npeak,fpk(1),bwpk(1),ampk(1)
c	if(1) stop
c*************
C	FCE = .028*BMAG
C	BETA = .402*(TEMPE+TEMPI)/BMAG**2
C
c	IF(FAVR(1).GT..8*FP) THEN
	  PRINT*,'EVENT NO.',NO_EVT,' FREQ', FAVR(1),' MAX',
     1		MAXCH1,MAXCH2
c	  IFOUND = 0
c	ENDIF
C
C*****	CHECKING
c	PRINT*,'VMATRIX',VMATRIX(1,1),VMATRIX(1,2),VMATRIX(2,2)
c	PRINT*,'EVALS,AXS,PX,PY',EVAL(1),EVAL(2),AXMAJOR,AXMINOR,
c     1		PX,PY
C*****
C
C	PLASMA PARAMETERS, ETC.

C	   write(s,'(i8.8,i6.6)',iostat=ios) s_scet(1), s_scet(2)
C	   s_scet = s(1:4)//'/'//s(5:6)//'/'//s(7:8)//' '//
C	1	s(9:10)//':'//s(11:12)//':'//s(13:14)
C
c	print date, time, event no., freq, Px, Py, Vy/Vx, angles of
c		antenna,B, phiB, mag B or fce, fp
c
c	PRINT*, 'PXTOT,PYTOT,RATIO',PXTOT,PYTOT,XYRATIO
C
	print*,'write 90'
C	WRITE FREQ,AMP,BW FOR CH 1 AND 2 FUNDAMENTAL, AMP AND BW FOR
C	    HARMONIC, DIFF, DISTANCE TO TANGENT PT., FP, ANGLE TO B,
	WRITE(90,1011) sceti4,NO_EVT,EVENT(1:1),.001*FAVR(1),XMAX,
     1    .001*FBW(1),.001*FAVR(2),YMAX,.001*FBW(2),.001*FP,
     2	  ANGTOB,DDIFF,TDIST,FCOUNT(1),FCOUNT(2),F2COUNT(1),F2COUNT(2)
C
 1011	format(I10,I7,I10,A2,F6.2,E9.2,F6.3,F6.2,E9.2,F6.3,F6.2,F7.1,
     1	F6.1,F6.1,4E9.2)
C
        WRITE(80,1080) sceti4,NO_EVT,ERTDAY,EVENT(1:1),.001*FAVR(1),
     1    XRMS,.001*FBW(1),.001*FAVR(2),YRMS,.001*FBW(2),.001*FP,
     2    ANGTOB,AXRATIO,THETAB,DDIFF,XRE,YRE,ZRE
C
 1080   format(I10,I7,I10,I3,A2,F6.2,E9.2,F6.3,F6.2,E9.2,F6.3,F6.2,F7.1,
     1  F6.2,F6.1,4F7.1)
C
	IF(FAVR(1).LT.25000.) THEN
	  TPWR = FCOUNT(1) + FCOUNT(2)
	  T2PWR = F2COUNT(1) + F2COUNT(2)
	  WRITE(91,1012) sceti4,NO_EVT,EVENT(1:1),.001*FAVR(1),
     1	   .001*F2AVR(1),F2AVR(1)/FAVR(1),FCOUNT(1),F2COUNT(1),
     2	   TPWR,T2PWR
 1012	format(I10,I7,I10,A2,2F6.2,F6.3,4E11.3)
	ENDIF
C
	IF(OK.NE.82) GO TO 100
	CLOSE(UNIT=90)
	return
	end
	SUBROUTINE FANDBW(CH,FREQ,DBX,DBY,NPEAK,FAVR,FBW,SKEW)
C
	COMMON /NOISE_LIM/ XNS
	INTEGER CH
	DATA SPS /120000./
	REAL AMP(2),DBX(1025),DBY(1025),DBSUM(1025),FREQ(1025)
	REAL FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	REAL X(25)
	DATA XNS /-120./
	DATA NFLLIMP,NFUPLIMP /4, 1024/
	DATA SPS /120000./
C
C	ROUTINE TO FIT FUNCTION TO LOG OF SPECTRAL POWER
C
C	CONST1 = (1./41.1)**2		! DO IN (VOLTS/M)**2
C	CONST2 = (1./3.79)**2
	LEFFX = 41.1
	LEFFY = 3.79
C
C	CALCULATE MAGNITUDE OF FIELD
C
	DO N = 1,1025
	  TPWR = 10.**(.1*DBX(N))/LEFFX**2 +
     1		10.**(.1*DBY(N))/LEFFY**2
	  DBSUM(N) = 10.*ALOG10(TPWR)
	ENDDO
C
C	FIND FREQUENCY OF HIGHEST PEAK BETWEEN N1 AND N2
C
	N1 = 4			! LOWEST FREQ 177 HZ
	N2 = 1024
	INC = 1
C
	SMAX = DBSUM(N1)
	SSAVE = DBSUM(N1)
	NPEAK = N1
	SLIM = AMAX1(XNS,SSAVE-80.)
	DO N = N1,N2
	    IF(DBSUM(N).GT.SMAX) THEN
	        SMAX = DBSUM(N)
	        NPEAK = N
	    ENDIF
	ENDDO
cdiag	PRINT*,'in FANDBW, MAX: N,F,SMAX',NPEAK,FREQ(NPEAK),SMAX
	FPK = FREQ(NPEAK)
	print*,'freq check',(npeak-1)*sps/2048.,fpk
C
C	LIMIT TO FAIRLY LARGE PEAKS
C
c	IF(DBSUM(NPEAK).LT.(XNS+20.)) GO TO 100
c	IF(NPK.LE.10) GO TO 100
C
C	DETERMINE DATA TO USE
C
	N2 = NPEAK+1
	SLIM = AMAX1(XNS,SMAX-25.)
	DO N = NPEAK+1,1024
	  IF(DBSUM(N).LT.SLIM) GO TO 22
	  N2 = N
	ENDDO
C
 22	N1 = NPEAK-1
	DO N = NPEAK-1,3,-1
	  IF(DBSUM(N).LT.SLIM) GO TO 24
	  N1 = N
	ENDDO
 24	CONTINUE
	print*,'fit range',n1,n2
C	if(1) stop
C
C	FIT
C
	RETURN
	END
	SUBROUTINE FANDBW2(N1,N2,FREQ,DBX,DBY,NPEAK,FAVR,FBW,SKEW)
C
	DATA SPS /120000./
	REAL AMP(2),DBX(1025),DBY(1025),DBSUM(1025),FREQ(1025)
	REAL FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	DATA NFLLIMP,NFUPLIMP /4, 1024/
C
	CONST1 = (1./41.1)**2		! DO IN (VOLTS/M)**2
	CONST2 = (1./3.79)**2
C
C	CALCULATE BANDWIDTH, ETC.
C
	NFLLIM = NPEAK/2
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFUPLIM = (3*NPEAK)/2
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	IF(N1.GT.N2) THEN
	   N2T = N1
	   N1T = N2
	ELSE
	   N1T = N1
	   N2T = N2
	ENDIF
C*********
C	THIS CHOICE DID NOT WORK WELL.  IT GAVE NARROW BANDWIDTHS
C	FOR SOME CASES WHEN THE BANDWIDTH WAS WIDE.  FOR EXAMPLE, IT
C	GAVE A BANDWIDTH OF 63.7 HZ FOR 19960421 # 13042040
C*********
C
	NPEAK = N1T
	VMAX = 0.
	DO NC = 1,2
	  FCOUNT(NC) = 1.E-20
	  FAVR(NC) = 0.
	  FSTD(NC) = 0.
	  F3MOM(NC) = 0.
	  DO N = N1T,N2T
	    IF(NC.EQ.1) THEN
		VLT = CONST1*10.**(.1*DBX(N))
	    ELSE
		VLT = CONST2*10.**(.1*DBY(N))
	    ENDIF
	    IF(VLT.GT.VMAX) THEN
	       NPEAK = N
	       VMAX = VLT
	    ENDIF
	    FCOUNT(NC) = FCOUNT(NC) + VLT
	    FAVR(NC) = FAVR(NC) + VLT*FREQ(N)
	    FSTD(NC) = FSTD(NC) + VLT*FREQ(N)**2
	    F3MOM(NC) = F3MOM(NC) + VLT*FREQ(N)**3
	  ENDDO
	  FN = FCOUNT(NC)
c	print*,'sk',fn,f3mom(nc)/fn,3.*fstd(nc)/fn**2,2.*(favr(nc)/fn)**2
c	  SKEW(NC) = F3MOM(NC)/FN - 3.*FSTD(NC)/FN**2
c     1		 + 2.*(FAVR(NC)/FN)**3
	  FAVR(NC) = FAVR(NC)/FCOUNT(NC)
	  FSTD(NC) = FSTD(NC)/FCOUNT(NC) - FAVR(NC)**2
	  FBW(NC) = SQRT(AMAX1(FSTD(NC),0.))
	  SKEW(NC) = F3MOM(NC)/FCOUNT(NC) - 3.*FAVR(NC)*FSTD(NC)
     1		 - FAVR(NC)**3
	ENDDO
C
	RETURN
	END


text:
      INCLUDE 'VICMAIN_FOR'

C VICAR program GRIDGEN
C Generates a perfect grid: either a MARK format grid location file, or an
C image.
C
      SUBROUTINE MAIN44

!     Program GRIDGEN local variables
      IMPLICIT INTEGER(A-Z)
      COMMON/C1/OBUF(1200)

      INTEGER*4 GSL,GSS,DNG,DNB
      LOGICAL XVPTST,IMAGE

      INTEGER INC, CNT, NCOL, NROW, DUM, OUNIT, STAT

!     Initialize GRIDGEN local variables
      DATA INC /0/, CNT /0/, GSL /0/, GSS /0/, NROW /0/, DNG /0/
     &     NCOL /0/, DNB /0/, DUM /0/, OUNIT /0/, STAT /0/
      DATA OBUF /1200*0/

!     Begin program GRIDGEN

!     14apr2012 -lwk- changed IFMESSAGE to XVMESSAGE to avoid garbage output
!     on Linux
      CALL XVMESSAGE ('GRIDGEN version 1 July 1994', ' ')

!     Set default error handling action
      CALL XVEACTION ('SA', ' ')

!     Determine if the keyword 'IMAGE' was set on input parameters
      IMAGE = XVPTST('IMAGE')

!     Obtain information about values for the indicated parameters
      CALL XVP('INC',INC,CNT)
      CALL XVP('GSL',GSL,CNT)
      CALL XVP('GSS',GSS,CNT)
      CALL XVP('NROW',NROW,CNT)
      CALL XVP('NCOL',NCOL,CNT)
      CALL XVP('DNGRID',DNG,CNT)
      CALL XVP('DNBACK',DNB,CNT)

!     If the keyword 'IMAGE' was set in the input parameters
      IF (IMAGE) GOTO 100
      NSO = NROW*NCOL*2
      IF (NSO.GT.1200) THEN
           CALL XVMESSAGE ('***Too many grid intersections', ' ')
           CALL XVMESSAGE ('***Reduce NROW or NCOL', ' ')
           GOTO 999
      ENDIF

      CALL GENGRID (OBUF,NROW,NCOL,GSL,GSS,INC)
      CALL PGRID   (OBUF,NROW,NCOL,DUM,0)
      CALL XVUNIT  (OUNIT,'OUT',1,STAT, ' ')
      CALL XVOPEN(OUNIT,STAT,'U_NL',1,'U_NS',NSO,
     *     'O_FORMAT','REAL','OP','WRITE', ' ')
      CALL XVWRIT(OUNIT,OBUF,STAT, ' ')
      GOTO 200

  100 NLO = GSL + NROW*INC - INC/2
      NSO = GSS + NCOL*INC - INC/2
      CALL XVUNIT (OUNIT,'OUT',1,STAT, ' ')
      CALL XVOPEN (OUNIT,STAT,'U_NL',NLO,'U_NS',NSO,
     *            'O_FORMAT','BYTE','OP','WRITE', ' ')

      CALL GENIMAGE (OUNIT,DNB,DNG,NLO,NSO,NROW,NCOL,GSL,GSS,INC)

  200 CALL XVMESSAGE ('GRIDGEN task completed', ' ')
      RETURN
  999 CALL XVMESSAGE ('***GRIDGEN task canceled', ' ')
      CALL ABEND
      END
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!                 Subroutine to generate grid locations...
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
C
      SUBROUTINE GENGRID (OBUF,NROW,NCOL,GSL,GSS,INC)

      IMPLICIT INTEGER(A-Z)

!     Subroutine GENGRID passed parameteres
      REAL*4 OBUF(1200)
      INTEGER*4 GSL,  GSS
      INTEGER   NROW, NCOL, INC

!     Subroutine GENGRID local parameteres
      INTEGER*4 GSS0
      INTEGER   I, J, K

!     Subroutine GENGRID local parameter initialization
      DATA      I /0/, J /0/, K /0/
      DATA      GSS0 /0/

!     Begin Subroutine GENGRID
      GSS0 = GSS
      K = 1

      DO J=1,NROW	!ROWS
          DO I=1,NCOL   !COLUMNS
              OBUF(K) = GSL
              OBUF(K+1) = GSS
              K = K + 2
              GSS = GSS + INC
          ENDDO
          GSS = GSS0
          GSL = GSL + INC
      ENDDO
      RETURN
      END
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!                 Subroutine GENIMAGE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
C
      SUBROUTINE GENIMAGE(OUNIT,DNB,DNG,NLO,NSO,NROW,NCOL,GSL,GSS,INC)

      IMPLICIT INTEGER(A-Z)

!     Subroutine GENIMAGE passed parameteres
      INTEGER*4 GSL,GSS,DNG,DNB
      INTEGER   OUNIT, NLO, NSO, NROW, NCOL, INC

!     Subroutine GENIMAGE local parameteres
      COMMON/C2/PIC(1200),PIC2(1200)
      BYTE PIC,PIC2
      INTEGER   L, STAT

!     Initialize Subroutine GENIMAGE local parameters
      DATA      STAT /0/, L /0/
      DATA      PIC /1200*0/, PIC2 /1200*0/

!     Begin Subroutine GENIMAGE
      CALL ITLA(DNB,PIC,NSO)			! Start with background
      CALL MVE(-5,NCOL,DNG,PIC(GSS),0,INC)	! Add vertical grid rulings
      CALL ITLA(DNG,PIC2,NSO)			! Generate horizontal grid

      DO L=1,NLO
          IF (L.EQ.GSL) THEN
              CALL XVWRIT(OUNIT,PIC2,STAT, ' ')
              GSL = GSL + INC
          ELSE
              CALL XVWRIT(OUNIT,PIC,STAT, ' ')
          ENDIF
      ENDDO

      RETURN
      END

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
  Coalition
  Consortium
  Concord
  Collaboration
  American Brontosaurus Cabal of Davis Eristopia Cabal
  Conspiracy
  Cooperative
  Competition
  Committee
  Commission
  Chairmanship
  Connection
  Corporation
  Company
  Crusade
  Conference
  Conjunction
  Caucus
  Cartel
  Cabinet
  Collective
  Coven
  Covenant

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


program slice
   type coarr
      real, allocatable :: a(:, :, :)[:, :, :]
   end type

   type(coarr) :: co
   integer :: nimg, me, z, nx, ny, nz, north, south, mex, mey, mez, coords(3)
   integer :: shape2d(2), shape3d(3)
   real, allocatable :: buf3d(:, :, :) ! a plane slice as a rank 3 array with a single transverse layer
   real, allocatable :: buf2d(:, :) ! a plane (2d) slice, normal in the y direction

   nx = 6
   ny = 4
   nz = 2

   me = this_image()
   nimg = num_images()

   if (nimg /= 8) stop

   allocate(co % a(nx, ny, nz)[1:2, 1:2, *])
   allocate(buf2d(nx, nz), buf3d(nx, 1, nz))

   ! this example should NOT reallocate buf2d nor buf3d
   ! compare shapes before and after syncing
   shape2d = shape(buf2d)
   shape3d = shape(buf3d)

   co % a = reshape([(z, z=1, nx * ny * nz)], shape(co % a))

   coords = this_image(co % a)
   mex = coords(1)
   mey = coords(2)
   mez = coords(3)

   north = mey + 1
   south = mey - 1

   sync all
   if (north <= 2) then
      z = image_index(co % a, [mex, north, mez])
      sync images(z)
      ! no reduction on rank
      buf3d = co % a(1:nx, 1:1, 1:nz)[mex, north, mez]
      co % a(1:nx, ny:ny, 1:nz) = buf3d

      ! reduction along dim 2
      buf2d = co % a(1:nx, 1, 1:nz)[mex, north, mez]
      co % a(1:nx, ny, 1:nz) = buf2d
   end if
   if (south >= 1) then
      z = image_index(co % a, [mex, south, mez])
      sync images(z)
      buf3d = co % a(1:nx, ny:ny, 1:nz)[mex, south, mez]
      co % a(1:nx, 1:1, 1:nz) = buf3d

      buf2d = co % a(1:nx, ny, 1:nz)[mex, south, mez]
      co % a(1:nx, 1, 1:nz) = buf2d
   end if
   sync all

   deallocate(co % a, buf2d, buf3d)

   if (any(shape2d /= shape(buf2d)) .or. any(shape3d /= shape(buf3d))) then
      write(*, *) 'Test failed!'
      error stop 5
   else
      write(*, *) 'Test passed.'
   end if

   ! Regression would cause error message:
   ! Fortran runtime error on image <...>: libcaf_mpi::caf_get_by_ref(): rank out of range.
end program


text:
The http://www.city.davis.ca.us/meetings/agenda.cfm?c17 Tree Commission is a City Commissions City Commission that establishes rules and regulations relating to the planting, care and maintenance of trees and other plants which overhang public streets and makes recommendations to the City Council on all Street Trees street tree removal requests.

The Tree Commission provides guidance to the Urban Forest Manager and to the City Council regarding tree removal and replacement requests.

The Tree Commission is charged to recommend the removal of a City tree on a casebycase basis for the following reasons: Poor health, disease, exceedingly slow growth, large scale limb failure, and decay.

The Commission also considers hazardous conditions that are caused by the street tree and cannot be mitigated without the removal of the tree.

The Tree Commission does not have the authority to recommend the removal of a City Tree for its tree litter (leaves, needles, pine cones, etc.).

All Tree Commission decisions can be appealed to the City Council for their consideration.


text:
! enddo-plus-label.f90 -- statement labels on 'end do's move to subsequent assignment statements

program prog

      integer :: x, y, z

100   do x = 1, 2
101     y = 1
102       do y = 3, 4
103         z = 2
104       end do
105     y = 3
106   end do

      z = 4
      z = 5

end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE ASPRO (DMAP,VAR,NVAR,OBITS,SOL)
C
C     THIS CODE  PERFORMS THE  ROUTINE PROCESSING OF THE  DMAP ALTERS
C     FOR ASDMAP.  KEY TABLES ARE-
C
C      DMAP  -  RAW  18 WORD PER CARD BCD DATA ON INPUT, VARIABLE
C               CHARACTERS ARE ADDED AND  FIELDS AND CARDS ARE DELETED
C               DEPENDING ON USER INPUT IN  VAR(IABLE) AND OPTION FLAGS.
C
C      VAR      CONTROL DATA AND USER INPUT DATA, 3 WORDS, BCD + DATA
C
C      PTBS     POSITIONS-TO-BE-SET TABLE, CONTENTS-PER ENTRY
C
C                    1   CARD NUMBER IN DMAP
C                    2   FIRST CHARACTER OF MODIFIED FIELD
C                    3   FIRST CHARACTER FOR ADDED VARIABLE
C                    4   NUMBER OF VARIABLE CHARACTERS
C                    5   KEY OF VARIABLE TO BE INSERTED
C                    6   MATRIX OPTION FLAG , 1= K, 2=M, 4=P  ETC
C                    7   OUTPUT CONTROL FLAG, AVOIDS SAME DATA BLOCK
C                        OUTPUT FROM TWO MODULES
C
C      OCT      OPTIONAL CARDS TABLE - EACH ENTRY =
C                   DMAP CARD NO. , DELETE BITS ,  KEEP BITS
C
C      OBITS -  BITS ARE ON FOR REQUIRED MATRICES  =  SUM OF NUMBERS
C                   K=1 , M=2 , P=4 , PA=8 , B=16 , K4=32
C
      EXTERNAL        ANDF
      LOGICAL         RMV,RMVALL
      INTEGER         ANDF,DBS(2,50),DMAP(18,1),FLAG,II(2),NAME(4),
     1                OBITS,OCT(3,50),PTBS(7,200),VAR(3,200),VWORD,SOL,
     2                ALTER,BLANK,AST,SLAS,OBALL,RFMASK(40)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25
      COMMON /XMSSG / UFM,UWM,UIM,SFM
      COMMON /ASDBD / IRDM,NRDM,IXTRA,NXTRA,IOCT,NOCT,IPTBS,NPTBS,
     1                IPH,NPH,IDAT(1)
      COMMON /ZZZZZZ/ SBD(2)
      COMMON /SYSTEM/ IDUM1,IOUT,NOGO
      EQUIVALENCE     (NDBS,SBD(1)),(DBS(1,1),SBD(2))
      DATA    ALTER / 4HALTE /,      BLANK  / 4H     /
      DATA    AST   / 4H*    /,      SLAS   / 4H/    /
      DATA    RFMASK/ 65536,131072,262144,0,0,0,0,524288,1048576,31*0 /
      DATA    OBALL / 63     /
C
      RMVALL = .TRUE.
      NXDEL  = 0
      NOLD   = 0
C
C     DELETE CARDS USING OCT TABLE.
C
      IF (NOCT .EQ. 0) GO TO 45
      M = IOCT - 1
      DO 10 J = 1,NOCT
      DO 10 I = 1,3
      M = M + 1
      OCT(I,J) = IDAT(M)
   10 CONTINUE
      DO 40 I = 1,NOCT
      ICD = OCT(1,I)
      IF (OCT(3,I) .EQ. 0) GO TO 20
      IF (ANDF(OCT(3,I),OBITS)       .EQ. 0) GO TO 35
   20 IF (ANDF(OCT(2,I),RFMASK(SOL)) .EQ. 0) GO TO 40
   35 DMAP(1,ICD) = -1
   40 CONTINUE
   45 IF (NPTBS .EQ. 0) GO TO 2000
      M = IPTBS - 1
      DO 46 J = 1,NPTBS
      DO 46 I = 1,7
      M = M + 1
      PTBS(I,J) = IDAT(M)
   46 CONTINUE
      DO 1000 I = 1,NPTBS
      ICD = PTBS(1,I)
      IF (DMAP(1,ICD) .EQ. -1) GO TO 1000
      IF (ICD .EQ. 0) GO TO 1000
      RMV = .FALSE.
C
C     CHECK IF  OPTION IS ON
C
      KOPT = PTBS(6,I)
      IF (ANDF(KOPT,OBITS) .EQ. 0) RMV = .TRUE.
      IF (ANDF(KOPT,OBALL) .EQ. 0) RMV = .FALSE.
      IF (ANDF(KOPT,RFMASK(SOL)) .NE. 0) RMV = .TRUE.
      NCHAR = PTBS(4,I)
      NC    = 0
      FLAG  = 0
      VWORD = BLANK
      ICOL  = PTBS(3,I)
C
C     FIND  VARIABLE  IF  REQUIRED
C
      IF (RMV) GO TO 300
      KEY = PTBS(5,I)
      I3  = NVAR/3
      DO 60 J = 1,I3
      IF (VAR(1,J) .EQ. KEY) GO TO 70
      IF (KEY.EQ.J .AND. VAR(1,J).EQ.ALTER) GO TO 450
   60 CONTINUE
C
C     VARIABLE HAS NOT BEEN SET
C
      VWORD = BLANK
      RMV   = .TRUE.
      GO TO 300
C
C     VARIABLE  IS FOUND , IT IS IN VAR(2,J) AND/OR VAR(3,J)
C
   70 VWORD   = VAR(2,J)
      NAME(1) = VAR(2,J)
      NAME(2) = VAR(3,J)
C
C     TEST FOR REAL OR INTEGER
C
      IF (VWORD .EQ.  0) GO TO 300
      IF (VWORD .EQ. -1) GO TO 74
      IF (VWORD .EQ. -2) GO TO 3010
C
C     WORD IS REAL (TEMPORARY ERROR)
C
      GO TO 75
C
C     WORD IS AN INTEGER
C
   74 NAME(1) = NAME (2)
      NAME(2) = 0
      FLAG    = 1
   75 IF (PTBS(7,I) .EQ. 0) GO TO 500
      NC    = PTBS(3,I) - PTBS(2,I)
      II(1) = BLANK
      II(2) = BLANK
      IF (NC .GT. 0) GO TO 80
      IF (NC .LT. 0) GO TO 3010
      II(1) = NAME(1)
      II(2) = NAME(2)
      GO TO 100
C
C     CONSTRUCT WHOLE DATA BLOCK NAME
C
   80 CALL PULL (DMAP(1,ICD),II,PTBS(2,I),NC,0)
      CALL PUSH (NAME,II,NC+1,NCHAR,FLAG)
C
C     CHECK OUTPUT DATA BLOCKS AGAINST PREVIOUS OUTPUTS
C
  100 IF (NDBS .EQ. 0) GO TO 142
C
      DO 140 L = 1, NDBS
      IF (II(1).EQ.DBS(1,L) .AND. II(2).EQ.DBS(2,L)) GO TO 150
  140 CONTINUE
  142 IF (PTBS(7,I) .GT. 0) GO TO 200
C
C     VARIABLE IS OK , ADD NAME TO LIST
C
      NDBS = NDBS + 1
      DBS(1,NDBS) = II(1)
      DBS(2,NDBS) = II(2)
      GO TO 500
  150 IF (PTBS(7,I) .GT. 0) GO TO 500
C
C     DATA BLOCK IS OUTPUT, REMOVE IF ALLREADY DEFINED.
C
  200 RMV  =.TRUE.
C
C     REMOVE WHOLE  NAME HERE  , CHECK FOR PARAMETER
C
  300 II(1)   = 0
      NAME(1) = BLANK
      NAME(2) = BLANK
      NAME(3) = BLANK
      NAME(4) = BLANK
      FLAG    = 0
      CALL PULL (DMAP(1,ICD),II,PTBS(2,I),1,0)
      IF (II(1).EQ.SLAS .OR. II(1).EQ.AST) GO TO 500
      ICOL  = PTBS(2,I)
      NCHAR = NCHAR + PTBS(3,I) - PTBS(2,I)
      GO TO 500
C
C     CHECK IF ALTER CARD, OUTPUT AS BCD AND TWO INTEGERS
C
  450 DMAP(1,ICD) = ALTER
      DMAP(2,ICD) = VAR(2,J)
      DMAP(3,ICD) = VAR(3,J)
      RMVALL = .FALSE.
      NXDEL  = 0
      IF (VAR(2,J) .EQ. 0) RMVALL = .TRUE.
      GO TO 910
C
C     ADD VARIABLES TO BCD DMAP
C
  500 CALL PUSH (NAME,DMAP(1,ICD),ICOL,NCHAR,FLAG)
C
      IF (.NOT.RMV) RMVALL = .FALSE.
C
C     IF ALL VARIABLES ARE REMOVED FROM ONE CARD, DELETE THE CARD
C
      NNEW = PTBS(1,I+1)
      IF (ICD .EQ. NNEW) GO TO 1000
C
C     NEXT  COMMAND GOES TO NEW CARD,  CHECK IF CONTINUATION
C
  905 CALL PULL (DMAP(1,ICD+1),II,1,4,0)
C
      IF (II(1) .NE. BLANK) GO TO 910
C
C     CONTINUATION FOUND
C
      NXDEL = NXDEL + 1
      IF (NNEW .EQ. ICD+1) GO TO 1000
      ICD = ICD+1
      GO TO 905
C
C     FINISHED WITH  LOGICAL CARD
C
  910 IF (.NOT.RMVALL) GO TO 940
      DMAP(1,ICD) = -1
      IF (NXDEL .LE. 0) GO TO 1000
      DO 930 L = 1,NXDEL
      J = ICD-L
      DMAP(1,J) = -1
  930 CONTINUE
  940 RMVALL = .TRUE.
C
C     END OF LOOP ON VARIABLE CHARACTERS
C
      NXDEL = 0
C
 1000 CONTINUE
C
C     PROCESS CARDS TO BE DELETED FROM SEQUENCE
C
 2000 IKEEP = 0
      DO 2500 ICD = 1,NRDM
C
      IF (DMAP(1,ICD) .EQ. -1) GO TO 2500
C
C     KEEP CARD
C
      IKEEP = IKEEP + 1
      DO 2450 J = 1,18
      DMAP(J,IKEEP) = DMAP(J,ICD)
 2450 CONTINUE
 2500 CONTINUE
      NRDM = IKEEP
      RETURN
 3010 WRITE  (IOUT,3020) SFM,DMAP(1,ICD)
 3020 FORMAT (A25,' 6010, ILLEGAL VARIABLE TO BE SET IN DMAP STATEMENT',
     1        3X,A4)
C
      NOGO = 1
      RETURN
      END


text:
!-------------------------------------------------------------------------------
!
!+  Program ll2ico
!
!-------------------------------------------------------------------------------
program ll2ico
  !-----------------------------------------------------------------------------
  !
  !++ Description:
  !       This program is the conversion program from lat-lon data
  !       to icosahedral data.
  !
  !++ Current Code Owner: H.Tomita
  !
  !++ History:
  !      Version   Date       Comment
  !      -----------------------------------------------------------------------
  !      0.00      09-??-??   H.Tomita  : Add this program.
  !                12-06-04   R.Yoshida : add direct access output, and add "kmax=k+2"
  !                12-06-07   R.Yoshida : allow large number of record length
  !      -----------------------------------------------------------------------
  !
  !-----------------------------------------------------------------------------
  !
  !++ Used modules (shared)
  !
  use mod_adm, only : &
       ADM_NSYS,      &
       ADM_GALL_PL,   &
       ADM_LALL_PL
  use mod_grd, only : &
       GRD_XDIR,      &
       GRD_YDIR,      &
       GRD_ZDIR
  use mod_misc, only : &
       misc_get_available_fid,&
       misc_make_idstr, &
       MISC_get_latlon

  implicit none

  character(128) :: lldata_fname = 'restart_LLGDATA.dat'
  character(128) :: lldata_dim_fname = 'LLGDATA_dimension.info'
  !
  character(128) :: hgrid_fname = 'grid'
  character(128) :: output_data = 'rst'
  integer :: glevel = 5
  integer :: rlevel = 1
  integer :: num_of_rgn, num_of_gp

  integer :: ctl_fid, fid, fid_grd,fid_odat
  integer :: i,j,k,l,n,t
  character(ADM_NSYS) :: z_coordinate_type
  real(8), allocatable :: lon(:), lat(:), zh(:), z(:), dz(:)
  real(8), allocatable :: data(:,:,:)
  integer :: imax,jmax,kall, tmax, time_increment, kmax
  character(128) :: fname
  !
  integer :: ierr
  real(8),allocatable :: x_ico(:,:),lat_ico(:), lon_ico(:), data_ico(:,:)
  real(8),allocatable :: x_ico_pl(:,:,:),lat_ico_pl(:,:), lon_ico_pl(:,:), data_ico_pl(:,:,:)
  integer :: gall_1d
  real(8) :: dlon, dlat
  real(8) :: a,b,c,d
  !
  integer(8) :: irecl                        ! for long record length ! 12/06/07 R.Yoshida
  character(20) :: access_outdata = 'direct' ! 'direct' or 'sequential' ! 12/06/01 R.Yoshida
  logical :: read_conf = .false.

  namelist / ll2icoparam / &
       lldata_fname,       &
       lldata_dim_fname,   &
       hgrid_fname,        &
       output_data,        &
       access_outdata,     &
       glevel,             &
       rlevel,             &
       read_conf
  !
  !--- open and read configuration file
  ctl_fid = 10
  open(CTL_FID,               &
       file='ll2ico.cnf',     &
       form='formatted',      &
       status='old'           &
       )
  read(CTL_FID,nml=ll2icoparam,iostat=ierr)
  write(*,nml=ll2icoparam)
  if(ierr<0) then
     write(*,*) &
          'Msg : Prg[ll2ico]'
     write(*,*) &
          ' *** Not found namelist.'
     write(*,*) &
          ' *** Use default values.'
  else if(ierr>0) then
     write(*,*) &
          'Msg : Prg[ll2ico]'
     write(*,*) &
          ' *** WARNING : Not appropriate names in namelist!! CHECK!!'
  end if
  close(ctl_fid)
  !
  !---- read lat-lon-z-t dimension data
  fid = MISC_get_available_fid()
  open(fid,file=lldata_dim_fname,form='unformatted',status='old')
  !
  !------ longitude
  read(fid) imax
  allocate(lon(-1:imax+2))
  read(fid) lon(1:imax)
  !
  !------ latitude
  read(fid) jmax
  allocate(lat(jmax))
  read(fid) lat(:)
  !
  !------ z
  read(fid) z_coordinate_type
  read(fid) kall
  kmax = kall + 2       ! 12/06/04 R.Yoshida
  !allocate(zh(kall+1))  ! 12/06/04 R.Yoshida
  allocate(z(kmax))
  allocate(zh(kmax))
  read(fid) z(:)
  read(fid) zh(:)
  !
  allocate(dz(kall+1))
  !allocate(z(kall))  ! 12/06/04 R.Yoshida
  do k=1,kall+1
     !z(k) = 0.5D0*(zh(k+1)+zh(k))  ! 12/06/04 R.Yoshida
     dz(k) = abs(zh(k+1)-zh(k))
  end do
  !
  !------ time
  read(fid) tmax
  read(fid) time_increment
  !
  close(fid)
  !
  !--- debug info
  write(*,*) 'data structure :',imax,jmax,kall,kmax
  write(*,*) 'number of data :',tmax
  !
  !--- set incremental value and longitudinal hallo(2)
  dlon=(lon(imax)-lon(1))/(dble(imax)-1.0D0)
  dlat=(lat(jmax)-lat(1))/(dble(jmax)-1.0D0)
  lon(0) = lon(1)-dlon
  lon(-1) = lon(0)-dlon
  lon(imax+1) = lon(imax)+dlon
  lon(imax+2) = lon(imax+1)+dlon
  !
  !--- memory allocation for latlon data
  !allocate(data(-1:imax+2,jmax,kall))
  allocate(data(-1:imax+2,jmax,kmax))    ! 12/06/04 R.Yoshida
  !
  !--- setup icosahedral info and memory allocation
  num_of_rgn=(2**rlevel)*(2**rlevel)*10
  num_of_gp=(2**(glevel-rlevel)+2)*(2**(glevel-rlevel)+2)
  allocate(x_ico(num_of_gp,GRD_XDIR:GRD_ZDIR))
  allocate(lat_ico(num_of_gp))
  allocate(lon_ico(num_of_gp))
!  allocate(data_ico(num_of_gp,kall))
  allocate(data_ico(num_of_gp,kmax))    ! 12/06/04 R.Yoshida

  allocate(x_ico_pl(ADM_GALL_PL,ADM_LALL_PL,GRD_XDIR:GRD_ZDIR))
  allocate(lat_ico_pl(ADM_GALL_PL,ADM_LALL_PL))
  allocate(lon_ico_pl(ADM_GALL_PL,ADM_LALL_PL))
!  allocate(data_ico_pl(ADM_GALL_PL,kall,ADM_LALL_PL))
  allocate(data_ico_pl(ADM_GALL_PL,kmax,ADM_LALL_PL))    ! 12/06/04 R.Yoshida
  !
  !--- open latlon data
  fid = MISC_get_available_fid()
  irecl = int(imax,kind=8)*int(jmax,kind=8)*int(kmax,kind=8)*8_8    ! 12/06/07 R.Yoshida
!  open(fid,file=lldata_fname,form='unformatted',status='old')
!  open(fid,file=lldata_fname,form='unformatted',status='old',access='direct',recl=8*imax*jmax*kall)
  open(fid,file=lldata_fname,form='unformatted',status='old',access='direct',recl=irecl)    ! 12/06/04 R.Yoshida
  !
  write (*,'(A)') "---------------------------"  ! 12/06/01 R.Yoshida
  !
  !--- main loop
  do t=1,tmax
     !
     write (*,'(" Time Loop:",I4)') t   ! 12/06/07 R.Yoshida
     !
     !--- read data of 1 layer and fix to cyclic condition (longitudinal)
     read(fid,rec=t) data(1:imax,:,:)                                 ! 12/06/04 R.Yoshida
     if( read_conf ) write (*,'(" Read Conf. (t=",I4,", k=1)  Max:",F16.8," Min:",F16.8)') t, &  ! 12/06/01 R.Yoshida
     maxval(data(1:imax,:,:)), minval(data(1:imax,:,:))    ! 12/06/04 R.Yoshida
!     do k=1,kall
!        read(fid) data(1:imax,1:jmax,k)
!        if(k == 1)then  ! 12/06/01 R.Yoshida
!           write (*,'(" Read Conf. (t=",I3,", k=1)  Max:",F10.4," Min:",F10.4)') t, maxval(data), minval(data)
!        endif
!     end do
     data(0,:,:) = data(imax,:,:)
     data(-1,:,:) = data(imax-1,:,:)
     data(imax+1,:,:) = data(1,:,:)
     data(imax+2,:,:) = data(2,:,:)
     !
     !--- regular region
     do l = 1, num_of_rgn
        call MISC_make_idstr(fname,trim(hgrid_fname),'rgn',l)
        fid_grd = MISC_get_available_fid()
        open(fid_grd,file=trim(fname),status='old',form='unformatted',iostat=ierr)
        if(ierr/=0) then
           write(*,*) &
                '   *** No grid file.',trim(fname)
           stop
        end if
        read(fid_grd) gall_1d
        if(gall_1d**2/=num_of_gp) then
           write(*,*) 'inconsistent in number of gridpoint!'
           stop
        end if
        !
        read(fid_grd) x_ico(:,GRD_XDIR)
        read(fid_grd) x_ico(:,GRD_YDIR)
        read(fid_grd) x_ico(:,GRD_ZDIR)
        !
        !
        close(fid_grd)
        !
        !------ setup lat-lon on icosahedral gridpoints
        do n=1,num_of_gp
           call MISC_get_latlon(     &
                lat_ico(n), lon_ico(n),    &
                x_ico(n,GRD_XDIR), x_ico(n,GRD_YDIR), x_ico(n,GRD_ZDIR) )
        end do
        !
        !------ bi-linear interpolation
        do n=1, num_of_gp
           i = floor((lon_ico(n)-lon(1))/dlon)+1
           j = min(max(floor((lat_ico(n)-lat(1))/dlat)+1,1),jmax-1)
           a = (lon(i+1)-lon_ico(n))/dlon
           b = (lon_ico(n)-lon(i))/dlon
           c = (lat(j+1)-lat_ico(n))/dlat
           d = (lat_ico(n)-lat(j))/dlat
           !do k=1,kall
           do k=1,kmax       ! 12/06/04 R.Yoshida
              data_ico(n,k) &
                   = c*(a*data(i,j  ,k)+b*data(i+1,j  ,k)) &
                   + d*(a*data(i,j+1,k)+b*data(i+1,j+1,k))
           end do
        end do
        !
        !------ output
        call MISC_make_idstr(fname,trim(output_data),'rgn',l)
        fid_odat = MISC_get_available_fid()
!        open(fid_odat,file=trim(fname),access='append',form='unformatted')
        if ( access_outdata == 'sequential' ) then ! 12/06/01 R.Yoshida
           open(fid_odat,file=trim(fname),access='sequential',position='append', form='unformatted')
           write(fid_odat) data_ico(:,:)
        else
           !open(fid_odat,file=trim(fname),access='direct',recl=8*num_of_gp*kall,form='unformatted')
           open(fid_odat,file=trim(fname),access='direct',recl=8*num_of_gp*kmax,form='unformatted')   ! 12/06/04 R.Yoshida
           write(fid_odat, rec=t) data_ico(:,:)
        endif
        !
        close(fid_odat)
     end do
     !
     !--- pole region
     fname=trim(hgrid_fname)//'.pl'
     fid_grd = MISC_get_available_fid()
     open(fid_grd,file=trim(fname),status='old',form='unformatted',iostat=ierr)
     if(ierr/=0) then
        write(*,*) &
             '   *** No grid file.',trim(fname)
        stop
     end if
     !
     read(fid_grd) x_ico_pl(:,:,:)
     !
     close(fid_grd)
     !
     !------ setup lat-lon on icosahedral gridpoints
     do l=1,ADM_LALL_PL
        do n=1,ADM_GALL_PL
           call MISC_get_latlon(                    &
                lat_ico_pl(n,l), lon_ico_pl(n,l),    &
                x_ico_pl(n,l,GRD_XDIR), x_ico_pl(n,l,GRD_YDIR), x_ico_pl(n,l,GRD_ZDIR) )
        end do
     end do
     !
     !------ bi-linear interpolation
     do l=1,ADM_LALL_PL
        do n=1,ADM_GALL_PL
           i = floor((lon_ico_pl(n,l)-lon(1))/dlon)+1
           j = min(max(floor((lat_ico_pl(n,l)-lat(1))/dlat)+1,1),jmax-1)
           a = (lon(i+1)-lon_ico_pl(n,l))/dlon
           b = (lon_ico_pl(n,l)-lon(i))/dlon
           c = (lat(j+1)-lat_ico_pl(n,l))/dlat
           d = (lat_ico_pl(n,l)-lat(j))/dlat
           !do k=1,kall
           do k=1,kmax      ! 12/06/04 R.Yoshida
              data_ico_pl(n,k,l) &
                   = c*(a*data(i,j  ,k)+b*data(i+1,j  ,k)) &
                   + d*(a*data(i,j+1,k)+b*data(i+1,j+1,k))
           end do
        end do
     end do
     !
     !------ output
     fname=trim(output_data)//'.pl'
     fid_odat = MISC_get_available_fid()
!     open(fid_odat,file=trim(fname),access='append',form='unformatted')
     open(fid_odat,file=trim(fname),access='sequential',position='append',form='unformatted')
     write(fid_odat) data_ico_pl(:,:,:)
     close(fid_odat)
  end do
  !
  close(fid)

end program ll2ico

max_stars_repo_name: fiber-miniapp/nicam-dc-mini

max_issues_repo_name: fiber-miniapp/nicam-dc-mini

max_forks_repo_name: fiber-miniapp/nicam-dc-mini


text:
#include <default_macro.h>
 module utility_mod

!***********************************************************************
!
!     LBsoft module containing generic supporting subroutines called
!     by different modules
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

 use version_mod, only : idrank,or_world_larr
 use error_mod

 implicit none

 private

 integer, public, save :: nlbuffservice=0
 integer, public, save :: nibuffservice=0
 integer, public, save :: nbuffservice=0
 integer, public, save :: nbdf=0
 integer, public, save, protected :: nbuffservice3d=0
 integer, public, save :: nbuffservice_x=0
 integer, public, save :: nbuffservice_y=0
 integer, public, save :: nbuffservice_z=0
 logical, public, save :: linit_seed=.false.
 logical, public, save :: ltest_mode=.false.
 logical, public, allocatable, save :: lbuffservice(:)
 integer, public, allocatable, save :: ibuffservice(:)
 real(kind=PRC), public, allocatable, save :: buffservice(:)
 real(kind=PRC), public, allocatable, save :: xdf(:),ydf(:),zdf(:)
 real(kind=PRC), public, allocatable, save :: buffservice3d(:,:,:)

 real(kind=PRC), public, parameter :: &
  Pi=real(3.141592653589793238462643383279502884d0,kind=PRC)

 real(kind=PRC), public, parameter :: conv_rad=Pi/real(180.d0,kind=PRC)

 public :: allocate_array_lbuffservice
 public :: allocate_array_ibuffservice
 public :: allocate_array_buffservice
 public :: allocate_array_bdf
 public :: allocate_array_buffservice3d
 public :: init_random_seed
 public :: gauss
 public :: modulvec
 public :: dot
 public :: cross
 public :: xcross
 public :: ycross
 public :: zcross
 public :: sig
 public :: write_fmtnumb
 public :: write_fmtnumb8
 public :: get_prntime
 public :: rand_noseeded
 public :: gauss_noseeded
 public :: dcell
 public :: invert
 public :: int_cube_sphere
 public :: openLogFile
 public :: fcut
 public :: space_fmtnumb
 public :: space_fmtnumb12
 public :: test_little_endian

 contains

 subroutine allocate_array_lbuffservice(imiomax)

!***********************************************************************
!
!     LBsoft subroutine for reallocating the service array buff
!     which is used within this module
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification july 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: imiomax

  if(nlbuffservice/=0)then
    if(imiomax>nlbuffservice)then
      deallocate(lbuffservice)
      nlbuffservice=imiomax+100
      allocate(lbuffservice(0:nlbuffservice))
    endif
  else
    nlbuffservice=imiomax+100
    allocate(lbuffservice(0:nlbuffservice))
  endif

  return

 end subroutine allocate_array_lbuffservice

 subroutine allocate_array_ibuffservice(imiomax)

!***********************************************************************
!
!     LBsoft subroutine for reallocating the service array buff
!     which is used within this module
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: imiomax

  if(nibuffservice/=0)then
    if(imiomax>nibuffservice)then
      deallocate(ibuffservice)
      nibuffservice=imiomax+100
      allocate(ibuffservice(0:nibuffservice))
    endif
  else
    nibuffservice=imiomax+100
    allocate(ibuffservice(0:nibuffservice))
  endif

  return

 end subroutine allocate_array_ibuffservice

 subroutine allocate_array_buffservice(imiomax)

!***********************************************************************
!
!     LBsoft subroutine for reallocating the service array buff
!     which is used within this module
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: imiomax

  if(nbuffservice/=0)then
    if(imiomax>nbuffservice)then
      deallocate(buffservice)
      nbuffservice=imiomax+100
      allocate(buffservice(0:nbuffservice))
    endif
  else
    nbuffservice=imiomax+100
    allocate(buffservice(0:nbuffservice))
  endif

  return

 end subroutine allocate_array_buffservice

 subroutine allocate_array_bdf(imiomax)

!***********************************************************************
!
!     LBsoft subroutine for reallocating the service array buff
!     which is used outside this module for particles
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification October 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: imiomax

  if(nbdf/=0)then
    if(imiomax>nbdf)then
      deallocate(xdf,ydf,zdf)
      nbdf=imiomax+100
      allocate(xdf(1:nbdf),ydf(1:nbdf),zdf(1:nbdf))
    endif
  else
    nbdf=imiomax+100
    allocate(xdf(1:nbdf),ydf(1:nbdf),zdf(1:nbdf))
  endif

  return

 end subroutine allocate_array_bdf

 subroutine allocate_array_buffservice3d(imiomin_x,imiomax_x,imiomin_y, &
  imiomax_y,imiomin_z,imiomax_z)

!***********************************************************************
!
!     LBsoft subroutine for reallocating the service array buff
!     which is used within this module
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: imiomin_x,imiomax_x,imiomin_y, &
   imiomax_y,imiomin_z,imiomax_z

  integer :: ierr
  logical, dimension(1) :: ltest=.false.

  if(allocated(buffservice3d))then
    if(imiomax_x>nbuffservice_x .or. imiomax_y>nbuffservice_y .or. &
     imiomax_z>nbuffservice_z)then
      deallocate(buffservice3d)
      nbuffservice_x=imiomax_x+100
      nbuffservice_y=imiomax_y+100
      nbuffservice_z=imiomax_z+100
      allocate(buffservice3d(imiomin_x:nbuffservice_x,imiomin_y:nbuffservice_y, &
       imiomin_z:nbuffservice_z),stat=ierr)
      ltest=.false.
      if(ierr.ne.0)then
        call warning(2,dble(1))
        ltest=.true.
      endif
    endif
  else
    nbuffservice_x=imiomax_x+100
    nbuffservice_y=imiomax_y+100
    nbuffservice_z=imiomax_z+100
    allocate(buffservice3d(imiomin_x:nbuffservice_x,imiomin_y:nbuffservice_y, &
     imiomin_z:nbuffservice_z),stat=ierr)
    ltest=.false.
    if(ierr.ne.0)then
      call warning(2,dble(1))
      ltest=.true.
    endif
  endif

  nbuffservice3d=(nbuffservice_x-imiomin_x+1)* &
   (nbuffservice_y-imiomin_y+1)*(nbuffservice_z-imiomin_z+1)

  call or_world_larr(ltest,1)
  if(ltest(1))call error(8)

  return

 end subroutine allocate_array_buffservice3d

 subroutine init_random_seed(myseed)

!***********************************************************************
!
!     LBsoft subroutine for initialising the random generator
!     by the seed given in input file or by a random seed
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer,intent(in),optional :: myseed
  integer :: i, n, clock

  integer, allocatable :: seed(:)

  call random_seed(size = n)

  allocate(seed(n))

  if(present(myseed))then
!   If the seed is given in input
    seed = myseed*(idrank+1) + 37 * (/ (i - 1, i = 1, n) /)

  else
!   If the seed is not given in input it is generated by the clock
    call system_clock(count=clock)

    seed = clock*(idrank+1) + 37 * (/ (i - 1, i = 1, n) /)

  endif

  call random_seed(put = seed)

  deallocate(seed)

  return

 end subroutine init_random_seed

 pure function rand_noseeded(i,j,k,l)

!***********************************************************************
!
!     LBsoft random number generator based on the universal
!     random number generator of marsaglia et Al.
!     it is delivered in a form which is depending on a seed squence
!     provided at each call in order to ensure the MPI independence
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer, intent(in) :: i,j,k,l
  integer :: isub,jsub,ksub,lsub,msub
  integer ::ii,jj
  real(4) :: s,t,u33,u97,csub,uni
  real(kind=PRC) :: rand_noseeded

  real(4), parameter :: c =  362436.0/16777216.0
  real(4), parameter :: cd= 7654321.0/16777216.0
  real(4), parameter :: cm=16777213.0/16777216.0

! initial values of i,j,k must be in range 1 to 178 (not all 1)
! initial value of l must be in range 0 to 168.

  isub=mod(i,178)
  isub=isub+1

  jsub=mod(j,178)
  jsub=jsub+1

  ksub=mod(k,178)
  ksub=ksub+1

  lsub=mod(l,169)

! initialization on fly

  ii=97
  s=0.0
  t=0.5
  do jj=1,24
    msub=mod(mod(isub*jsub,179)*ksub,179)
    isub=jsub
    jsub=ksub
    ksub=msub
    lsub=mod(53*lsub+1,169)
    if(mod(lsub*msub,64).ge.32)s=s+t
    t=0.5*t
  enddo
  u97=s

  ii=33
  s=0.0
  t=0.5
  do jj=1,24
    msub=mod(mod(isub*jsub,179)*ksub,179)
    isub=jsub
    jsub=ksub
    ksub=msub
    lsub=mod(53*lsub+1,169)
    if(mod(lsub*msub,64).ge.32)s=s+t
    t=0.5*t
  enddo
  u33=s
  uni=u97-u33
  if(uni.lt.0.0)uni=uni+1.0
  csub=c-cd
  if(csub.lt.0.0)csub=csub+cm
  uni=uni-csub
  if(uni.lt.0.0)uni=uni+1.0
  rand_noseeded=real(uni,kind=PRC)

  return

 end function rand_noseeded

 function gauss()

!***********************************************************************
!
!     LBsoft subroutine for generating random number normally
!     distributed by the Box-Muller transformation
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification January 2020
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: gauss
  real(kind=PRC) :: dtemp1,dtemp2

#if PRC==8
  real(kind=PRC), parameter :: mylimit=1.d-100
#elif PRC==4
  real(kind=PRC), parameter :: mylimit=1.e-50
#endif

  call random_number(dtemp1)
  call random_number(dtemp2)

  dtemp1=dtemp1*(ONE-mylimit)+mylimit

! Box-Muller transformation
  gauss=sqrt(- TWO *log(dtemp1))*cos(TWO*pi*dtemp2)

  return

 end function gauss

 pure function gauss_noseeded(i,j,k,l)

!***********************************************************************
!
!     LBsoft subroutine for generating random number normally
!     distributed by the Box-Muller transformation without a seed
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification January 2020
!
!***********************************************************************

  implicit none
  integer, intent(in) :: i,j,k,l
  integer :: kk
  real(kind=PRC) :: gauss_noseeded
  real(kind=PRC) :: dtemp1,dtemp2

#if PRC==8
  real(kind=PRC), parameter :: mylimit=1.d-100
#elif PRC==4
  real(kind=PRC), parameter :: mylimit=1.e-50
#endif

  dtemp1=rand_noseeded(i,j,k,l)
  kk=nint(dtemp1*FIFTY)
  dtemp2=rand_noseeded(i,j,k,l+kk)

  dtemp1=dtemp1*(ONE-mylimit)+mylimit

! Box-Muller transformation
  gauss_noseeded=sqrt(- TWO *log(dtemp1))*cos(TWO*pi*dtemp2)

  return

 end function gauss_noseeded

 pure function modulvec(a)

!***********************************************************************
!
!     LBsoft function for computing the module of a vector
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: modulvec
  real(kind=PRC), dimension(3), intent(in) :: a

  modulvec = dsqrt(a(1)** TWO + a(2)** TWO + a(3)** TWO )

  return

 end function modulvec

 pure function cross(a,b)

!***********************************************************************
!
!     LBsoft function for computing the cross product of two vectors
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC), dimension(3) :: cross
  real(kind=PRC), dimension(3), intent(in) :: a, b

  cross(1) = a(2) * b(3) - a(3) * b(2)
  cross(2) = a(3) * b(1) - a(1) * b(3)
  cross(3) = a(1) * b(2) - a(2) * b(1)

  return

 end function cross

 pure function xcross(a,b)

!***********************************************************************
!
!     LBsoft function for computing the cross product of two vectors
!     along x
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification November 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: xcross
  real(kind=PRC), dimension(3), intent(in) :: a, b

  xcross = a(2) * b(3) - a(3) * b(2)

  return

 end function xcross

 pure function ycross(a,b)

!***********************************************************************
!
!     LBsoft function for computing the cross product of two vectors
!     along y
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification November 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: ycross
  real(kind=PRC), dimension(3), intent(in) :: a, b

  ycross = a(3) * b(1) - a(1) * b(3)

  return

 end function ycross

 pure function zcross(a,b)

!***********************************************************************
!
!     LBsoft function for computing the cross product of two vectors
!     along z
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification November 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: zcross
  real(kind=PRC), dimension(3), intent(in) :: a, b

  zcross = a(1) * b(2) - a(2) * b(1)

  return

 end function zcross

 pure function int_cube_sphere(rdim,xo,yo,zo,side,xc,yc,zc,ires)

!***********************************************************************
!
!     LBsoft function for computing the volume intersection
!     between a cube and a sphere
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification November 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: int_cube_sphere
  real(kind=PRC), intent(in) :: rdim,xo,yo,zo
  real(kind=PRC), intent(in) :: side,xc,yc,zc
  integer, intent(in) :: ires

  real(kind=PRC) :: myres,sqrdim,xmin,ymin,zmin,xx,yy,zz,rdist,mysum
  integer :: i,j,k,l

  sqrdim=rdim**TWO

  myres=side/real(ires,kind=PRC)
  xmin=xc+(myres-side)/TWO
  ymin=yc+(myres-side)/TWO
  zmin=zc+(myres-side)/TWO
  mysum=ZERO
  do k=1,ires
    zz=zmin+real(k-1,kind=PRC)*myres
    do j=1,ires
      yy=ymin+real(j-1,kind=PRC)*myres
      do i=1,ires
        xx=xmin+real(i-1,kind=PRC)*myres
        rdist=(xx-xo)**TWO+(yy-yo)**TWO+(zz-zo)**TWO
        if(rdist>sqrdim)then
          mysum=mysum+ONE
        endif
      enddo
    enddo
  enddo

  int_cube_sphere=mysum/real(ires*ires*ires,kind=PRC)

  return

 end function int_cube_sphere

 pure function dot(a,b)

!***********************************************************************
!
!     LBsoft function for computing the dot product of two vectors
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC) :: dot
  real(kind=PRC), dimension(3), intent(in) :: a, b

  dot=dot_product(a,b)

  return

 end function dot

 pure function sig(num)

!***********************************************************************
!
!     LBsoft function for returning the sign of a floating number
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC), intent(in):: num

  real(kind=PRC) :: sig

  if(num>0)then
    sig= ONE
  elseif(num==0)then
    sig= ZERO
  else
    sig= - ONE
  endif

  return

 end function sig

 function dimenumb(inum)

!***********************************************************************
!
!     LBsoft function for returning the number of digits
!     of an integer number
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer,intent(in) :: inum
  integer :: dimenumb
  integer :: i
  real(kind=PRC) :: tmp

  i=1
  tmp=real(inum,kind=PRC)
  do
    if(tmp< TEN )exit
    i=i+1
    tmp=tmp/ TEN
  enddo

  dimenumb=i

  return

 end function dimenumb

 function write_fmtnumb(inum)

!***********************************************************************
!
!     LBsoft function for returning the string of six characters
!     with integer digits and leading zeros to the left
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer,intent(in) :: inum
  character(len=6) :: write_fmtnumb
  integer :: numdigit,irest
  real(kind=PRC) :: tmp
  character(len=22) :: cnumberlabel

  numdigit=dimenumb(inum)
  irest=6-numdigit
  if(irest>0)then
    write(cnumberlabel,"(a,i8,a,i8,a)")"(a",irest,",i",numdigit,")"
    write(write_fmtnumb,fmt=cnumberlabel)repeat('0',irest),inum
  else
    write(cnumberlabel,"(a,i8,a)")"(i",numdigit,")"
    write(write_fmtnumb,fmt=cnumberlabel)inum
  endif

  return

 end function write_fmtnumb

 function write_fmtnumb8(inum)

!***********************************************************************
!
!     LBsoft function for returning the string of eight characters
!     with integer digits and leading zeros to the left
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  integer,intent(in) :: inum
  character(len=8) :: write_fmtnumb8
  integer :: numdigit,irest
  real(kind=PRC) :: tmp
  character(len=22) :: cnumberlabel

  numdigit=dimenumb(inum)
  irest=8-numdigit
  if(irest>0)then
    write(cnumberlabel,"(a,i8,a,i8,a)")"(a",irest,",i",numdigit,")"
    write(write_fmtnumb8,fmt=cnumberlabel)repeat('0',irest),inum
  else
    write(cnumberlabel,"(a,i8,a)")"(i",numdigit,")"
    write(write_fmtnumb8,fmt=cnumberlabel)inum
  endif

  return

 end function write_fmtnumb8

 subroutine get_prntime(hms,timelp,prntim)

!***********************************************************************
!
!     LBsoft subroutine for casting cpu elapsed time into days, hours,
!     minutes and seconds for printing (input timelp in seconds)
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification July 2018
!
!***********************************************************************

  implicit none

  character(len=1), intent(out) :: hms
  real(kind=PRC), intent(in) :: timelp
  real(kind=PRC), intent(out) :: prntim

  if(timelp.ge. real(8.64d4,kind=PRC))then
    hms='d'
    prntim=timelp/real(8.64d4,kind=PRC)
  elseif(timelp.ge.real(3.6d3,kind=PRC))then
    hms='h'
    prntim=timelp/real(3.6d3,kind=PRC)
  elseif(timelp.ge.real(6.0d1,kind=PRC))then
    hms='m'
    prntim=timelp/real(6.0d1,kind=PRC)
  else
    hms='s'
    prntim=timelp
  endif

  return

 end subroutine get_prntime

 subroutine dcell(aaa,bbb)

!***********************************************************************
!
!     LBsoft subroutine to calculate the dimensional properties of
!     a simulation cell specified by the input matrix aaa.
!     the results are returned in the array bbb, with :
!
!     bbb(1 to 3) - lengths of cell vectors
!     bbb(4 to 6) - cosines of cell angles
!     bbb(7 to 9) - perpendicular cell widths
!     bbb(10)     - cell volume
!
!     modified from DL Protein
!     original copyright - daresbury laboratory 1992
!     original author    - w. smith       april 1992
!
!***********************************************************************

  real(kind=PRC), intent(in) :: aaa(9)
  real(kind=PRC), intent(out) :: bbb(10)
  real(kind=PRC) :: axb1,axb2,axb3,bxc1,bxc2,bxc3,cxa1,cxa2,cxa3

! calculate lengths of cell vectors

  bbb(1)=sqrt(aaa(1)*aaa(1)+aaa(2)*aaa(2)+aaa(3)*aaa(3))
  bbb(2)=sqrt(aaa(4)*aaa(4)+aaa(5)*aaa(5)+aaa(6)*aaa(6))
  bbb(3)=sqrt(aaa(7)*aaa(7)+aaa(8)*aaa(8)+aaa(9)*aaa(9))

! calculate cosines of cell angles

  bbb(4)=(aaa(1)*aaa(4)+aaa(2)*aaa(5)+aaa(3)*aaa(6))/(bbb(1)*bbb(2))
  bbb(5)=(aaa(1)*aaa(7)+aaa(2)*aaa(8)+aaa(3)*aaa(9))/(bbb(1)*bbb(3))
  bbb(6)=(aaa(4)*aaa(7)+aaa(5)*aaa(8)+aaa(6)*aaa(9))/(bbb(2)*bbb(3))

! calculate vector products of cell vectors

  axb1=aaa(2)*aaa(6)-aaa(3)*aaa(5)
  axb2=aaa(3)*aaa(4)-aaa(1)*aaa(6)
  axb3=aaa(1)*aaa(5)-aaa(2)*aaa(4)
  bxc1=aaa(5)*aaa(9)-aaa(6)*aaa(8)
  bxc2=aaa(6)*aaa(7)-aaa(4)*aaa(9)
  bxc3=aaa(4)*aaa(8)-aaa(5)*aaa(7)
  cxa1=aaa(8)*aaa(3)-aaa(2)*aaa(9)
  cxa2=aaa(1)*aaa(9)-aaa(3)*aaa(7)
  cxa3=aaa(2)*aaa(7)-aaa(1)*aaa(8)

! calculate volume of cell

  bbb(10)=abs(aaa(1)*bxc1+aaa(2)*bxc2+aaa(3)*bxc3)

!  calculate cell perpendicular widths
  bbb(7)=bbb(10)/sqrt(bxc1*bxc1+bxc2*bxc2+bxc3*bxc3)
  bbb(8)=bbb(10)/sqrt(cxa1*cxa1+cxa2*cxa2+cxa3*cxa3)
  bbb(9)=bbb(10)/sqrt(axb1*axb1+axb2*axb2+axb3*axb3)

  return

 end subroutine dcell

 subroutine invert(a,b,d)

!***********************************************************************
!
!     LBsoft subroutine to invert a 3 * 3 matrix using cofactors
!
!     modified from DL Protein
!     original copyright - daresbury laboratory 1992
!     original author    - w. smith       april 1992
!
!***********************************************************************

  implicit none

  real(kind=PRC), intent(in) :: a(9)

  real(kind=PRC), intent(out) :: b(9),d
  real(kind=PRC) :: r

! calculate adjoint matrix
  b(1)=a(5)*a(9)-a(6)*a(8)
  b(2)=a(3)*a(8)-a(2)*a(9)
  b(3)=a(2)*a(6)-a(3)*a(5)
  b(4)=a(6)*a(7)-a(4)*a(9)
  b(5)=a(1)*a(9)-a(3)*a(7)
  b(6)=a(3)*a(4)-a(1)*a(6)
  b(7)=a(4)*a(8)-a(5)*a(7)
  b(8)=a(2)*a(7)-a(1)*a(8)
  b(9)=a(1)*a(5)-a(2)*a(4)

! calculate determinant
  d=a(1)*b(1)+a(4)*b(2)+a(7)*b(3)
  r=ZERO
  if(abs(d).gt.ZERO)r=ONE/d

! complete inverse matrix
  b(1)=r*b(1)
  b(2)=r*b(2)
  b(3)=r*b(3)
  b(4)=r*b(4)
  b(5)=r*b(5)
  b(6)=r*b(6)
  b(7)=r*b(7)
  b(8)=r*b(8)
  b(9)=r*b(9)

  return

 end subroutine invert

 subroutine openLogFile(nstep, msg, iounit)
  implicit none
  integer, intent(in) :: nstep,iounit
  character(len=*), intent(in) :: msg
  character(len=120) :: mynamefile

  mynamefile=repeat(' ',120)
  mynamefile=trim(msg)// '.iter'//write_fmtnumb(nstep)//'.'//write_fmtnumb(idrank)//'.dat'
  open(unit=iounit, file=trim(mynamefile), status='replace')
 end subroutine openLogFile

 pure function fcut(r,inner_cut,outer_cut)

!***********************************************************************
!
!     LBsoft function for fading an observable (r) within a given
!     interval
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification January 2018
!
!***********************************************************************

  implicit none

  real(kind=PRC), intent(in) :: r,inner_cut,outer_cut
  real(kind=PRC) :: fcut

  if ( r <= inner_cut ) then
    fcut = ONE
  elseif ( r > outer_cut ) then
      fcut = ZERO
  else
      fcut = HALF*cos((r-inner_cut)*Pi/(outer_cut-inner_cut))+HALF
  endif

  return

 end function fcut

 function space_fmtnumb(inum)

!***********************************************************************
!
!     LBsoft function for returning the string of six characters
!     with integer digits and leading spaces to the left
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification October 2019
!
!***********************************************************************

  implicit none

  integer,intent(in) :: inum
  character(len=6) :: space_fmtnumb
  integer :: numdigit,irest
  real(kind=8) :: tmp
  character(len=22) :: cnumberlabel

  numdigit=dimenumb(inum)
  irest=6-numdigit
  if(irest>0)then
    write(cnumberlabel,"(a,i8,a,i8,a)")"(a",irest,",i",numdigit,")"
    write(space_fmtnumb,fmt=cnumberlabel)repeat(' ',irest),inum
  else
    write(cnumberlabel,"(a,i8,a)")"(i",numdigit,")"
    write(space_fmtnumb,fmt=cnumberlabel)inum
  endif

  return

 end function space_fmtnumb

 function space_fmtnumb12(inum)

!***********************************************************************
!
!     LBsoft function for returning the string of six characters
!     with integer digits and leading TWELVE spaces to the left
!     originally written in JETSPIN by M. Lauricella et al.
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification October 2019
!
!***********************************************************************

  implicit none

  integer,intent(in) :: inum
  character(len=12) :: space_fmtnumb12
  integer :: numdigit,irest
  real(kind=8) :: tmp
  character(len=22) :: cnumberlabel

  numdigit=dimenumb(inum)
  irest=12-numdigit
  if(irest>0)then
    write(cnumberlabel,"(a,i8,a,i8,a)")"(a",irest,",i",numdigit,")"
    write(space_fmtnumb12,fmt=cnumberlabel)repeat(' ',irest),inum
  else
    write(cnumberlabel,"(a,i8,a)")"(i",numdigit,")"
    write(space_fmtnumb12,fmt=cnumberlabel)inum
  endif

  return

 end function space_fmtnumb12

 subroutine test_little_endian(ltest)

!***********************************************************************
!
!     LBsoft subroutine for checking if the computing architecture
!     is working in little-endian or big-endian
!
!     licensed under the 3-Clause BSD License (BSD-3-Clause)
!     author: M. Lauricella
!     last modification October 2019
!
!***********************************************************************

  implicit none
  integer, parameter :: ik1 = selected_int_kind(2)
  integer, parameter :: ik4 = selected_int_kind(9)

  logical, intent(out) :: ltest

  if(btest(transfer(int((/1,0,0,0/),ik1),1_ik4),0)) then
    !it is little endian
    ltest=.true.
  else
    !it is big endian
    ltest=.false.
  end if

  return

 end subroutine test_little_endian

 end module utility_mod

max_stars_repo_name: copmat/LBsoft

max_issues_repo_name: copmat/LBsoft

max_forks_repo_name: copmat/LBsoft


text:
!------------------------------------------------------------------------------
!M+
! NAME:
!       parameters
!
! PURPOSE:
!       Module to hold parameter constants
!
! CATEGORY:
!       NCEP RTM
!
! CALLING SEQUENCE:
!       USE parameters
!
! CREATION HISTORY:
!       Written by:     Yong Chen, CIRA/NESDIS/JCSDA 10-Jan-2008
!                       Yong.Chen@noaa.gov
!
!
!  Copyright (C) 2008 Paul van Delst
!
!M-
!------------------------------------------------------------------------------

MODULE Parameters

  ! ---------------------
  ! Module use statements
  ! ---------------------

  USE  Type_Kinds, ONLY : fp_kind

  ! ------------------
  ! Default visibility
  ! ------------------

  PRIVATE

  ! ----------------------------
  !   modlecule set list
  ! ----------------------------
  ! -------------------------------
  ! The molecular set specification
  ! -------------------------------

  ! -- The number of "molecular sets" recognised
  INTEGER, PUBLIC, PARAMETER :: N_MOLECULE_SETS = 23

  ! -- The names of the allowed molecular sets.
  ! -- These values are used in filenames and
  ! -- other transmittance production program inputs
  CHARACTER( * ), PUBLIC, PARAMETER, DIMENSION( N_MOLECULE_SETS ) :: MOLECULE_SET_TAG = &
    (/ 'mol1            ', &  !   1
       'mol2            ', &  !   2
       'mol3            ', &  !   3
       'mol4            ', &  !   4
       'mol5            ', &  !   5
       'mol6            ', &  !   6
       'mol7            ', &  !   7
       'all_nocontinuum ', &  !   8
       'continua_only   ', &  !   9
       'all_withcontinua', &  !  10
       'wvo             ', &  !  11
       'wet             ', &  !  12
       'dry             ', &  !  13
       'ozo             ', &  !  14
       'wco             ', &  !  15
       'effective_mol1  ', &  ! 101
       'effective_wet   ', &  ! 112
       'effective_dry   ', &  ! 113
       'effective_ozo   ', &  ! 114
       'effective_co    ', &  ! 115
       'effective_ch4   ', &  ! 116
       'effective_co2   ', &  ! 117
       'effective_n2o   ' /)  ! 118

  ! -- The ID values associated with the allowed moleculer sets
  INTEGER, PUBLIC, PARAMETER, DIMENSION( N_MOLECULE_SETS ) :: MOLECULE_SET_TAG_ID = &
    (/  1, &  !  mol1
        2, &  !  mol2
        3, &  !  mol3
        4, &  !  mol4
        5, &  !  mol5
        6, &  !  mol6
        7, &  !  mol7
        8, &  !  all_nocontinuum
        9, &  !  continua_only
       10, &  !  all_withcontinua
       11, &  !  wvo
       12, &  !  wet
       13, &  !  dry
       14, &  !  ozo
       15, &  !  wco
      101, &  !  effective_wet (line only)
      112, &  !  effective_wet
      113, &  !  effective_dry
      114, &  !  effective_ozo
      115, &  !  effective_co
      116, &  !  effective_ch4
      117, &  !  effective_co2
      118 /)  !  effective_n2o

  ! ------------------------------------------------------------
  ! The absorber IDs. Use HITRAN definitions
  ! -----------------------------------------------------------
  INTEGER, PUBLIC, PARAMETER :: H2O_ID  = 1
  INTEGER, PUBLIC, PARAMETER :: CO2_ID  = 2
  INTEGER, PUBLIC, PARAMETER :: O3_ID   = 3
  INTEGER, PUBLIC, PARAMETER :: N2O_ID  = 4
  INTEGER, PUBLIC, PARAMETER :: CO_ID   = 5
  INTEGER, PUBLIC, PARAMETER :: CH4_ID  = 6

  ! -- Define the secant of the zenith angles to be used
  INTEGER, PUBLIC,PARAMETER :: N_ZENITH_ANGLES = 7

  ! --------------------------------------
  ! Number of absorber layers in algorithm
  ! --------------------------------------

  INTEGER, PUBLIC, PARAMETER :: MAX_N_ABSORBER_LAYERS = 300

 ! ----------------------------------------------------------
 ! Cosmic background temperature. Taken from
 ! Numerical limits based on experiment.
  REAL(fp_kind), PUBLIC, PARAMETER :: LIMIT_EXP = 20.0_fp_kind
  REAL(fp_kind), PUBLIC, PARAMETER :: LIMIT_LOG = 4.8e+08_fp_kind   ! EXP( LIMIT_EXP )

  ! --------------------
  ! Numerical parameters
  ! --------------------

  ! -- Numbers
  REAL( fp_kind ), PUBLIC, PARAMETER :: ZERO      = 0.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: ONE       = 1.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: TWO       = 2.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: THREE     = 3.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: FOUR      = 4.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: POINT_25   = 0.25_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: POINT_5   = 0.5_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: POINT_75  = 0.75_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: ONE_POINT_5   = 1.5_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: ONE_POINT_25  = 1.25_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: ONE_POINT_75  = 1.75_fp_kind

  ! -- Precision/tolerance
  REAL( fp_kind ), PUBLIC, PARAMETER :: TOLERANCE = EPSILON( ONE )
  REAL( fp_kind ), PUBLIC, PARAMETER :: INFINITE  = HUGE( ONE )

  ! -- Constant to allow degrees->radians conversion
  REAL( fp_kind ), PUBLIC, PARAMETER :: PI = 3.14159265358979323_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: DEGREES_TO_RADIANS = PI / 180.0_fp_kind

  ! -- Top-Of-Atmosphere pressure in hPa
  REAL( fp_kind ), PUBLIC, PARAMETER :: TOA_PRESSURE = 0.005_fp_kind

  ! -- Reciprocal gravity (scaled by 100 for use with pressure in hPa)
  REAL( fp_kind ), PUBLIC, PARAMETER :: RECIPROCAL_GRAVITY = ONE / 980.665_fp_kind

  ! -- Diffusivity angle secant = ACOS( 3/5 ) in degrees (~53.13)
  REAL( fp_kind ), PUBLIC, PARAMETER :: SECANT_DIFFUSIVITY_ANGLE = 5.0_fp_kind / THREE

  ! -- Maximum solar zenith angle secant definition. Should be determined
  ! -- by the maximum angle secant used in generating the transmittance
  ! -- model coefficients, i.e. a secant of 2.25 => 63.6deg. Users have
  ! -- requested the value be 85deg => secant of ~11.47.
  REAL( fp_kind ), PUBLIC, PARAMETER :: MAX_SOLAR_ANGLE = 85.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: MAX_SECANT_SOLAR_ANGLE = 11.473711738554476_fp_kind

  ! Minimum component and total tau, below which the transmittance is considered
  ! insignificant. Maximum optical path and minimum layer optical depth
  REAL( fp_kind ), PUBLIC, PARAMETER :: TAU_MIN  = 1.0e-10_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: TOTAL_TAU_MIN  = 1.0e-5_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: PATH_MAX = 7.0_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: OD_MIN = 0.000001_fp_kind
  REAL( fp_kind ), PUBLIC, PARAMETER :: OD_MAX = 20.0_fp_kind

  REAL( fp_kind ), PUBLIC, PARAMETER :: RMISS = -2**15 * ONE
  REAL( fp_kind ), PUBLIC, PARAMETER :: HMISS = RMISS / TWO

  REAL( fp_kind ), PUBLIC, PARAMETER :: VirtEmiss_mw  = 0.0_fp_kind   ! MW virtial emissivity
  REAL( fp_kind ), PUBLIC, PARAMETER :: VirtEmiss_ir  = 0.98_fp_kind  ! IR virtial emissivity
  LOGICAL, PUBLIC, PARAMETER         :: Tb_OPTIMAL = .true.

  ! Interferometer WMO_Sensor_IDs
  INTEGER, PUBLIC, PARAMETER :: INTERFEROMETER_LIST(2) = (/221, 620/)  !IASI and CrIS

END MODULE Parameters

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm