SUBROUTINE COSZ0
C****
C**** COSZ0 calculates the Earth's zenith angle, weighted either
C**** by time or by sun light.
C****
      IMPLICIT REAL*8 (A-H,M,O-Z)
      PARAMETER (IM=72,JM=46, TWOPI=6.283185307179586477d0)
      REAL*8 COSZ(IM,JM),COSZA(IM,JM),
     *       SINJ(JM),COSJ(JM),RI(IM),SINI(IM),COSI(IM),LT1,LT2
      COMMON /WORK02/ LT1(IM),SLT1(IM),S2LT1(IM),
     *                LT2(IM),SLT2(IM),S2LT2(IM)
C**** ZERO1 must equal the cut-off value for COSZ used in SOLAR
      PARAMETER (ZERO1=1.d-2)
C**** Compute area weighted latitudes and their sines and cosines
      DLAT  = TWOPI*NINT(360./(JM-1))/720.
      FJEQ  = (1+JM)/2.
      RLATS = -.25*TWOPI
      SLATS = -1.
      CLATS = 0.
      DO 20 J=1,JM
      RLATN = DLAT*(J+.5-FJEQ)
      IF(J.eq.JM)  RLATN = .25*TWOPI
      SLATN = SIN(RLATN)
      CLATN = COS(RLATN)
      RLATM = (RLATN*SLATN+CLATN-RLATS*SLATS-CLATS)/(SLATN-SLATS)
      SINJ(J) = SIN(RLATM)
      COSJ(J) = COS(RLATM)
      RLATS = RLATN
      SLATS = SLATN
   20 CLATS = CLATN
C**** Compute sines and cosines of longitude.  The left edge of
C**** grid box I=1 is assumed to on the International Date Line.
      DO 40 I=1,IM
      RI(I)   = (TWOPI/IM)*(I-.5)
      SINI(I) = SIN(RI(I))
   40 COSI(I) = COS(RI(I))
      RETURN
C****
C****
      ENTRY COSZT (SIND,COSD,ROT1,ROT2,COSZ)
C****
C**** COSZT computes the zenith angle weighted by time from ROT1
C**** to ROT2, Greenwich Mean Time in radians.
C**** Input: SIND,COSD = sine and cosine of the declination angle
C****             ROT1 = initial time,  0 < ROT1 < 2*PI
C****             ROT2 = final time, ROT1 < ROT2 < ROT1+2*PI
C**** Output: COSZ = S(cos(Z)*dT) / S(dT)
C****
      DROT = ROT2-ROT1
C**** Compute sines and cosines of initial and final times
  100 SR1 = SIN(ROT1)
      CR1 = COS(ROT1)
      SR2 = SIN(ROT2)
      CR2 = COS(ROT2)
C**** Compute initial and final local times (measured from noon)
C**** and their sines and cosines
      DO 120 I=1,IM
      LT1(I)  = ROT1+RI(I)
      SLT1(I) = SR1*COSI(I)+CR1*SINI(I)
      LT2(I)  = ROT2+RI(I)
  120 SLT2(I) = SR2*COSI(I)+CR2*SINI(I)
C****
C**** CALCULATION FOR POLAR GRID BOXES
C****
      DO 200 J=1,JM,JM-1
      SJSD = SINJ(J)*SIND
      CJCD = COSJ(J)*COSD
      IF(SJSD+CJCD.LE.ZERO1)  GO TO 180
      IF(SJSD-CJCD.ge.0.)  GO TO 160
C**** AVERAGE COSZ FROM DAWN TO DUSK NEAR THE POLES
      DUSK  = ACOS(-SJSD/CJCD)
      SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD
      DAWN  = -DUSK
      SDAWN = -SDUSK
      COSZ(1,J) = (SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN))/TWOPI
      GO TO 200
C**** CONSTANT DAYLIGHT NEAR THE POLES
  160 COSZ(1,J) = SJSD
      GO TO 200
C**** CONSTANT NIGHTIME NEAR THE POLES
  180 COSZ(1,J) = 0.
  200 continue
C****
C**** LOOP OVER NON-POLAR LATITUDES
C****
      DO 500 J=2,JM-1
      SJSD = SINJ(J)*SIND
      CJCD = COSJ(J)*COSD
      IF(SJSD+CJCD.LE.ZERO1)  GO TO 460
      IF(SJSD-CJCD.ge.0.)  GO TO 420
C**** COMPUTE DAWN AND DUSK (AT LOCAL TIME) AND THEIR SINES
      DUSK  = ACOS(-SJSD/CJCD)
      SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD
      DAWN  = -DUSK
      SDAWN = -SDUSK
C**** NEITHER CONSTANT DAYTIME NOR CONSTANT NIGHTIME AT THIS LATITUDE,
C**** LOOP OVER LONGITUDES
      ZERO2 = ZERO1/CJCD
      DO 400 I=1,IM
C**** FORCE DUSK TO LIE BETWEEN LT1 AND LT1+2*PI
      IF(DUSK.gt.LT1(I)+ZERO2)  GO TO 220
      DUSK = DUSK+TWOPI
      DAWN = DAWN+TWOPI
  220 IF(DAWN.LT.LT2(I)-ZERO2)  GO TO 240
C**** CONTINUOUS NIGHTIME FROM INITIAL TO FINAL TIME
      COSZ(I,J) = 0.
      GO TO 400
  240 IF(DAWN.ge.LT1(I))  GO TO 300
      IF(DUSK.LT.LT2(I))  GO TO 260
C**** CONTINUOUS DAYLIGHT FROM INITIAL TIME TO FINAL TIME
      COSZ(I,J) = SJSD+CJCD*(SLT2(I)-SLT1(I))/DROT
      GO TO 400
  260 IF(DAWN+TWOPI.LT.LT2(I)-ZERO2)  GO TO 280
C**** DAYLIGHT AT INITIAL TIME AND NIGHT AT FINAL TIME
      COSZ(I,J) = (SJSD*(DUSK-LT1(I))+CJCD*(SDUSK-SLT1(I)))/DROT
      GO TO 400
C**** DAYLIGHT AT INITIAL AND FINAL TIMES WITH NIGHTIME IN BETWEEN
  280 COSZ(I,J) = (SJSD*(LT2(I)-DAWN-TWOPI+DUSK-LT1(I))+CJCD*
     *  (SLT2(I)-SDAWN+SDUSK-SLT1(I)))/DROT
      GO TO 400
  300 IF(DUSK.LT.LT2(I))  GO TO 320
C**** NIGHT AT INITIAL TIME AND DAYLIGHT AT FINAL TIME
      COSZ(I,J) = (SJSD*(LT2(I)-DAWN)+CJCD*(SLT2(I)-SDAWN))/DROT
      GO TO 400
C**** NIGHTIME AT INITIAL AND FINAL TIMES WITH DAYLIGHT IN BETWEEN
  320 COSZ(I,J) = (SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN))/DROT
  400 continue
      GO TO 500
C**** CONSTANT DAYLIGHT AT THIS LATITUDE
  420 DO 440 I=1,IM
  440 COSZ(I,J) = SJSD+CJCD*(SLT2(I)-SLT1(I))/DROT
      GO TO 500
C**** CONSTANT NIGHTIME AT THIS LATITUDE
  460 DO 480 I=1,IM
  480 COSZ(I,J) = 0.
  500 continue
      RETURN
C****
C****
      ENTRY COSZS (SIND,COSD,ROT1,ROT2,COSZ,COSZA)
C****
C**** COSZS computes the zenith angle from ROT1 to ROT2 twice,
C**** weighted by time and weighted by sun light.
C**** Input: SIND,COSD = sine and cosine of the declination angle
C****             ROT1 = initial time,  0 < ROT1 < 2*PI
C****             ROT2 = final time, ROT1 < ROT2 < ROT1+2*PI
C**** Output: COSZ  = S(cos(Z)*dT) / S(dT)
C****         COSZA = S(cos(Z)*cos(Z)*dT) / S(cos(Z)*dT)
C****
      DROT = ROT2-ROT1
C**** COMPUTE THE SINES AND COSINES OF THE INITIAL AND FINAL GMT'S
      SR1 = SIN(ROT1)
      CR1 = COS(ROT1)
      SR2 = SIN(ROT2)
      CR2 = COS(ROT2)
C**** COMPUTE THE INITIAL AND FINAL LOCAL TIMES (MEASURED FROM NOON TO
C****   NOON) AND THEIR SINES AND COSINES
      DO 520 I=1,IM
      LT1(I)  = ROT1+RI(I)
      SLT1(I) = SR1*COSI(I)+CR1*SINI(I)
      CLT1    = CR1*COSI(I)-SR1*SINI(I)
      S2LT1(I)= 2.*SLT1(I)*CLT1
      LT2(I)  = ROT2+RI(I)
      SLT2(I) = SR2*COSI(I)+CR2*SINI(I)
      CLT2    = CR2*COSI(I)-SR2*SINI(I)
  520 S2LT2(I)= 2.*SLT2(I)*CLT2
C****
C**** CALCULATION FOR POLAR GRID BOXES
C****
      DO 600 J=1,JM,JM-1
      SJSD = SINJ(J)*SIND
      CJCD = COSJ(J)*COSD
      IF(SJSD+CJCD.LE.ZERO1) GO TO 580
      IF(SJSD-CJCD.ge.0.) GO TO 560
C**** AVERAGE COSZ FROM DAWN TO DUSK NEAR THE POLES
      CDUSK = -SJSD/CJCD
      DUSK  = ACOS(CDUSK)
      SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD
      S2DUSK= 2.*SDUSK*CDUSK
      DAWN  = -DUSK
      SDAWN = -SDUSK
      S2DAWN= -S2DUSK
      ECOSZ = SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN)
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SDAWN)+
     *  .5*CJCD*(DUSK-DAWN+.5*(S2DUSK-S2DAWN)))
      COSZ(1,J)  = ECOSZ/TWOPI
      COSZA(1,J) = QCOSZ/ECOSZ
      GO TO 600
C**** CONSTANT DAYLIGHT NEAR THE POLES
  560 ECOSZ = SJSD*TWOPI
      QCOSZ = SJSD*ECOSZ+.5*CJCD*CJCD*TWOPI
      COSZ(1,J)  = ECOSZ/TWOPI
      COSZA(1,J) = QCOSZ/ECOSZ
      GO TO 600
C**** CONSTANT NIGHTIME NEAR THE POLES
  580 COSZ(1,J)  = 0.
      COSZA(1,J) = 0.
  600 continue
C****
C**** LOOP OVER NON-POLAR LATITUDES
C****
      DO 900 J=2,JM-1
      SJSD = SINJ(J)*SIND
      CJCD = COSJ(J)*COSD
      IF(SJSD+CJCD.LE.ZERO1)  GO TO 860
      IF(SJSD-CJCD.ge.0.)  GO TO 820
C**** COMPUTE DAWN AND DUSK (AT LOCAL TIME) AND THEIR SINES
      CDUSK = -SJSD/CJCD
      DUSK  = ACOS(CDUSK)
      SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD
      S2DUSK= 2.*SDUSK*CDUSK
      DAWN  = -DUSK
      SDAWN = -SDUSK
      S2DAWN= -S2DUSK
C**** NEITHER CONSTANT DAYTIME NOR CONSTANT NIGHTIME AT THIS LATITUDE,
C**** LOOP OVER LONGITUDES
      ZERO2 = ZERO1/CJCD
      DO 800 I=1,IM
C**** FORCE DUSK TO LIE BETWEEN LT1 AND LT1+2*PI
      IF(DUSK.gt.LT1(I)+ZERO2) GO TO 620
      DUSK = DUSK+TWOPI
      DAWN = DAWN+TWOPI
  620 IF(DAWN.LT.LT2(I)-ZERO2) GO TO 640
C**** CONTINUOUS NIGHTIME FROM INITIAL TO FINAL TIME
      COSZ(I,J)  = 0.
      COSZA(I,J) = 0.
      GO TO 800
  640 IF(DAWN.ge.LT1(I))  GO TO 700
      IF(DUSK.LT.LT2(I))  GO TO 660
C**** CONTINUOUS DAYLIGHT FROM INITIAL TIME TO FINAL TIME
      ECOSZ = SJSD*DROT+CJCD*(SLT2(I)-SLT1(I))
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SLT1(I))+
     *  .5*CJCD*(DROT+.5*(S2LT2(I)-S2LT1(I))))
      COSZ(I,J)  = ECOSZ/DROT
      COSZA(I,J) = QCOSZ/ECOSZ
      GO TO 800
  660 IF(DAWN+TWOPI.LT.LT2(I)-ZERO2)  GO TO 680
C**** DAYLIGHT AT INITIAL TIME AND NIGHT AT FINAL TIME
      ECOSZ = SJSD*(DUSK-LT1(I))+CJCD*(SDUSK-SLT1(I))
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SLT1(I))+
     *  .5*CJCD*(DUSK-LT1(I)+.5*(S2DUSK-S2LT1(I))))
      COSZ(I,J)  = ECOSZ/DROT
      COSZA(I,J) = QCOSZ/ECOSZ
      GO TO 800
C**** DAYLIGHT AT INITIAL AND FINAL TIMES WITH NIGHTIME IN BETWEEN
  680 ECOSZ = SJSD*(DROT-DAWN-TWOPI+DUSK)+
     *  CJCD*(SLT2(I)-SDAWN+SDUSK-SLT1(I))
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SLT1(I)+SLT2(I)-SDAWN)+
     *  .5*CJCD*(DUSK+DROT-DAWN-TWOPI+
     *  .5*(S2DUSK-S2LT1(I)+S2LT2(I)-S2DAWN)))
      COSZ(I,J)  = ECOSZ/DROT
      COSZA(I,J) = QCOSZ/ECOSZ
      GO TO 800
  700 IF(DUSK.LT.LT2(I))  GO TO 720
C**** NIGHT AT INITIAL TIME AND DAYLIGHT AT FINAL TIME
      ECOSZ = SJSD*(LT2(I)-DAWN)+CJCD*(SLT2(I)-SDAWN)
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SDAWN)+
     *  .5*CJCD*(LT2(I)-DAWN+.5*(S2LT2(I)-S2DAWN)))
      COSZ(I,J)  = ECOSZ/DROT
      COSZA(I,J) = QCOSZ/ECOSZ
      GO TO 800
C**** NIGHTIME AT INITIAL AND FINAL TIMES WITH DAYLIGHT IN BETWEEN
  720 ECOSZ = SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN)
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SDAWN)+
     *  .5*CJCD*(DUSK-DAWN+.5*(S2DUSK-S2DAWN)))
      COSZ(I,J)  = ECOSZ/DROT
      COSZA(I,J) = QCOSZ/ECOSZ
  800 continue
      GO TO 900
C**** CONSTANT DAYLIGHT AT THIS LATITUDE
  820 DO 840 I=1,IM
      ECOSZ = SJSD*DROT+CJCD*(SLT2(I)-SLT1(I))
      QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SLT1(I))+
     *  .5*CJCD*(DROT+.5*(S2LT2(I)-S2LT1(I))))
      COSZ(I,J)  = ECOSZ/DROT
  840 COSZA(I,J) = QCOSZ/ECOSZ
      GO TO 900
C**** CONSTANT NIGHTIME AT THIS LATITUDE
  860 DO 880 I=1,IM
      COSZ(I,J)  = 0.
  880 COSZA(I,J) = 0.
  900 continue
      RETURN
      END