SUBROUTINE CON implicit none C **** DEFINES THE MODEL CONSTANTS include "params.h" include "blnk.h" include "vscr.h" include "buff.h" include "rfft.h" include "index.h" include "cons.h" include "strt.h" include "unit.h" include "sechis.h" ! ! Local: real :: prndtl,pi,dlamda,dphi,ds,delt,omega,phi,delta,alpha,s, | expds integer :: n,k,len,i,nbaphy,lim,ifld,maxfld,maxphy,jmaxh,js,j ! DATA ISAVE/3/ ! ! Local: integer :: kphys(nzphys) DATA KPHYS/nzphys*3/ C C **** C **** SET PRANDTL NUMBER C **** DATA PRNDTL/1./ c C **** UNIT NUMBERS IHIST = 12 ISORC = 11 ! IUNIN=8 ! IUNOT=9 IMAG=17 IUAMI = 16 C **** IF FORECAST START THEN SET FLAG TO 0 MODEFC=-1 IF(MODEST.EQ.3) MODEFC=0 C **** INITIALIZE CONTROL PARAMETERS IFRST=0 NERR=0 C **** RESOLUTION PARAMETERS IMAX=ZIMX JMAX=ZJMX KMAX=ZKMX C **** DERIVED CONSTANTS PI=3.14159265358979 DLAMDA=2.*PI/IMAX DPHI=PI/JMAX IF(NHEMI.NE.0) DPHI=PI/(2*JMAX) DS=(ZST-ZSB)/KMAX IMAXH=IMAX/2 IMAXP2=IMAX+2 IMAXP4=IMAX+4 JMAXM1=JMAX-1 JMAXM2=JMAX-2 KMAXM1=KMAX-1 KMAXP1=KMAX+1 LAST=IMAXP4*KMAX-2 NPROG=3*KMAXP1 NPROGP=NPROG+1 C **** BUFFER INFORMATION C **** COMPUTE NUMBER OF FIELDS ON TAPE C **** AND INITIALIZE POINTERS LEN1=IMAXP4 LEN2=IMAXP4*KMAX LEN3=IMAXP4*KMAXP1 LONG=LEN2-4 LBDSK=0 LBSCM=0 LBTAP=1 NCOLS=0 NCBUF = 0 DO 15 N=1,NFLDS IF(KFLDS(N)-2) 11,12,13 11 K=1 LEN=LEN1 GO TO 14 12 K=KMAX LEN=LEN2 GO TO 14 13 K=KMAXP1 LEN=LEN3 14 NDEXA(N+1) = NCBUF LBSCM=LBSCM+LEN NCBUF = NCBUF+K IF(N.LE.NDISK) LBDSK=LBDSK+LEN IF(N.LE.NTAPE) NCOLS=NCOLS+K IF(N.LE.NTAPE) LBTAP=LBTAP+IMAXP2*K 15 CONTINUE c c lbsech = j-slice buffer length for secondary histories c (dependent on fields selected in secflds input) c lbsech = 1 ! secondary hist buffer length per lat if (lusech.gt.0) then do n=1,mxfsech if (len_trim(secflds(n)).gt.0) lbsech = lbsech + imaxp2*kmaxp1 enddo write(6,"('Buffer length for secondary ', + 'histories (lbsech) = ',i6,' (per latitude)')") lbsech endif LBDSKI = ((LBDSK+511)/512)*512 LBPHY=0 NCPHYS=0 DO 20 N=1,NPHYS IF(KPHYS(N)-2) 16,17,18 16 LEN=LEN1 K=1 GO TO 19 17 LEN=LEN2 K=KMAX GO TO 19 18 LEN=LEN3 K=KMAXP1 19 NDEXB(N+1) = NCPHYS NCPHYS = NCPHYS+K LBPHY=LBPHY+LEN 20 CONTINUE C **** COMPUTE COLUMN POINTERS FOR SCM BUFFERS NBADDR(1)=1 DO 23 I=2,8 23 NBADDR(I)=NBADDR(I-1)+NCBUF NBAPHY=NBADDR(8)+NCBUF DO 24 N=1,NPHYS 24 NDEXB(N+1) = NDEXB(N+1)+NBAPHY C **** OUTPUT POINTER INFORMATION LIM=NFLDS+NPHYS+2 write(6,"(/72('-'))") write(6,"('POINTER VALUES FOR EACH FIELD IN COMMON BLOCK ', + 'NDEX:')") write(6,"('KMAX=',i3,' KMAXP1=',i3,' ZSB=',f8.3)") + kmax,kmaxp1,zsb write(6,"('NDEXA: nflds=',i2)") nflds write(6,"('IFLD NAME VALUE KMAX(P1) KPOLE KEQTR')") write(6,"('-------------------------------------------')") do ifld=1,nflds if (kflds(ifld).eq.3) then write(6,"(i4,2x,a8,i7,2x,'KMAXP1',2x,2i6)") ifld, + nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld) elseif (kflds(ifld).eq.2) then write(6,"(i4,2x,a8,i7,2x,'KMAX ',2x,2i6)") ifld, + nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld) else write(6,"(i4,2x,a8,i7,2x,'1 ',2x,2i6)") ifld, + nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld) endif enddo write(6,"(/'NDEXB: nphys=',i2)") nphys write(6,"('IFLD NAME VALUE KMAX(P1)')") write(6,"('-------------------------------')") do ifld=1,nphys if (kphys(ifld).eq.3) then write(6,"(i4,2x,a8,i7,2x,'KMAXP1')") ifld, + nphys_lab(ifld),ndexb(ifld+1) elseif (kphys(ifld).eq.2) then write(6,"(i4,2x,a8,i7,2x,'KMAX ')") ifld, + nphys_lab(ifld),ndexb(ifld+1) else write(6,"(i4,2x,a8,i7,2x,'1 ')") ifld, + nphys_lab(ifld),ndexb(ifld+1) endif enddo write(6,"(72('-')/)") C **** CHECK IF F DIMENSION IS APPROPRIATE MAXFLD=ZFLDX MAXPHY=ZPHYX WRITE(6,1010)MAXFLD,MAXPHY,NCBUF,NCPHYS 1010 FORMAT(//* F ARRAY DIMENSIONED FOR FLDX=*I5,2X,*PHYX=*I5, A*.*/* MODEL NEEDS NCBUF=*I5,2X,*NCPHYS=*I5, B*.*/) IF(NCBUF.LE.MAXFLD.AND.NCPHYS.LE.MAXPHY) AGO TO 35 WRITE(6,1020)NCBUF,NCPHYS 1020 FORMAT(//* EXIT CALLED BY CON*/* NEED FLDX.GE.*I5,2X,*PHYX.GE. A*I5) STOP 35 CONTINUE C **** INITIALIZE FFT ROUTINE CALL SET99(TRIGS,IFAX,IMAX) C **** COMPUTE PHI VARIATION ARRAYS DELT=86400./FLOAT(IMAX) JMAXH=JMAX/2 OMEGA=7.292E-5 JS=-JMAXH IF(NHEMI.NE.0) JS=0 DO 40 J=1, JMAX PHI=(J+JS-.5)*DPHI CS(J)=COS(PHI) SN(J)=SIN(PHI) TN(J)=TAN(PHI) COR(J)=2.*OMEGA*SN(J) 40 CONTINUE CSSP(1)=-CS(2) CSSP(2)=-CS(1) C **** IF HEMI CASE, COS HAS SAME SIGN ACROSS EQUA IF(NHEMI.NE.0) CSSP(1)=CS(2) IF(NHEMI.NE.0) CSSP(2)=CS(1) CSNP(1)=-CS(JMAX) CSNP(2)=-CS(JMAXM1) KUTSP(1) = KUT(2) KUTSP(2) = KUT(1) KUTNP(1) = KUT(JMAX) KUTNP(2) = KUT(JMAXM1) C **** CALCULATE SUN'S DECLINATION C **** (WILL HAVE TO RECALCULATE THIS IF DAY IS ADVANCED) DELTA=ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IIDAY-80)/365.)) C **** GRID CONSTANTS C(1)=DLAMDA C(2)=DPHI C(3)=DS C(6)=2.*C(4) C(7)=1./C(6) C **** FINITE DIFFERENCE OPERATOR CONSTANTS C(10)=2./(3.*DLAMDA) C(11)=1./(12.*DLAMDA) C(12)=2./(3.*DPHI) C(13)=1./(12.*DPHI) C(14)=1./(2.*DLAMDA) C(15)=1./(2.*DPHI) C **** DIFFUSION COEFFICIENT C(16)=0. C(16)=C(16)/6. C(17)=0. C **** INTEGRATION CONSTANTS C(18)=C(3)/(1.+.5*C(3)) C(19)=(1.-.5*C(3))/(1.+.5*C(3)) C(20)=1./DS C(22)=COS(PI/3.) C(23)=4.*PI/(24.*60.*60.) C(24)=28.9 C **** ANNUAL FREQUENCY C(25) = C(23)/(2.*365.25) C **** SHAPIRO SMOOTHER CONSTANT C(26)=3.0E-2 C **** EDDY DIFFUSION C C(27)=2.5E-6 C(27)=5.0E-6 C **** EDDY THERMAL CONDUCTION C(28)=C(27)/PRNDTL C **** TIME SMOOTHING CONSTANT ALPHA=.95 C(30)=ALPHA C(31)=.5*(1.-ALPHA) C **** RATE CONSTANT FOR ATOMIC OXYGEN RECOMBINATION C(44)=0. C ***** SURPLUS HEAT PER EVENT CONVERTED FROM E.V. TO ERGS C(45)=5.11*1.602E-12 C **** PHYSICAL CONSTANTS C **** C EARTHS RADIUS C(51)=6.37122E8 C(52)=1./C(51) C(53)=C(51)*C(51) C GRAVITATIONAL CONSTANT C(54)=870. C(55)=1./C(54) C EARTHS ANGULAR VELOCITY C(56)=OMEGA C GAS CONSTANT C(57)=8.314E7 C **** H*C C(60)=1.9845E-16 C **** F107, F107A C C(61) = F107 C C(62) = F107A C **** G/R C(65)=C(54)/C(57) C **** STANDARD PRESSURE C(81)=5.0E-4 C **** EDDY VISCOSITY C(82)=0. C(83)=C(81)*C(82)/C(57)*EXP(C(3)/2.) C **** BOLTZMANN CONST C(84)=1.38E-16 C **** AVOGRADO NUMBER C(85)=6.023E23 C **** EXP(-.5*DS) C(86)=EXP(-.5*DS) C **** EXP(.5*DS) C(87)=1./C(86) C **** SUN'S DECLINATION: C **** (WILL HAVE TO RECALCULATE IF DAY IS ADVANCED) C(95)=SIN(DELTA) C(96)=COS(DELTA) C **** HEIGHT OF BOTTOM OF MODEL C(97)=97.E5 C **** O+ EDDY DIFFUSION C(100)=0. C **** KO COEFFICIENT OF EDDY VISCOSITY C(109)=0.2 C(110) = PI C **** SAVE ORIGINAL DT VALUE C(111) = C(4) C(112)=SQRT(PI) ! ! Intitialize SFEPS: Default solar flux variation due to orbital ! eccentricity assumes no calendar day advance, thus no variation. ! If calendar day *is* being advanced, SFEPS will be calculated ! for each new day in advnce.f. SFEPS is in common in cons.h. ! SFEPS = 1. C **** DUMP CONSTANTS WRITE(6,70) IMAX,JMAX,KMAX,NFLDS,NPHYS,C(4),DS 70 FORMAT(*1IMAX=*I3,2X,*JMAX=*I2,2X,*KMAX=*I2,2X,*NFLDS=*I3,2X, A*NPHYS=*I3,2X,*DT=*E10.3,2X,*DS=*E10.3) WRITE(6,80) LBDSK,LBSCM,LBPHY,LBTAP 80 FORMAT(* LBDSK=*I6,2X,*LBSCM=*I6,2X,*LBPHY=*I6,2X,*LBTAP=*I6) WRITE(6,95) KUT 95 FORMAT(1X,*KUT = *20I5/(7X,20I5)) WRITE(6,90) CSSP 90 FORMAT(*0 J*,4X,*CS(J)*,7X,*F (J)*7X,*SN(J)*,7X,*TN(J)*// A(3X,E12.3)) WRITE(6,100) (J,CS(J),COR(J),SN(J),TN(J),J=1,JMAX) 100 FORMAT(1H ,I2,4E12.3) WRITE(6,105) CSNP 105 FORMAT(3X,E12.3) WRITE(6,130) C 130 FORMAT(*1C=*10E12.3/(3X,10E12.3)) C **** CALCULATE EXPS AND A ARRAYS S=ZSB+.5*DS EXPS(1)=EXP(-S) EXPDS=EXP(-DS) A(1)=0. DIFK(1)=C(27) DIFT(1)=C(28) XMUE(1) = C(27) DO 1 K=2,KMAX EXPS(K)=EXPS(K-1)*EXPDS A(K)=A(K-1)*EXPDS DIFK(K)=DIFK(K-1)*EXPDS DIFT(K) = DIFT(K-1)*EXPDS XMUE(K) = DIFK(K) 1 CONTINUE DIFK(KMAXP1)=DIFK(KMAX)*EXPDS DIFT(KMAXP1) = DIFT(KMAX)*EXPDS XMUE(KMAXP1) = DIFK(KMAXP1) C **** WRITE(6,135)(K,EXPS(K),A(K),DIFK(K),DIFT(K),K=1,KMAX) 135 FORMAT(*0 K*,5X,*EXPS(K)*,8X,*A(K)*,5X,*DIFK(K)*,5X,*DIFT(K)*// A(1H ,I2,4E12.3)) CALL RATE RETURN END C