SUBROUTINE START implicit none C **** STARTS THE MODEL FROM A SAVE TAPE, C **** HISTORY TAPE, OR DEADSTART include "params.h" include "blnk.h" include "vscr.h" include "buff.h" include "cons.h" include "index.h" include "strt.h" include "unit.h" include "fieldz.h" include "cterp.h" include "trig.h" include "dynphi.h" include "sechis.h" include "fgcom.h" real :: ttbound integer :: nerd COMMON/RUNTIM/TTBOUND,NERD(2*ZJMX+1) include "amie.h" ! ! Local: integer :: stop(3,2),lim,i,k,ios1,istat,ios3,ios4,ios5,ios6,ios7, | ios9,ios,iyda1,nsect,ndayad,iyra,ida1,ida2,iyr4,iyr100,lpyr, | ienda,iyda2,iyd1,iutsmod,iutsmend,nrd,n,istrt,isearch,iend_iyd, | iend_day,iend_hr,iend_min,idelmin,iendonly,ibeg_iyd,ibeg_day, | ibeg_hr,ibeg_min,j,nsin,lenh,num,njm1k,nxk,ixk,nb,iyk,index, | iindex,kk,len,llen,nnxk,nnnxk,kkk,nphik,iunit,nsot,ios2,jamie, | nsut integer,external :: unlink EQUIVALENCE (INPT(106),STOP) real :: dum,tgrad real,external :: vsum real :: summary(100) ! ! Cray file routines for SGI (-L/usr/local/lib64/r8i4 -lncaru) integer,external :: readirec,readfrec integer :: nwds C **** READ CONTROL PARAMETERS CALL INPUT C **** DEFINE MODEL CONSTANTS CALL CON ! ! ixtimep is 4th dimension index to fg-array for previous time step ! ixtimec is 4th dimension index to fg-array for current time step ! (see fogcm.f) ! ixtimep = 1 ixtimec = 1 ! C **** C **** ZERO PERIODIC POINTS IN F ARRAY AND WORKING ARRAYS C **** LIM = 15*KMAXP1+9 DO 55 I = 1,2 DO 55 K = 1,LIM S15(I,K) = 0.0 S15(I+IMAXP2,K) = 0.0 55 CONTINUE LIM = 8*NCBUF+NCPHYS DO 57 I = 1,IMAXP4 DO 57 K = 1,LIM F(I,K) = 0.0 57 CONTINUE C **** CALCULATE DYNAMO CONSTANTS CALL CNSTNT C **** SET TAPE CONTROL PARAMETERS ITAP=1 IFIL=0 IHIS=0 ISAV=MSAV if (lusech.gt.0) then itapsech = 1 ifilsech = 0 ihissech = isecstart-iter isavsech = msavsech + (isecstart-iter) endif C **** ACCESS NEEDED TAPES CALL ACCSTP C **** C **** OPEN ALL DATA SETS C **** IOS1 = 0 IF (MODEST.EQ.1)THEN OPEN(UNIT=ISORC,FILE='SOURCE',STATUS='OLD',IOSTAT=IOS1,ERR=1234, 1 FORM='UNFORMATTED') 1234 CONTINUE istat = unlink('HISTCPY') OPEN(UNIT=IHIST,FILE='HISTCPY',STATUS='NEW',IOSTAT=IOS3, 1 ERR=3456,FORM='UNFORMATTED') 3456 CONTINUE call mkvolinfo(-1,-1,-1,-1, -1,-1,-1,-1, 0,0,tgcm_version, + 0,volinfo) ELSE OPEN(UNIT=IHIST,FILE='HISTCPY',STATUS='OLD',IOSTAT=IOS3, 1 ERR=7890,FORM='UNFORMATTED') 7890 CONTINUE ENDIF OPEN(UNIT=IMAG,FILE='MAGNET',STATUS='OLD',IOSTAT=IOS2,ERR=2345, 1 FORM='UNFORMATTED') 2345 CONTINUE IOS4 = 0 ! OPEN(8,FORM='UNFORMATTED') IOS5 = 0 ! OPEN(9,FORM='UNFORMATTED') C IOS6 = 0 IOS7 = 0 IF(IAMIE.EQ.1)THEN OPEN(UNIT=IUAMI,FILE='AMIE',STATUS='OLD',IOSTAT=IOS7,ERR=8765, 1 FORM='UNFORMATTED') 8765 CONTINUE ENDIF c c Open unit for secondary history volumes if lusech > 0: c if (lusech.gt.0) then istat = unlink('SECHIST') open(unit=lusech,file='SECHIST',status='NEW',iostat=ios9, + form='UNFORMATTED') if (ios9.eq.0) then write(6,"('start: opened unit ',i2,' for secondary ', + 'histories')") lusech else write(6,"('>>> start: error opening SECHIST', + ' (status=NEW) lusech=',i2,' ios9=',i4)") lusech,ios9 endif call mkvolinfo(-1,-1,-1,-1, -1,-1,-1,-1, 0,0,tgcm_version, + 1,volinfo_sech) c c Open direct access unit for non-primary secondary history fields: c (note record length is in bytes. See addfsech.f for write and mkfsech.f c for read of this file) c istat = unlink('SECHIST1') open(lusech1,file="SECHIST1",form='UNFORMATTED',access="DIRECT", + recl=zimxp*zkmxp*8,iostat=ios) ! recl=76*45*8=3420*8=15200 if (ios.eq.0) then write(6,"('Opened direct access unit ',i2,' for diagnostic', + ' secondary history fields')") lusech1 else write(6,"('>>> INPUT WARNING: error opening diagnostic', + ' secondary history field unit ',i2)") lusech1 endif endif IF(IOS1+IOS2+IOS3+IOS4+IOS5+IOS6+IOS7.NE.0)THEN WRITE(6,8901)IOS1,IOS2,IOS3,IOS4,IOS5,IOS6,IOS7 8901 FORMAT(* START --- PROBLEM OPENING DATA FILES, IOS=*,7I5) STOP ENDIF C **** CHECK AMIE VOLUME INFORMATION IF (IAMIE .EQ. 1) THEN C **** IYDA1 = YYDDD YEAR AND DAYNUMBER OF THE FIRST AMIE READ C **** IUTS1,2 = FIRST AND LAST UT SEC OF AMIE READS C **** NSPSEC = SPACING IN SECONDS BETWEEN AMIE READS C **** NTIMS = NUMBER OF TIMES WITH AMIE PATTERNS C **** JAMIE = NUMBER OF LATITUDE SLICES WITH AMIE RESULTS C **** (SHOULD BE SAME AS TIGCM, I.E. 36) READ (IUAMI) IYDA1, IUTS1, IUTS2, NSPSEC, NTIMS, JAMIE NSECT = (NTIMS-1) * NSPSEC + MOD(IUTS1,86400) NDAYAD = NSECT / 86400 IYRA = IYDA1 / 1000 IDA1 = IYDA1 - IYRA*1000 IDA2 = IDA1 + NDAYAD IYR4 = IYRA/4 IYR100 = IYRA/100 LPYR = 0 IF (IYR4*4 .EQ. IYRA .AND. IYR100*100 .NE. IYRA) LPYR=1 IENDA = 365 + LPYR IF (IDA2 .GT. IENDA) THEN IYRA = IYRA + 1 IDA2 = IDA2 - IENDA ENDIF ! FOR PAST YEAR'S END IYDA2 = IYRA*1000 + IDA2 WRITE (6,"(1X,' READ AMIE VOLUME FOR IYDA1,2 IUTS1,2 NSPSEC ', 1 'NTIMS JAMIE =', 7I8)") IYDA1,IYDA2,IUTS1,IUTS2,NSPSEC,NTIMS, 2 JAMIE IYD1 = (IYEAR-1900)*1000 + IIDAY IUTSMOD = ITER*C(4) SECTGCM = IUTSMOD IUTSMEND = (STOP(1,1)*24*60 + STOP(2,1)*60 + STOP(3,1))*60 WRITE (6,"(1X,' TIGCM RUN BEGINS AT IYD IUTS =', 2I8,' AND ', 1 'ENDS AT IUTS =',I8)") IYD1,IUTSMOD,IUTSMEND IF (IUTSMOD.LT.IUTS1-NSPSEC .OR. IUTSMEND.GT.IUTS2+NSPSEC) THEN WRITE (6,"(1X,' TIGCM begin IUTS or end IUTS is out of ', | 'range (+/-NSPSEC) of the AMIE volume --- STOP.')") STOP ENDIF IF (ZIMX .NE. 72 .OR. ZJMX .NE. JAMIE) THEN WRITE (6,"(1X,'DIMENSIONS OF IMX AND JMX OF TIGCM ARE ', 1 'INCOMPATIBLE WITH THE DIMENSIONS OF 73 BY',I3,' OF ', 2 'AMIE --- STOP.')") JAMIE STOP ENDIF C **** READ UP TO START TIME NRD = MAX0(1,(IUTSMOD-IUTS1)/NSPSEC + 1) NRD = MIN0(NTIMS-1,NRD) DO 60 N=1,NRD READ (IUAMI) SECUTA(1),HPSHA(1),HPNHA(1),CPSHA(1),CPNHA(1), 1 CUSPSTA(1),CUSPNTA(1),CUSPSLA(1),CUSPNLA(1),POTKVA1,VIXYZA1, 2 EKEVA1,EFLXA1 60 CONTINUE C **** READ ONE MORE TIME ON AMIE VOLUME NTST = NRD + 1 READ (IUAMI) SECUTA(2),HPSHA(2),HPNHA(2),CPSHA(2),CPNHA(2), 1 CUSPSTA(2),CUSPNTA(2),CUSPSLA(2),CUSPNLA(2),POTKVA2,VIXYZA2, 2 EKEVA2,EFLXA2 WRITE (6,"(1X,'READ NRD=',I4,' TIMES FROM AMIE VOLUME. LAST ', 1 'IS IN NO 1. TIMES IN STORAGE ARE:',2F9.0)") NRD,SECUTA(1), 2 SECUTA(2) WRITE (6,"(1X,'TIME=1: P CP MLT LAT (S,N) =', 4F7.2,2X,4F7.2)") 1 HPSHA(1),CPSHA(1),CUSPSTA(1),CUSPSLA(1),HPNHA(1),CPNHA(1), 2 CUSPNTA(1),CUSPNLA(1) WRITE (6,"(1X,'TIME=2: P CP MLT LAT (S,N) =', 4F7.2,2X,4F7.2)") 1 HPSHA(2),CPSHA(2),CUSPSTA(2),CUSPSLA(2),HPNHA(2),CPNHA(2), 2 CUSPNTA(2),CUSPNLA(2) ENDIF ! FOR IAMIE=1 C **** HISTORY TAPE START - SEARCH FOR FILE 70 IDATE(1)=OR(DATE(1.),0) ISTRT=ISORC IF(MODEST.LT.0) ISTRT=IHIST IF(MODEST.EQ.-1) THEN INPT(147)=INPT(103) INPT(148)=INPT(104) INPT(149)=INPT(105) END IF ! ! Search for source history: ! isearch = 1 75 continue ! ! Read header: nwds = readirec(istrt,itert,512,3,0) if (nwds /= 512) then write(6,"('>>> START: ERROR reading history header from', + ' unit ',i2,' nwds=',i5)") istrt,nwds stop 'i/o' endif ! ! Read summary: nwds = readfrec(istrt,summary,100,2,0) if (nwds /= 100) then write(6,"('>>> Error reading summary: istrt=',i3,' nwds=',i6)") + istrt,nwds stop 'i/o' endif ! ! Report time read: write(6,"('Reading time ',i3,':',i2,':',i2,' Seeking time ', + i3,':',i2,':',i2)") itert(2:4),(inpt(n+146),n=1,3) c c If searching OUTPUT volume, construct volinfo each time history is read c (use input DATE for iyd rather than header just read). c volinfo is written to mss file comment field. c if (modest.lt.0) then iend_iyd = (inpt(41)-1900)*1000+inpt(42) iend_day = itert(2) iend_hr = itert(3) iend_min = itert(4) idelmin = 0 iendonly = 0 if (isearch.eq.1) then ! first history found ibeg_iyd = (inpt(41)-1900)*1000+inpt(42) ibeg_day = itert(2) ibeg_hr = itert(3) ibeg_min = itert(4) else ! not first history found idelmin = (iend_hr-ibeg_hr)*60 + (iend_min-ibeg_min) iendonly = 1 endif call mkvolinfo(ibeg_iyd,ibeg_day,ibeg_hr,ibeg_min, + iend_iyd,iend_day,iend_hr,iend_min, + idelmin,iendonly,tgcm_version,0,volinfo) endif IFIL=IFIL+1 IF(IFIL.NE.1) GO TO 116 IDAY(1) = ITERT(2) IHR(1) = ITERT(3) 116 CONTINUE IF (ITERT(2).EQ.INPT(147).AND.ITERT(3).EQ.INPT(148).AND.ITERT(4) 1.EQ.INPT(149)) GO TO 140 ! ! Read past unwanted history: do j=1,jmax nwds = readfrec(istrt,dum,1,2,0) enddo isearch = isearch+1 GO TO 75 C **** LOAD DISK FROM HISTORY TAPE 140 continue ! ! History found -- read into global shared fg-array: ! (retain moulo, init, and flip temporarilly, for testing, ! see also putf call below) ! MODULO=3 CALL INIT do 240 j=1,jmax call flip ! keep flipping f-array temporarilly for testing nwds = readfrec(istrt,fg(1,1,j,ixtimep),lbscm,2,0) if (nwds <= 0) then write(6,"('>>> start: error reading latitude slice j=', | i3)") j stop 'i/o' endif ! ! Change from "external" to "internal" format: nxk = ncols+1 ixk = lbtap do nb=1,ncols iyk = imaxp4 nxk = nxk-1 do i=1,imaxp2 fg(iyk,nxk,j,ixtimep) = fg(ixk,1,j,ixtimep) ixk = ixk-1 iyk = iyk-1 enddo enddo C **** ADD PERIODIC POINTS do nxk = 1,ncols-1 fg( 1,nxk,j,ixtimep) = fg(imax+1,nxk,j,ixtimep) fg( 2,nxk,j,ixtimep) = fg(imax+2,nxk,j,ixtimep) fg(imax+3,nxk,j,ixtimep) = fg( 3,nxk,j,ixtimep) fg(imax+4,nxk,j,ixtimep) = fg( 4,nxk,j,ixtimep) enddo ! ! Model tgcm13 accepted histories that did not have dynamo potential. ! In this case, the following conditional was true, and the data was ! necessarilly shifted to make room for the potential. Because now ! "all" histories currently used contain potential, this ability is ! removed from tgcm13mt. ! if (nwds < lbtap) then write(6,"(/'>>> start: nwds < lbtap: lenh=',i5,' lbtap=',i5)") | nwds,lbtap write(6,"(4x,'tgcm13mt does not accept histories without', | ' dynamo potential.',/,4x,'Please use tgcm13.')") stop 'start' endif ! ! Define dynpot: nphik = nphi do k = 1,kmaxp1 nphik = nphik+1 do i = 1,imax+1 dynpot(i,j,k) = fg(i+2,nphik,j,ixtimep) enddo enddo ! ! Copy fg(..j..) into f-array (at nj) so testing can be continued: call putf(nj,j,ixtimep,1,ntape) ! ! End j read loop: 240 CONTINUE ! ISAV=MSAV IHIS=MHIS ITAP=1 IF(IFIL.LT.MFIL)GO TO 330 IFIL=0 ITAP=2 ! REWIND ISTRT 330 CONTINUE ! ! Get geographic coords at ends of magnetic field lines: ! (subroutine othend defines /othend_com/ in othend.h) ! call othend ! ! ! Insert polar points in dynpot: ! do k = 1,kmaxp1 dynpot(1,0,k) = (9.*vsum(imaxg,dynpot(1,1,k),1)- | vsum(imaxg,dynpot(1,2,k),1))/(8.*float(imaxg)) dynpot(1,jmaxgp,k) = (9.*vsum(imaxg,dynpot(1,jmaxg,k),1) | -vsum(imaxg,dynpot(1,jmaxg-1,k),1))/(8.*float(imaxg)) do i = 2,imaxgp dynpot(i,0,k) = dynpot(1,0,k) dynpot(i,jmaxgp,k) = dynpot(1,jmaxgp,k) enddo enddo C **** C **** VALUE OF T AT LOWER BOUNDARY AND GRADIENT C **** TBOUND = TTBOUND TGRAD = 6. T0(1) = TBOUND T0(2) = TBOUND+C(3)*TGRAD CALL BNDRY CALL BNDRY2 CALL BNDRYA ! ! Init sigmas, segeuv, etc (see qrj.f): call init_sigmas ! C **** SET T0=0. DO 261 K=1,KMAX T0(K)=0. 261 CONTINUE T0(KMAXP1)=0. C **** READ MAGNETIC DATA FOR DYNAMO BUFFER IN(IMAG,1)(ALATM(1,0),RJACD(ZIMXP1,ZJMXP1)) IF(UNIT(IMAG))3,4,4 4 IUNIT = IMAG WRITE(6,100)IUNIT 100 FORMAT(* PROBLEM READING TGCM FIELDS FROM UNIT *,I5) CALL EXIT 3 CONTINUE BUFFER IN(IMAG,1)(IG(1,1),DJM(ZIMXP1,ZJMXP1)) IF(UNIT(IMAG))5,6,6 6 IUNIT = IMAG WRITE(6,100)IUNIT CALL EXIT 5 CONTINUE BUFFER IN(IMAG,1)(CSLATM(1,1),SNLONG(ZIMXP1)) IF(UNIT(IMAG))7,8,8 8 IUNIT = IMAG WRITE(6,100)IUNIT CALL EXIT 7 CONTINUE CALL MAGDYN C **** C **** TRANSFORM DYNPOT TO GEOMAGNETIC COORDINATES IN C **** PHIM3D(IMAXMP,JMAXM,-2:ZKMXP) C **** DO 343 K = 1,KMAXP1 DO 343 J = 1,JMAXM CALL GRDINT(PHIM3D(1,J,K),DYNPOT(1,0,K),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) 343 CONTINUE C **** C **** PERIODIC POINTS C **** DO 348 K = 1,KMAXP1 DO 348 J = 1,JMAXM PHIM3D(IMAXMP,J,K) = PHIM3D(1,J,K) 348 CONTINUE C **** C **** DELETE SOURCE AND MAGNETIC DATA FILES C **** IF(ISTRT.EQ.ISORC)THEN CLOSE(ISORC) ISTAT = UNLINK('SOURCE') IF(ISTAT.NE.0)THEN WRITE(6,9012)ISTAT 9012 FORMAT(* START --- PROBLEM DELETING SOURCE FILE, ISTAT=*,I2) STOP ENDIF ENDIF C **** C **** DELETE MAGNETIC DATA FILE C **** CLOSE(IMAG) ISTAT = UNLINK('MAGNET') IF(ISTAT.NE.0)THEN WRITE(6,9876)ISTAT 9876 FORMAT(* START --- PROBLEM DELETING MAGNETIC FILE, ISTAT=*,I2) STOP ENDIF IF(MODEST.LT.0) RETURN C **** INITIALIZE FOR NEW HISTORY TAPE IFIL=0 IHIS=0 ITAP=1 RETURN 4000 WRITE(6,4010) NUM,LEN,NSOT,J 4010 FORMAT(*1DISK PROBLEM IN START AT NUM=*I5,2X,*LEN=*I5,2X,*NSOT=* AI1,2X,*J=*I2) CALL EXIT 5000 WRITE(6,5010) NUM,LEN,NSIN,J 5010 FORMAT(*1HIST PROBLEM IN START AT NUM=*I5,2X,*LEN=*I5,2X,*NSIN=* AI1,2X,*J=*I2) CALL EXIT 6000 WRITE(6,6010)NUM,LEN,NSUT,J 6010 FORMAT(*1DISK PROBLEM IN START. NUM=*I5,2X,*LEN=*I5,2X,*NSUT=*I5, A2X,*J=*I2) CALL EXIT 1020 FORMAT(//* HISTORY TAPE LABELS*/(10A12)) 1030 FORMAT(//* SAVE TAPE LABELS*/(10A12)) END C