T_var(1) = SIN(lat(2)*deg2rad) - SIN(lat(1)*deg2rad)

        T_var(2) = SIN(lat(3)*deg2rad) - SIN(lat(2)*deg2rad)

        T_var(3) = SIN(lat(2)*deg2rad)**3. - SIN(lat(1)*deg2rad)**3.

        T_var(4) = SIN(lat(3)*deg2rad)**3. - SIN(lat(2)*deg2rad)**3.

     <        ((3.*SIN(lat(2)*deg2rad)**2. - 1.))

     <        SIN(lat(2)*deg2rad)/rSphere

          Hd35(no_so) = 2.*PI*rSphere*COS(lat(2)*deg2rad)

          Fw35(no_so) = 2.*PI*rSphere*COS(lat(2)*deg2rad)

      _RL lat(3)

      DATA lat / 0.0,  35.0,  85.0 /

          ELSEIF ( yC(i,j,bi,bj) .GE. lat(2) .AND. 

     &       yC(i,j,bi,bj) .LT. lat(3)) THEN

          IF ( yC(i,j,bi,bj) .GE. -lat(3) .AND. 

     &         yC(i,j,bi,bj) .LE. -lat(2)) THEN

          ELSEIF ( yC(i,j,bi,bj) .GT. -lat(2) .AND. 

     &       yC(i,j,bi,bj) .LT. -lat(1)) THEN 

          ELSEIF ( yC(i,j,bi,bj) .GT. lat(1) .AND. 

     &       yC(i,j,bi,bj) .LT. lat(2)) THEN

      _RL alpha, wvNum, lat

       _RL  lat, sun1, dayhrs

          lat=YC(1,j,1,bj)/180.d0*3.1416

          sun1 = -sin(delta)/cos(delta) * sin(lat)/cos(lat)

          cosz = ( sin(delta)*sin(lat)+              !average zenith angle

     &            (cos(delta)*cos(lat)*sin(dayhrs)/dayhrs) )

      do  lat = 1,nlat

      qz(lat,lev,time0) = qz(lat,lev,time0)*voltomas

      integer lat

      do   lat = 1,nlat

      lats(lat) = -90. + (lat-1)*5.

           do lat=1,nlat

           oz(lat,nlev-lev+1,time0) = o3(lat)

      integer lat

      do   lat = 1,nlat

      lats(lat) = -85. + (lat-1)*10.

      do  lat=1,nlat

      read(ku,1000) (qz(lat,lev,time0),lev=1,nlev)

          lat = yC(I,J,bi,bj)

          dXF(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad

          lat = 0.5*(yGloc(I,J)+yGloc(I+1,J))

          dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad

          lat=0.5*(yGloc(I,J)+yGloc(I+1,J))

     &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )

          lat=yC(I,J,bi,bj)-delY( jGl(J,bj) )*0.5 _d 0

     &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))

          lat=yC(I,J,bi,bj)

     &        *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )

          lat=yC(I,J,bi,bj)

     &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))

          IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAs(I,J,bi,bj)=0.

          lat =0.5 _d 0*(yGloc(I,J)+yGloc(I,J+1))

     &     *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )

          IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.

          lat=0.5*(yGloc(I,J)+yGloc(I,J+1))

          tanPhiAtU(I,J,bi,bj)=tan(lat*deg2rad)

          lat=0.5*(yGloc(I,J)+yGloc(I+1,J))

          tanPhiAtV(I,J,bi,bj)=tan(lat*deg2rad)

      _RL lat, dlat, dlon, xG0, yG0

      _RL lat, lat_deg, lon, radius, darea, dradius, dvolume, depth

                  lat = yC(I,J,bi,bj) * pi / 180.0

                xcom = xcom + density * cos(lat) * cos(lon)

                ycom = ycom + density * cos(lat) * sin(lon)

                zcom = zcom + density * sin(lat) *

                xoamc = xoamc + ( v*sin(lon)-u*sin(lat)*cos(lon))

                yoamc = yoamc + (-v*cos(lon)-u*sin(lat)*sin(lon))

                zoamc = zoamc +   u*cos(lat)

                xoamp = xoamp - sin(lat) * cos(lat) * cos(lon)

                yoamp = yoamp - sin(lat) * cos(lat) * sin(lon)

                zoamp = zoamp + cos(lat)**2