
      SUBROUTINE BASIN
C
C     BASIC PARAMETERS DERIVING
C
      INCLUDE 'basp.h'
      INCLUDE 'garr.h'
      INCLUDE 'din.h'
      INCLUDE 'wrm.h'  ! NGZ

      CHARACTER*8 ATE(3)
      DIMENSION PAG(46)         ! smotri DTG and PAG and ATE and PABK

      DP=0.2844   ! ALTITUDE STEP (PART OF H)
      RA=0.1      ! BACKGROUND AEROSOL RADIUS (um)
      ITM=5000    ! MAX NUMBER OF ITERATIONS
Cgal      ITM=50    ! MAX NUMBER OF ITERATIONS
      EKZ=2.E3    ! BACKGROUND KZZ (KIDA,1983)
C     Variables IRD and IWR are not in use.
C     It's not possible to continue a previous run anymore. GM2/00
C
C     ***KEEP IRD=2 AND IWR=2 -OTHER OPTIONS DISABLED- GM1/00***
C      IRD=2       ! READNG OF INPUT CONCENTRATIONS WHEN CONTINUING
C                  ! A PREVIOUS RUN (1-FORMAT,2-UNFORMAT.)
C      IWR=2       ! WRITING IN CASE THE RUN IS CONTINUED AT A LATER TIME
C                  ! (1-FORMAT,2-UNFORMAT.)
      RELM=50.
      UGMA=97.5
C******************     MODEL PARAMETERS READING  ********************
C     DTF - TIME STEP FOR SHORT-LIVED GASES (S)
C     UGMA  MAX.ZENITH ANGLE FOR DAY
C     RELM- MAX VARIATION OF OPTICAL MASS
C     DYM - MAX REL.VARIATION
C     TPM - PERCENT FOR TIME STEP MISSING
C     DT  - TIME  STEP IN DIURNAL CALCULATIONS
C           (i.e., PHOTOLYSIS RATE CALC.) FOR SHORT-LIVED SPECIES (S)
C     NCH - TIME STEP NUMBER (FOR LONG-LIVED)
C     ITM - MAX.NUMBER OF ITERATIONS
C     NOEB- SWITCH FOR OFF-LINE RUNNING
C     FI  - LATITUDE
C
      OPEN(UNIT=16,FILE=FLN6,STATUS='OLD',ERR=30)
      GOTO 35
 30   WRITE(*,*) 'File NOT opened: ',FLN6
 35   CONTINUE

      READ(16,105)DTF,DTG,UGMA,G_LIG,DYM,NDLON,ISEM,ITI,NOEB,IPS,
     $IPSC1,IPSC2,NVA 
C    'NPE','IPS','ITD','NGOD_B','NVA','IWB','ITMD' and 'IKT' not used  GM2/00
C    'ISEM' introduced by GM3/01
      DYM=DYM*0.01
      TPM=0.01
C      Ti=iti
C      WRITE(*,'(A14,F6.1)') ' LATITUDE = ',fi  GM3/01

C    Next block is not used GM2/00
c--------------------------------------------------------------
C      KG=0
C 108  READ(16,112)ATE  !Species for printing. ATE(2) and ATE(3) are empty.
C
C      CALL COMKOM(AY,KFE,ATE(1),NGRR)
C      IF(NGRR.GT.NY)GOTO106
C      KG=KG+1
C      NGR(KG)=NGRR
C      GOTO108
C--------------------------------------------------------------
 112  FORMAT(4X,A8,3X,A8,1X,A8)

C     NOW READS THE 38 FIXED PRESSURE LEVELS IN FLN6
Cgal                                 error in PA, nado PAG
 106  READ(16,'(7F10.4)') PAG
      READ(16,'(10X)')
Cgal      DO i=1,L
Cgal        PABK(i)=PA(i)   !Backup of the pressure levels read  GM2/00
Cgal      ENDDO

      NZAP=0
 102  READ(16,112)ATE   !Species for writing in FTA3
      CALL COMKOM(AY,KFE,ATE(1),NZRR)    !COMKOM IS IN COMKOM.FOR
      IF(NZRR.GT.NY)GOTO101
      NZAP=NZAP+1
      NGZ(NZAP)=NZRR
      GOTO102

 101  CLOSE(16)
 105  FORMAT(4(50X,G8.2/),50X,G8.2/,8(50X,I8/))
 100  RETURN
      END
C
C ************************************************************
C
      SUBROUTINE PH_VDV
      character*64 AAA
      character*32 react
      INCLUDE 'basp.h'
      INCLUDE 'din.h'
      INCLUDE 'pharr.h'
      INCLUDE 'scfl.h'
      DO 1 I=1,L
         DO 1 J=1,MNabs
            CONC(I,J)=0.
            DO 1 K=1,Len_tm
   1           SIG(J,I,K)=0.
      do ir = 1, Ngas
         K_T(ir) = 1
      enddO
      NL = MNL
      NABS = MNABS
      Len_t = Len_tm
c     Read input parameters and photolysis reactions:
      OPEN ( UNIT=25,FILE=FLN3,STATUS='OLD',ERR=30)
      GOTO 35
 30   WRITE(*,*) 'File NOT opened: ',FLN3
 35   CONTINUE
      read ( 25,217 ) N_O3, N_O2, ALB, aerosol
 217  format ( 2(45x,I5/), 40x, f10.3/, 45x, i5/ )
      irr = 1
      do ir = 1, 200
         read ( 25, 517, ERR=518) Num_r, react
         if ( Num_r .ne. 0 ) then
           GAS(irr) = react
           irr = irr+1
         endif
         enddo
  518 continue
      CLOSE(25)
  517 format( i3, 2x, A32 )
      A_gas(1) = GAS(N_O2)
      A_gas(2) = 'O3tot   '
      n_abs(1) = N_O2
      n_abs(2) = N_O3
 520  format (2x, A8 )
c **********************************************************************
c    read NEW wavelength grid and TOA fluxes  8/17/97
      OPEN (UNIT=22, FILE=FLN7, STATUS='OLD',ERR=40)
      GOTO 45
 40   WRITE(*,*) 'File NOT opened: ',FLN7
 45   CONTINUE

      Nheader = 24
      read (22,'(A64)') ( AAA, i=1, Nheader )
      do jg=1,NGLIS
         read (22,'(A64)')AAA
C         print *,AAA ; pause
         read (22,'(2f8.3,12e12.4)')(wbd(i),we(i),f88(i,:,jg),
     $   i=1,Len_tm)
      enddo
      do i=1, Len_tm
         wlt(i) = .5*(wbd(i)+we(i))
      enddo

      CLOSE (22)
      OPEN (UNIT=57,FILE=FTA7,STATUS='OLD')
      OPEN (UNIT=58,FILE=FTA1,STATUS='OLD')
C         READ(57,'(A64)')(AAA,i=1,1)
C         print *,AAA ; pause
         DO JG=1,NGROZ
            DO JM=1,12
               READ(57,'(i4,i3,14E12.4)')IYER,IMON,
     $         RADLISMON(1:13,JM,JG),RADLISMON(16,JM,JG)
               READ(58,'(7x,2E12.4)')RADLISMON(14:15,JM,JG)
C               print *,JM,RADLISMON(:,JM,JG)
            ENDDO
CSOL            READ(57,'(A64)')AAA
C            print *,IYER,IMON ; pause
         ENDDO
      CLOSE(57)
      CLOSE(58)
c
c ================================================================
c  Aerosol parameters:
      if ( aerosol .eq. 0 ) then
         do ih = 1, L
            do il = 1, MNL
               ext(ih,il) = 0.
               ssam(ih,il) = 1.
               gm(ih,il) = 0.
            enddo
         enddo
      endif
      RETURN
      END
C
C ************************************************************
C
      SUBROUTINE NAM
C
C     GASES NUMBERING
C
      INCLUDE 'basp.h'
      INCLUDE 'garr.h' ! NC,NY,KFE,NYPR,NOE,NBN,KEN,SER,WD,NCA,NDA,NCD
      INCLUDE 'fam.h'  ! MS,NBS,NFS,NUS
      INCLUDE 'din.h'  ! L
      INCLUDE 'wrm.h'  ! NGZ
      CHARACTER SPPR*1, SPS*1, SPD*1
Cgal dobavka
      COMMON /CHMDOL/ NUMDOL(NC),AYDOL(NC),
     *                NBNDOL(NC),SERDOL(NC),WDDOL(NC),AMDOL(NC)
      CHARACTER*8 AYDOL
      COMMON /NCALDL/ NCAL,NDAL,NCDL
C
C     GASES NAMES READING
C
      OPEN(UNIT=11,FILE=FLN1,status='OLD',ERR=40)
      GOTO 45
 40   WRITE(*,*) 'File NOT opened: ',FLN1
 45   CONTINUE
      DO 203 I=1,MS
 203     NBS(I)=0
      NCAL=0
      NDAL=0
      NCDL=0
      ig=0                   ! Cgal dobavka
      DO 201 I=1,NC
         II=I
         KEN(I)=L
         READ(11,'(A1,4X,A8,I3,2(1X,A1),I3,2E10.2,1X,I3,F6.2)') SPPR,
     $   AY(I),NOE(I),SPD,SPS,NBN(I),SER(I),WD(I),NZIM(I),AMOL(I)
         if (NZIM(I).ne.0) then
            ig=ig+1
            NUMDOL(ig)=I
            AYDOL(ig) = AY(I)
            NBNDOL(ig) = NBN(I)
            SERDOL(ig) = SER(I)
            WDDOL(ig) =  WD(I)
            AMDOL(ig) =  AMOL(I)
         endif
         IF(SPPR.NE.'*')THEN
            NYPR(I)=4
            IF(NBN(I).LE.L)THEN           ! all calculated species
               NDAL=NDAL+1
               NDA(NDAL)=I
               NYPR(I)=1
            ENDIF
            IF(SPS.EQ.'S')THEN            ! short-lived species
               NCAL=NCAL+1
               NCA(NCAL)=I
               IF(SPD.EQ.'D')THEN         ! dynamical short-lived species
                  NYPR(I)=3
               ELSE
                  NYPR(I)=2
               ENDIF
            ENDIF
            IF(SPD.NE.'D')THEN            ! non dynamical species
               NCDL=NCDL+1
               NCD(NCDL)=I
            ENDIF
         ELSE
            NYPR(I)=0
         ENDIF
C
C     SS3/00  Only short-lived species make up families (semr is in this file)
C
         IF(SPS.NE.'S')GOTO201
 199     CALL SEMR(AY(I),II,NBS(2),NFS(1,2),NUS(1,2),NY,'N',1)
         CALL SEMR(AY(I),II,NBS(3),NFS(1,3),NUS(1,3),NY,'CL',2)
         CALL SEMR(AY(I),II,NBS(4),NFS(1,4),NUS(1,4),NY,'BR',2)
         CALL SEMR(AY(I),II,NBS(5),NFS(1,5),NUS(1,5),NY,'S',1)
 201     CONTINUE
Cgal      PRINT *,' SPECIES WITH DIURNAL CALCULATIONS, NUMBER ',NCAL
 177  FORMAT(8(2X,A8))
      NBS(1)=3
      CALL COMKOM(AY,KFE,'O3      ',NFS(1,1))   !COMKOM IS IN COMKOM.FOR
      CALL COMKOM(AY,KFE,'O       ',NFS(2,1))
      CALL COMKOM(AY,KFE,'O(1D)   ',NFS(3,1))
      DO 204 J=1,3
 204     NUS(J,1)=1
C
C     GASES NUMBERS DERIVING
C
      CALL COMKOM(AY,KFE,'O3      ',KO3)
      CALL COMKOM(AY,KFE,'NO2     ',NO2)
      CALL COMKOM(AY,KFE,'NO3     ',NNO3)
      CALL COMKOM(AY,KFE,'N2O5    ',NN2O5)
      CALL COMKOM(AY,KFE,'CLO     ',NCLO)
      CALL COMKOM(AY,KFE,'BRO     ',NBRO)
      CALL COMKOM(AY,KFE,'H2O     ',KH2O)
      CALL COMKOM(AY,KFE,'NOY     ',NNOY)
      CALL COMKOM(AY,KFE,'CLY     ',NCLY)
      CALL COMKOM(AY,KFE,'BRY     ',NBRY)
      CALL COMKOM(AY,KFE,'CO      ',NCO)
      CALL COMKOM(AY,KFE,'CH4     ',NCH4)
      CALL COMKOM(AY,KFE,'HNO3    ',NHNO3)
      CALL COMKOM(AY,KFE,'HCL     ',NHCL)
      CALL COMKOM(AY,KFE,'N2O     ',NN2O)
      CALL COMKOM(AY,KFE,'OCLO    ',NOCLO)
      CALL COMKOM(AY,KFE,'CLONO   ',NCLONO)
      CALL COMKOM(AY,KFE,'SOX     ',NSOX)
      CALL COMKOM(AY,KFE,'SO2     ',NSO2)
      CALL COMKOM(AY,KFE,'H2SO4   ',NH2SO4)
      CALL COMKOM(AY,KFE,'NO      ',NGE(2))
      CALL COMKOM(AY,KFE,'CL      ',NGE(3))
      CALL COMKOM(AY,KFE,'BR      ',NGE(4))
      CALL COMKOM(AY,KFE,'HBR     ',NGE(5))
      CALL COMKOM(AY,KFE,'HCL     ',NGE(6))
      CALL COMKOM(AY,KFE,'N       ',NGE(7))
      CALL COMKOM(AY,KFE,'CO2     ',NCO2)
      CALL COMKOM(AY,KFE,'O       ',NO)
      CLOSE(11)
      NGE(1)=KO3
      NGE(8)=NCLONO
      NGEN=8
      RETURN
      END

      SUBROUTINE PHOT_R(NCH,FLN,AY,KFE,KS,KSS,KSC,NS,NSC)
C
C     PHOTODISSOCIATION PROCESSES READING
C
      CHARACTER*60 FLN
      CHARACTER*8 ATE(3), AY(KFE)
      DIMENSION KS(1,NS),KSS(2,NSC),KSC(1,NSC)
      OPEN(UNIT=NCH,FILE=FLN,status='OLD',ERR=40)
      GOTO 45
 40   WRITE(*,*) 'File NOT opened: ',FLN
 45   CONTINUE
      READ(NCH,121)
      NCS=0
      DO 2 I=1,NSC
         READ(NCH,113)ATE
         CALL COMKOM(AY,KFE,ATE,KSC(1,I))  !COMKOM IS IN COMKOM.FOR
         IF(KSC(1,I).GT.KFE)GOTO28
         NCS=NCS+1
         KS(1,NCS)=KSC(1,I)
  28     DO 22 J=2,3
  22        CALL COMKOM(AY,KFE,ATE(J),KSS(J-1,I))
   2     continue
      CLOSE(NCH)
 121  FORMAT(4(45X/))
 113  FORMAT(5X,A8,3X,A8,1X,A8)
      RETURN
      END

      SUBROUTINE PARAM_D
C
C     PRESSURE (PA) AND LOG-ISOBARIC ALTITUDE CALCULATION
C
      INCLUDE 'basp.h'
      INCLUDE 'din.h'
      REAL SIGK(L)
      DATA PI/3.141592/, H/7.E5/, RZ/6.367E8/

      DATA SIGK /3.573e-008,7.267e-008,
     $     1.478e-007, 3.007e-007, 6.115e-007, 1.244e-006,
     *     0.00000253, 0.00000445, 0.00000780, 0.00001370,
     *     0.00002400, 0.00004210, 0.00007390, 0.00012400,
     *     0.00020300, 0.00031800, 0.00048600, 0.00073000,
     *     0.00110000, 0.00165000, 0.00248000, 0.00372000,
     *     0.00560000, 0.00842000, 0.01270000, 0.01900000,
     *     0.02860000, 0.04290000, 0.06460000, 0.09000000,
     *     0.12000000, 0.17000000, 0.23000000, 0.30000000,
     *     0.38000000, 0.46000000, 0.54000000, 0.62000000,
     *     0.70000000, 0.78000000, 0.85000000, 0.91000000,
     *     0.95000000, 0.97500000, 0.99300000 /
      DTR=PI/180.
      DO 501 I=1,L
         PA(I)=SIGK(L-I+1)*1000.
         Z(I)=H*ALOG(1000./PA(I))
         ZBK(I)=Z(I)      !Backup GM2/00
         ZKM(I)=Z(I)*1.E-5
Cgal         PRINT *,ZKM(I),PA(I)
 501     ZKMBK(I)=ZKM(I)  !Backup GM2/00
      ZST(1)=0.
      ZST(LL)=Z(L)
      DO 502 I=2,L
 502     ZST(I)=0.5*(Z(I)+Z(I-1))
      FCY(0)=-2.*DTR*RZ
      DO 510 K=1,M
         FI(K)=-92.+4.*K               ! *** Latitude (degree)
         FCY(K)=(FI(K)+90.)*DTR*RZ     ! *** Distance from SP (cm)
         NRM(K)=NVA
         IF(K.NE.1)CSI(K)=COS(0.5*(FI(K)+FI(K-1))*DTR) ! Cos at midpoints
 510     CS(K)=COS(FI(K)*DTR)          ! *** Cosine at gridpoints
      CSI(1)=0.                        ! *** Cosine at SP
      CSI(MMA)=0.                      ! *** Cosine at NP
      FST(1)=0.
      FCY(MMA)=182.*DTR*RZ
      DO 503 K=2,M
 503     FST(K)=(0.5*(FI(K)+FI(K-1))+90.)*DTR*RZ ! Distance at midpoints
      FST(MMA)=180.*DTR*RZ             ! *** Distance from SP to NP
      DO 504 J=1,NLO
         ALON(J)=-185.+5.*J
         JNLOH(J)=J+NLO/2 !+1            ! *** Opposite at the latitude circle
         IF(JNLOH(J).GT.NLO)JNLOH(J)=JNLOH(J)-NLO
         ALYM(J)=(-5.+5.*J)*DTR*RZ
 504     ALYMS(J)=ALYM(J)+2.5*DTR*RZ
      !do j=1,NLO
      !   print *,J,ALON(J),JNLOH(J),ALON(JNLOH(J))
      !enddo
      !stop
      JNLOH(0)=JNLOH(NLO)
      JNLOH(NLO1)=JNLOH(1)
      ALON(0)=-185.
      ALON(NLO1)=180.
      ALYMS(NLO1)=ALYMS(NLO)+5.*DTR*RZ
      ALYM(0)=-5.*DTR*RZ
      ALYM(NLO1)=360.*DTR*RZ
C
C     PARAMETER FOR DIFFUSIVE CALCULATION
C
  33  DO 3 I=1,L
         DO 3 K=1,M
            DO 3 J=1,NLO
               NCOU(I,K,J)=1
               IF(I.EQ.1)NCOU(I,K,J)=2
               IF(I.EQ.L)NCOU(I,K,J)=3
 3             CONTINUE
      RETURN
      END

C234567890123456789012345678901234567890123456789012345678901234567890
      SUBROUTINE INIT_R        ! 'NCDL' and 'NDAL' not used GM 1/00
      INCLUDE 'PRMT.INC'
      INCLUDE 'basp.h'
      INCLUDE 'garr.h' ! NC,NY,SER,NYPR,NBN,NCO,NAR
      INCLUDE 'din.h'  ! FI,PA,AM,NREL,DURO,OPMDO,NRM,EM,GCR
      INCLUDE 'fam.h'  ! MS,NBS,NFS,NUS
      INCLUDE 'yarr.h' ! Y,Y3,P,D
      INCLUDE 'wrm.h'
      INCLUDE 'trajarr.h'
      INCLUDE 'pscarr.h'
Cgal                            smotri TA and L,M,12 in TA ! otladka
      DIMENSION BNL(13,2),ABL(13),EMIS(13,2),Y2D(NC,L,M)
      INCLUDE 'xx.h'
      COMMON /HISDAY/ TXYZ(NLONP2,NLAT,NLEV),T_DAY(NLONP2,NLAT,NLEV)
      DATA YMIN/1.E-20/
      DATA BNL/2.,5.,10.,50.,100.,120.,100.,120.,200.,300.,200.,50.,15.,
     *70.,72.,74.,76.,78.,80.,90.,100.,120.,125.,120.,100.,90./
      DATA EMIS/0.,0., 0., 0.1,0.4,0.7,1.4,1.,1.7,4.0,0.9,0.1,0.,
     *          0.,0.1,0.2,0.5,1.1,2.8,2.4,3.,2.8,4.1,2.1,0.2,0./
      DATA ABL/-90.,-75.,-60.,-45.,-30.,-15.,0.,15.,30.,45.,60.,75.,90./
      PPA=1./1.3806E-19
      TI=1.
      DO 1 K=1,M
         DO 1 J=1,2
   1        BONL(K,J)=A18KWI(FI(K),ABL,BNL(1,J),1,13)
C
C     INITIAL NUMBER DENSITIES READING
C
  31  OPEN(UNIT=21,FILE=FLN5,STATUS='OLD',ERR=40)
      GOTO 45
 40   WRITE(*,*) 'File NOT opened: ',FLN5
 45   CONTINUE
      READ(21,703)Y2D(:,1:39,:)
      CLOSE(21)
      DO I=40,L
         DO J=1,NY
            Y2D(J,I,:)=Y2D(J,39,:)*0.489542**(I-39.)
         ENDDO
         Y2D(KO3,I,:)=Y2D(KO3,39,:)/Y2D(NC,39,:)*0.489542**(I-39.)*
     $   Y2D(NC,I,:)
      ENDDO
      DO 48 K=2,NLAT
         KK=K-1
         DO 48 ILON=1,NLONP2
            DO 48 I=1,L
               DO 49 J=1,NC
 49               Y(J,I)=AMAX1(Y2D(J,I,KK),YMIN)
C
C              Makes families equal to the sum of its component species 
C              before the run starts.
C     
               DO J=1,MS
                  YY=0.
                  DO JI=1,NBS(J)
                     YY=YY+Y(NFS(JI,J),I)*NUS(JI,J)
                  ENDDO
                  Y(J,I)=YY
               ENDDO
               SPSC1(I,ILON,K)=YMIN
               SPSC2(I,ILON,K)=YMIN
               SAR(I,ILON,K)=YMIN
               SHNO3(I,ILON,K)=YMIN
               SH2O(I,ILON,K)=YMIN
               SH2SO4(I,ILON,K)=YMIN
               IF(IPS.EQ.0)THEN
                  SHNO3(I,ILON,K)=XXALLY(ILON,K,I,NC-1)
                  SH2O(I,ILON,K)=XXALLY(ILON,K,I,NC-2)
	       ENDIF
               AM(I)=PA(I)/TXYZ(ILON,K,NLEV-I+1)*PPA
               Y(NC,I)=AM(I)
               Y(NC-1,I)=AM(I)*0.209476
               Y(NC-2,I)=AM(I)*0.78
               DO 318 J=1,NY
                  Y3(J,I)=Y(J,I)
 318              CONTINUE
!               print *,IPS,XXALLY(ILON,K,I,31)
               DO 48 J=1,NC
                  XXALLY2(ILON,K,I,J)=1.
                  IF(IPS.EQ.0)GOTO48
                  XXALLY(ILON,K,I,J)=AMAX1(Y(J,I),YMIN)
                  if(i.GT.l-10)XXALLY(ILON,K,I,J)=XXALLY(ILON,K,L-10,J)/
     $            XXALLY(ILON,K,l-10,NC)*XXALLY(ILON,K,I,NC)
  48              CONTINUE
!      pause
      DO 314 K=1,M
         NREL(K)=NRM(K)
         DURO(K)=86400.
         OPMDO(K)=1.
 314     CONTINUE
C      call minmax3(XXALLY(:,:,:,3)/XXALLY(:,:,:,NC),NLONP2,NLAT,NLEV,
C     $Ymi,Yma)
C       pause
C
c*******No continuation of a previous run is allowed any longer. GM2/00*******
C---------------------------------------------------------------------------
C      PRINT *,' NEW RUN: READING FROM ',FLN5
C      GOTO33
C
C     READING AFTER CALCULATION
C
C     IN CASE THE USER CHOOSES TO CONTINUE A PREVIOUS RUN
C
c  32  CONTINUE
c      GOTO(304,305),IRD  !If IRD=1 it goes to 304, if IRD=2 it goes to 305
c
C     THIS FORMATTED READING (WHEN CONTINUING A PREVIOUS RUN) IS DISABLED GM1/00
C-------------------------------------------------------------------------------
C     FORMATTED READING
C
C 304  OPEN(UNIT=21,FILE=FLN9,STATUS='OLD',ERR=50)
C      WRITE(*,*) 'File opened: ',FLN9
C      GOTO 55
C 50   WRITE(*,*) 'File NOT opened: ',FLN9
C 55   CONTINUE
C      READ(21,'(F7.2,3I6)')TI,NGOD,NCHM,NREL
C      READ(21,705,ERR=37)Y,Y3,HNO3S,H2OS,H2SO4
C  37  READ(21,705)P,D
C      CLOSE(21)
C      PRINT *,' READING FORMATTED '
C------------------------------------------------------------------------------
C
C    Also the unformatted reading below has been disabled. No continuation
C    of a previous run is allowed any longer. GM2/00
C------------------------------------------------------------------------
C     UNFORMATTED READING FOR CONTINUING A PREVIOUS RUN
C
c 305  OPEN(UNIT=21,FILE=FLN9,FORM='UNFORMATTED',STATUS='OLD',ERR=60)
c      GOTO 65
c 60   WRITE(*,*) 'File NOT opened: ',FLN9
c 65   CONTINUE
c         READ(21)TI,NGOD,NCHM,NREL
c         PRINT *,' CONTINUING RUN: READING (UNF.) FROM ',FLN9
c         READ(21)Y,Y3,HNO3S,H2OS,H2SO4S
c         READ(21)P,D
c      CLOSE(21)
c------------------------------------------------------------------------
  33  CONTINUE
 703  FORMAT(8(1PE10.4))
 705  FORMAT(6(1PE13.6))
      RETURN
      END

      SUBROUTINE GCR_R(NCH,FTA2,GCR,L,L1,M,ZKM,FI)
      CHARACTER*60 FTA2
      DIMENSION GCR(L,M),GCR74(37,36),ZKM37(37),ALAT(36),
     $ZKM(L),FI(M)
C
C     GALACTIC RAYS READING
C
      DO 1 I=1,37
 1       ZKM37(I)=2.*(I-1)
      DO 2 K=1,36
 2       ALAT(K)=-92.5+5.*K
      OPEN(UNIT=NCH,FILE=FTA2,STATUS='OLD',ERR=70)
      GOTO 75
 70   WRITE(*,*) 'File NOT opened: ',FTA2
 75   CONTINUE
      READ(NCH,'(161(1X/))')
      DO 3 K=19,36
         READ(NCH,'(1X)')
         READ(NCH,'(6e11.3)')(GCR74(I,K),I=1,37)
         READ(NCH,'(1X)')
         DO 3 I=1,37
 3          GCR74(I,37-K)=GCR74(I,K)
      CLOSE(NCH)
      print *,zkm
      print *,zkm37
      print *,FI
      print *,ALAT
      CALL LNT2(GCR74,ZKM37,ALAT,37,36,GCR,ZKM,FI,L,M,37,36)
      CALL PRNS(GCR,1,L,M,ZKM,FI,2,1,6,' GCR    ')
      END

      SUBROUTINE SEMY(Y,NC,MS,NB,NF,NU)
      DIMENSION Y(NC),NB(MS),NF(NC,MS),NU(NC,MS)
      DO 1 J=1,MS
         Y(J)=0.
         DO 1 I=1,NB(J)
   1        Y(J)=Y(J)+Y(NF(I,J))*NU(I,J)
      RETURN
      END

      SUBROUTINE SEMR(AY,I,NBS,NFS,NCS,ND,OB,NOB)
C     NBS - NUMBER GASES IN FAMILY
C     NFS - NO FAMILY MEMBER
C     NCS - PARAMETER TIMES IN FAMILY
      CHARACTER AY*8 !, OB*2  !OB definition changed by GM3/01
      CHARACTER*1 OB(NOB)
      DIMENSION NFS(ND),NCS(ND)
      NF=9-NOB
      DO 202 K=1,NF
         KF=K-1+NOB
         KFF=KF+1
         IF(AY(K:K).NE.OB(1))GOTO202  !from OB(1) to OB(1:1) GM3/01
         IF(NOB.EQ.2)THEN
            IF(AY(KF:KF).NE.OB(2))GOTO202 !from OB(2) to OB(2:2) GM3/01
         ENDIF
         NBS=NBS+1
         NFS(NBS)=I
         NCS(NBS)=1
         IF(K.NE.NF.AND.AY(KFF:KFF).EQ.'2')NCS(NBS)=2
 202     CONTINUE
      RETURN
      END

      SUBROUTINE LNT1(FIN,ZIN,NZIN,FOUT,ZOUT,NZOUT)
      DIMENSION FIN(NZIN),ZIN(NZIN),FOUT(NZOUT),ZOUT(NZOUT)
      DO 1 I=1,NZOUT
         II=1
         DO WHILE(ZIN(II)+0.1.LT.ZOUT(I))
            II=II+1
            IF(II.GT.NZIN)GOTO1
         ENDDO
         IF(II.EQ.1)II=II+1
         ZXP=ZIN(II)-ZOUT(I)
         ZXM=ZOUT(I)-ZIN(II-1)
         ZNZ=1./(ZIN(II)-ZIN(II-1))
         FOUT(I)=(FIN(II-1)*ZXP+FIN(II)*ZXM)*ZNZ
   1     CONTINUE
      RETURN
      END

      FUNCTION A18KWI(X,XN,YN,M,N)
      DIMENSION XN(N),YN(N)
      K=M
      L=N-1
 1    I=INT((K+L)/2.)
      IF(I.EQ.K) GOTO 4
      IF(X-XN(I)) 2,3,3
 2    L=I
      GOTO 1
 3    K=I
      GOTO 1
 4    CONTINUE
      Z=(X-XN(K+1))/(XN(K)-XN(K+2))
      A18KWI=YN(K)+((X-XN(K))/(XN(K)-XN(K+1)))*(1.+Z)*(YN(K)-YN(K+1))-
     *((X-XN(K))/(XN(K+1)-XN(K+2)))*Z*(YN(K+1)-YN(K+2))
      RETURN
      END
