      SUBROUTINE CSM2(TI,Z,UGM,UG,TWR,TZR)
      SUM=0.
      CALL TWZ(TI,Z,TW,TZ)
      OPMP=OPTM(UGM)
      CALL CSTU(TI,Z,OPMP,TZR,-1.)
      TWR=TW-TZR+TZ
      IF(TW.GT.TZ.AND.TZR.LT.0.)GOTO5
      T=TWR
      DTT=3600.
   2  T2=T+DTT
      IF(T2.LE.TZR)GOTO4
      DTT=TZR-T
   4  CALL CSM(TI,Z,U,T,DTT)
      T=T2
      SUM=SUM+U*DTT
      IF(T+0.1.GE.TZR)GOTO3
      GOTO2
   3  UG=SUM/(TZR-TWR)
      GOTO6
   5  UG=UGM+1.
   6  RETURN
      END

      SUBROUTINE TWZ(TI,Z,TW,TZ)
      CALL PAR(TI,SC,SV)
      GR=0.01745
      PM=0.2618
      AA=Z*GR
      SCA=SC*GR
      S1=SIN(SCA)*SIN(AA)
      S2=COS(SCA)*COS(AA)
      IF(ABS(S1).GT.S2)GOTO2
      TT=ACO(-S1/S2)/PM
      TW=(SV-TT)*3600.
      TZ=(SV+TT)*3600.
      GOTO10
   2  IF(S1)3,4,4
   3  TW=24.*3600.
      TZ=0.
      GOTO10
   4  TW=0.
      TZ=24.*3600.
  10  RETURN
      END

      SUBROUTINE TWZN(TI,Z,UGF,TW,TZ)
      CALL PAR(TI,SC,SV)
      GR=0.01745
      PM=0.2618
      AA=Z*GR
      SCA=SC*GR
      S1=SIN(SCA)*SIN(AA)-COS(GR*UGF)
      S2=COS(SCA)*COS(AA)
      IF(ABS(S1).GT.S2)GOTO2
      TT=ACO(-S1/S2)/PM
      TW=(SV-TT)*3600.
      TZ=(SV+TT)*3600.
      GOTO10
   2  IF(S1)3,4,4
   3  TW=24.*3600.
      TZ=0.
      GOTO10
   4  TW=0.
      TZ=24.*3600.
  10  RETURN
      END

      SUBROUTINE CSM(TI,Z,UG,T,DT)
      CALL PAR(TI,SC,SV)
      PM=0.2618
      GR=0.01745
      AA=GR*Z
      SCA=SC*GR
      T1=T/3600.
      T2=(T+DT)/3600.
      S1=SIN(AA)*SIN(SCA)
      S2=COS(AA)*COS(SCA)
      UG=ACO(S1+S2*(SIN(PM*(T2-SV))-SIN(PM*(T1-SV)))/(PM*(T2-T1)))/GR
   5  RETURN
      END

      FUNCTION ACO(A)
      IF(A.NE.0.)GOTO1
      AC=0.5*3.14159
      GOTO2
   1  B=SQRT(1.-A*A)
      AC=ATAN(B/A)
      IF(AC.LT.0.)AC=AC+3.14159
   2  ACO=AC
      RETURN
      END

      SUBROUTINE PAR(TI,SC,SV)
      T=2.*3.14159*TI/365.
      T2=2.*T
      T3=3.*T
      SC=0.396372-22.9133*COS(T)-0.387205*COS(T2)-0.154527*
     *COS(T3)+4.02543*SIN(T)+0.0519673*SIN(T2)+0.0847978*SIN(T3)
      SV=12.-0.0072*COS(T)-0.0528*COS(T2)-0.0012*COS(T3)-
     -0.1229*SIN(T)-0.1565*SIN(T2)-0.0041*SIN(T3)
      RETURN
      END

      SUBROUTINE CSM0(TI,Z,UG,T)
      CALL PAR(TI,SC,SV)
      PM=0.2618
      GR=0.01745
      AA=GR*Z
      SCA=SC*GR
      T1=T/3600.
      S1=SIN(AA)*SIN(SCA)
      S2=COS(AA)*COS(SCA)
      UG=ACO(S1+S2*COS(PM*(T1-SV)))/GR
   5  RETURN
      END

      SUBROUTINE CSTU(TI,Z,AM,T,APER)
      CSZ=SQRT((1225./(AM*AM)-1.)/1224.)
      CALL PAR(TI,SC,SV)
      PM=0.2618
      GR=0.01745
      AA=GR*Z
      SCA=SC*GR
      S1=SIN(AA)*SIN(SCA)
      S2=COS(AA)*COS(SCA)
      CZS=(CSZ*APER-S1)/S2
      IF(ABS(CZS).LE.1.)THEN
         T=3600.*(SV+ACO(CZS)/PM)
      ELSE
C         PRINT *,' CSZ=',CSZ,' CZS=',CZS
         T=-1.
      ENDIF
   5  RETURN
      END

      SUBROUTINE CSTM(TI,Z,AM)
      CALL PAR(TI,SC,SV)
      PM=0.2618
      GR=0.01745
      AA=GR*Z
      SCA=SC*GR
      S1=SIN(AA)*SIN(SCA)
      S2=COS(AA)*COS(SCA)
      CZS=S1-S2
      AM=35./(SQRT(1224.*CZS*CZS+1.))
      RETURN
      END

      FUNCTION COSD(Z)
      COSD=COS(3.14159/180.*Z)
      RETURN
      END

      FUNCTION OPTM(Z)
      SZA=COSD(Z)
      OPTM=35./(SQRT(1224.*SZA*SZA+1.))
      RETURN
      END


