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)

      _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

      _RS lat

      ampfact(lat,i) = MIN( oneRL, 

     &   ( COS( ABS(lat)*deg2rad )

          lat = yC(i,j,bi,bj)

          dxF(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad

          lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))

          dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad

          lat=0.5 _d 0*(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 _d 0*(yGloc(i,j)+yGloc(i,j+1))

          tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)

          lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))

          tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)

      _RL lat, dlat, dlon, xG0, yG0

       _RL  lat, sun1, dayhrs

          lat=YC(1,j,1,bj)*deg2rad

     &         lat = asin( fCori(1,j,1,bj)/(2. _d 0*omega) )

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

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

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

      _RL lat(3)

      DATA lat / 0. _d 0 ,  35. _d 0 ,  85. _d 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

        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)