! modified by Alex Pogoreltsev 11.05.16
!-----------------------------------
 program ave_QBO_JRA55
!

   integer,PARAMETER :: ILAT=36,ngit=17,nvar=5
   integer,PARAMETER :: ntime=744
   real :: ampl(nvar,ntime,ilat,ngit)
   real :: q1(ILAT,ngit),q2(ILAT,ngit),       &
  &        q3(ILAT,ngit),q4(ILAT,ngit),       &
  &        q5(ILAT,ngit),q6(ILAT,ngit),       &
  &        q7(ILAT,ngit),q8(ILAT,ngit),       &
  &        q9(ILAT,ngit),q10(ILAT,ngit)

   real :: uq1(ILAT,ngit),sko_uq1(ILAT,ngit), &
  &        uq2(ILAT,ngit),sko_uq2(ILAT,ngit), &
  &        uq3(ILAT,ngit),sko_uq3(ILAT,ngit), &
  &        uq4(ILAT,ngit),sko_uq4(ILAT,ngit), &
  &        uq5(ILAT,ngit),sko_uq5(ILAT,ngit), &
  &        uq6(ILAT,ngit),sko_uq6(ILAT,ngit), &
  &        uq7(ILAT,ngit),sko_uq7(ILAT,ngit), &
  &        uq8(ILAT,ngit),sko_uq8(ILAT,ngit), &
  &        uq9(ILAT,ngit),sko_uq9(ILAT,ngit), &
  &        uq10(ILAT,ngit),sko_uq10(ILAT,ngit)
   real :: uq_av(ILAT,ngit),uq_EL(ILAT,ngit),sdv_EL(ILAT,ngit)
                          
   real :: sdv_EL_sk(ILAT,ngit)
   real :: dsp_uq1(ILAT,ngit),dsp_uq2(ILAT,ngit),dsp_uq3(ILAT,ngit), &
  &        dsp_uq4(ILAT,ngit),dsp_uq5(ILAT,ngit),dsp_uq6(ILAT,ngit), &
  &        dsp_uq7(ILAT,ngit),dsp_uq8(ILAT,ngit),dsp_uq9(ILAT,ngit), &
  &        dsp_uq10(ILAT,ngit)
   real :: dsp_sk1(ILAT,ngit),dsp_sk2(ILAT,ngit),dsp_sk3(ILAT,ngit), &
  &        dsp_sk4(ILAT,ngit),dsp_sk5(ILAT,ngit),dsp_sk6(ILAT,ngit), &
  &        dsp_sk7(ILAT,ngit),dsp_sk8(ILAT,ngit),dsp_sk9(ILAT,ngit), &
  &        dsp_sk10(ILAT,ngit)
   real :: sko_EL(ILAT,ngit),sko_av(ILAT,ngit)



!-----------------------input file 1958-2019 monmean
          Open (11, file='tp_a0_grads.dx', &
  &   access='direct',status='unknown',recl=4*ilat*ngit)

!------------------------files for grads
          Open (22, file='tp_a0_jan.dx', &
  &   access='direct',status='unknown',recl=4*ilat*ngit)
          Open (33, file='tp_a0_sdv_sko_jan.dx', &
  &   access='direct',status='unknown',recl=4*ilat*ngit*2)

!----------------------read the data
      
      do l=1,ntime

          read(11,rec=625) q1
          read(11,rec=637) q2
          read(11,rec=649) q3
          read(11,rec=661) q4
          read(11,rec=673) q5
          read(11,rec=685) q6
          read(11,rec=697) q7
          read(11,rec=709) q8
          read(11,rec=721) q9
          read(11,rec=733) q10

!          read(1,rec=l+488) q1
!          read(2,rec=l+488) q2
!          read(3,rec=l+488) q3
!          read(4,rec=l+488) q4
!          read(5,rec=l+488) q5
!          read(6,rec=l+488) q6
!          read(7,rec=l+488) q7
!          read(8,rec=l+488) q8
!          read(9,rec=l+488) q9
!          read(10,rec=l+488) q10

      do k=1,ngit
       do j=1,ilat 
       
       ampl(1,l,j,k)=q1(j,k)
       ampl(2,l,j,k)=q2(j,k)
       ampl(3,l,j,k)=q3(j,k)
       ampl(4,l,j,k)=q4(j,k)
       ampl(5,l,j,k)=q5(j,k)
       ampl(6,l,j,k)=q6(j,k)
       ampl(7,l,j,k)=q7(j,k)
       ampl(8,l,j,k)=q8(j,k)
       ampl(9,l,j,k)=q9(j,k)
       ampl(10,l,j,k)=q10(j,k)

       end do
      end do
      end do !time

!-------------------averaging over time period
     do k=1,ngit
       do j=1,ilat

       uq1(j,k)=0.
       uq2(j,k)=0.
       uq3(j,k)=0.
       uq4(j,k)=0.
       uq5(j,k)=0.
       uq6(j,k)=0.
       uq7(j,k)=0.
       uq8(j,k)=0.
       uq9(j,k)=0.
       uq10(j,k)=0.
       
       do l=1,ntime
        
        uq1(j,k)=uq1(j,k)+ampl(1,l,j,k)/ntime
        uq2(j,k)=uq2(j,k)+ampl(2,l,j,k)/ntime
        uq3(j,k)=uq3(j,k)+ampl(3,l,j,k)/ntime
        uq4(j,k)=uq4(j,k)+ampl(4,l,j,k)/ntime
        uq5(j,k)=uq5(j,k)+ampl(5,l,j,k)/ntime
        uq6(j,k)=uq6(j,k)+ampl(6,l,j,k)/ntime
        uq7(j,k)=uq7(j,k)+ampl(7,l,j,k)/ntime
        uq8(j,k)=uq8(j,k)+ampl(8,l,j,k)/ntime
        uq9(j,k)=uq9(j,k)+ampl(9,l,j,k)/ntime
        uq10(j,k)=uq10(j,k)+ampl(10,l,j,k)/ntime 
       end do
      end do
     end do
!      print*, uq1(5,5)
!---------------------sko calculation for every year
     do k=1,ngit
       do j=1,ilat

       sko_uq1(j,k)=0.
       sko_uq2(j,k)=0.
       sko_uq3(j,k)=0.
       sko_uq4(j,k)=0.
       sko_uq5(j,k)=0.
       sko_uq6(j,k)=0.
       sko_uq7(j,k)=0.
       sko_uq8(j,k)=0.
       sko_uq9(j,k)=0.
       sko_uq10(j,k)=0.

       do l=1,ntime
       fl_ntm1=float(ntime-1) 
        sko_uq1(j,k)=sko_uq1(j,k)+ &
  &     (ampl(1,l,j,k)-uq1(j,k))*(ampl(1,l,j,k)-uq1(j,k))/fl_ntm1

        sko_uq2(j,k)=sko_uq2(j,k)+ &
  &     (ampl(2,l,j,k)-uq2(j,k))*(ampl(2,l,j,k)-uq2(j,k))/fl_ntm1

        sko_uq3(j,k)=sko_uq3(j,k)+ &
  &     (ampl(3,l,j,k)-uq3(j,k))*(ampl(3,l,j,k)-uq3(j,k))/fl_ntm1

        sko_uq4(j,k)=sko_uq4(j,k)+ &
  &     (ampl(4,l,j,k)-uq4(j,k))*(ampl(4,l,j,k)-uq4(j,k))/fl_ntm1

        sko_uq5(j,k)=sko_uq5(j,k)+ &
  &     (ampl(5,l,j,k)-uq5(j,k))*(ampl(5,l,j,k)-uq5(j,k))/fl_ntm1

        sko_uq6(j,k)=sko_uq6(j,k)+ &
  &     (ampl(6,l,j,k)-uq6(j,k))*(ampl(6,l,j,k)-uq6(j,k))/ntime

        sko_uq7(j,k)=sko_uq7(j,k)+ &
  &     (ampl(7,l,j,k)-uq7(j,k))*(ampl(7,l,j,k)-uq7(j,k))/ntime

        sko_uq8(j,k)=sko_uq8(j,k)+ &
  &     (ampl(8,l,j,k)-uq8(j,k))*(ampl(8,l,j,k)-uq8(j,k))/ntime

        sko_uq9(j,k)=sko_uq9(j,k)+ &
  &     (ampl(9,l,j,k)-uq9(j,k))*(ampl(9,l,j,k)-uq9(j,k))/ntime

        sko_uq10(j,k)=sko_uq10(j,k)+ &
  &     (ampl(10,l,j,k)-uq10(j,k))*(ampl(10,l,j,k)-uq10(j,k))/ntime

      end do ! time

      sko_uq1(j,k)=sqrt(sko_uq1(j,k))
      sko_uq2(j,k)=sqrt(sko_uq2(j,k))
      sko_uq3(j,k)=sqrt(sko_uq3(j,k))
      sko_uq4(j,k)=sqrt(sko_uq4(j,k))
      sko_uq5(j,k)=sqrt(sko_uq5(j,k))
      sko_uq6(j,k)=sqrt(sko_uq6(j,k))
      sko_uq7(j,k)=sqrt(sko_uq7(j,k))
      sko_uq8(j,k)=sqrt(sko_uq8(j,k))
      sko_uq9(j,k)=sqrt(sko_uq9(j,k))
      sko_uq10(j,k)=sqrt(sko_uq10(j,k))
     end do
     end do


!---------------------averaging over years

     do k=1,ngit
       do j=1,ilat   

       sko_av(j,k)=(sko_uq1(j,k)+sko_uq2(j,k)+sko_uq3(j,k)+ &
  &                 sko_uq4(j,k)+sko_uq5(j,k))/5.
     
       uq_av(j,k) =(uq1(j,k)+uq2(j,k)+uq3(j,k)+uq4(j,k)+ &
  &                 uq5(j,k))/5.

       sko_av(j,k)=(sko_uq1(j,k)+sko_uq2(j,k)+sko_uq3(j,k)+ &
  &                 sko_uq4(j,k)+sko_uq5(j,k)+sko_uq6(j,k)+ &
  &                 sko_uq7(j,k)+sko_uq8(j,k)+sko_uq9(j,k)+ &
  &                 sko_uq10(j,k))/10.
     
       uq_av(j,k) =(uq1(j,k)+uq2(j,k)+uq3(j,k)+uq4(j,k)+ &
  &           uq5(j,k)+uq6(j,k)+uq7(j,k)+uq8(j,k)+uq9(j,k)+ &
  &           uq10(j,k))/10.

       end do
     end do
!-------------------SDV interyear calculation
     do k=1,ngit
       do j=1,ilat
        dsp_uq1(j,k)=(uq1(j,k)-uq_av(j,k))*(uq1(j,k)-uq_av(j,k))
        dsp_uq2(j,k)=(uq2(j,k)-uq_av(j,k))*(uq2(j,k)-uq_av(j,k))
        dsp_uq3(j,k)=(uq3(j,k)-uq_av(j,k))*(uq3(j,k)-uq_av(j,k))
        dsp_uq4(j,k)=(uq4(j,k)-uq_av(j,k))*(uq4(j,k)-uq_av(j,k))
        dsp_uq5(j,k)=(uq5(j,k)-uq_av(j,k))*(uq5(j,k)-uq_av(j,k))
        dsp_uq6(j,k)=(uq6(j,k)-uq_av(j,k))*(uq6(j,k)-uq_av(j,k))
        dsp_uq7(j,k)=(uq7(j,k)-uq_av(j,k))*(uq7(j,k)-uq_av(j,k))
        dsp_uq8(j,k)=(uq8(j,k)-uq_av(j,k))*(uq8(j,k)-uq_av(j,k))
        dsp_uq9(j,k)=(uq9(j,k)-uq_av(j,k))*(uq9(j,k)-uq_av(j,k))
        dsp_uq10(j,k)=(uq10(j,k)-uq_av(j,k))*(uq10(j,k)-uq_av(j,k))
!-----------------------------------------------------------------
        dsp_sk1(j,k)=(sko_uq1(j,k)-sko_av(j,k))*(sko_uq1(j,k)-sko_av(j,k))
        dsp_sk2(j,k)=(sko_uq2(j,k)-sko_av(j,k))*(sko_uq2(j,k)-sko_av(j,k))
        dsp_sk3(j,k)=(sko_uq3(j,k)-sko_av(j,k))*(sko_uq3(j,k)-sko_av(j,k))
        dsp_sk4(j,k)=(sko_uq4(j,k)-sko_av(j,k))*(sko_uq4(j,k)-sko_av(j,k))
        dsp_sk5(j,k)=(sko_uq5(j,k)-sko_av(j,k))*(sko_uq5(j,k)-sko_av(j,k))
        dsp_sk6(j,k)=(sko_uq6(j,k)-sko_av(j,k))*(sko_uq6(j,k)-sko_av(j,k))
        dsp_sk7(j,k)=(sko_uq7(j,k)-sko_av(j,k))*(sko_uq7(j,k)-sko_av(j,k))
        dsp_sk8(j,k)=(sko_uq8(j,k)-sko_av(j,k))*(sko_uq8(j,k)-sko_av(j,k))
        dsp_sk9(j,k)=(sko_uq9(j,k)-sko_av(j,k))*(sko_uq9(j,k)-sko_av(j,k))
        dsp_sk10(j,k)=(sko_uq10(j,k)-sko_av(j,k))*(sko_uq10(j,k)-sko_av(j,k))
       end do
     end do
     do k=1,ngit
       do j=1,ilat
       uq_EL(j,k)=uq_av(j,k)
       sko_EL(j,k)=sko_av(j,k)

        sdv_EL(j,k)=sqrt((dsp_uq1(j,k)+dsp_uq2(j,k)+dsp_uq3(j,k)+ &
  &              dsp_uq4(j,k)+dsp_uq5(j,k)+dsp_uq6(j,k)+dsp_uq7(j,k)+ &
  &              dsp_uq8(j,k)+dsp_uq9(j,k)+dsp_uq10(j,k))/9.)
        sdv_EL_sk(j,k)=sqrt((dsp_sk1(j,k)+dsp_sk2(j,k)+dsp_sk3(j,k)+ &
  &               dsp_sk4(j,k)+dsp_sk5(j,k)+dsp_sk6(j,k)+dsp_sk7(j,k)+ &
  &               dsp_sk8(j,k)+dsp_sk9(j,k)+dsp_sk10(j,k))/9.)
       end do
    end do
          write(22,rec=1) uq_EL
          write(33,rec=1) sdv_EL, sko_EL

       STOP
       END     























!######################################################################

!---------- END OF AVERAGING Westerly QBO


!----------------------read the data

      do l=1,ntime

          read(11,rec=l+488) q1
          read(12,rec=l+488) q2
          read(13,rec=l+488) q3
          read(14,rec=l+488) q4
          read(15,rec=l+488) q5
          read(16,rec=l+488) q6
          read(17,rec=l+488) q7
          read(18,rec=l+488) q8
          read(19,rec=l+488) q9
          read(20,rec=l+488) q10

      do k=1,ngit
       do j=1,ilat

       ampl(1,l,j,k)=q1(j,k)
       ampl(2,l,j,k)=q2(j,k)
       ampl(3,l,j,k)=q3(j,k)
       ampl(4,l,j,k)=q4(j,k)
       ampl(5,l,j,k)=q5(j,k)
       ampl(6,l,j,k)=q6(j,k)
       ampl(7,l,j,k)=q7(j,k)
       ampl(8,l,j,k)=q8(j,k)
       ampl(9,l,j,k)=q9(j,k)
       ampl(10,l,j,k)=q10(j,k)

       end do
      end do
      end do !time

!---------------------sko calculation for every year
     do k=1,ngit
       do j=1,ilat

       uq1(j,k)=0.
       uq2(j,k)=0.
       uq3(j,k)=0.
       uq4(j,k)=0.
       uq5(j,k)=0.
       uq6(j,k)=0.
       uq7(j,k)=0.
       uq8(j,k)=0.
       uq9(j,k)=0.
       uq10(j,k)=0.
       
       do l=1,ntime
        
        uq1(j,k)=uq1(j,k)+ampl(1,l,j,k)/ntime
        uq2(j,k)=uq2(j,k)+ampl(2,l,j,k)/ntime
        uq3(j,k)=uq3(j,k)+ampl(3,l,j,k)/ntime
        uq4(j,k)=uq4(j,k)+ampl(4,l,j,k)/ntime
        uq5(j,k)=uq5(j,k)+ampl(5,l,j,k)/ntime
        uq6(j,k)=uq6(j,k)+ampl(6,l,j,k)/ntime
        uq7(j,k)=uq7(j,k)+ampl(7,l,j,k)/ntime
        uq8(j,k)=uq8(j,k)+ampl(8,l,j,k)/ntime
        uq9(j,k)=uq9(j,k)+ampl(9,l,j,k)/ntime
        uq10(j,k)=uq10(j,k)+ampl(10,l,j,k)/ntime 
       end do
      end do
     end do

     do k=1,ngit
       do j=1,ilat

       sko_uq1(j,k)=0.
       sko_uq2(j,k)=0.
       sko_uq3(j,k)=0.
       sko_uq4(j,k)=0.
       sko_uq5(j,k)=0.
       sko_uq6(j,k)=0.
       sko_uq7(j,k)=0.
       sko_uq8(j,k)=0.
       sko_uq9(j,k)=0.
       sko_uq10(j,k)=0.

       do l=1,ntime
       fl_ntm1=float(ntime-1) 
        sko_uq1(j,k)=sko_uq1(j,k)+ &
  &     (ampl(1,l,j,k)-uq1(j,k))*(ampl(1,l,j,k)-uq1(j,k))/fl_ntm1

        sko_uq2(j,k)=sko_uq2(j,k)+ &
  &     (ampl(2,l,j,k)-uq2(j,k))*(ampl(2,l,j,k)-uq2(j,k))/fl_ntm1

        sko_uq3(j,k)=sko_uq3(j,k)+ &
  &     (ampl(3,l,j,k)-uq3(j,k))*(ampl(3,l,j,k)-uq3(j,k))/fl_ntm1

        sko_uq4(j,k)=sko_uq4(j,k)+ &
  &     (ampl(4,l,j,k)-uq4(j,k))*(ampl(4,l,j,k)-uq4(j,k))/fl_ntm1

        sko_uq5(j,k)=sko_uq5(j,k)+ &
  &     (ampl(5,l,j,k)-uq5(j,k))*(ampl(5,l,j,k)-uq5(j,k))/fl_ntm1

        sko_uq6(j,k)=sko_uq6(j,k)+ &
  &     (ampl(6,l,j,k)-uq6(j,k))*(ampl(6,l,j,k)-uq6(j,k))/ntime

        sko_uq7(j,k)=sko_uq7(j,k)+ &
  &     (ampl(7,l,j,k)-uq7(j,k))*(ampl(7,l,j,k)-uq7(j,k))/ntime

        sko_uq8(j,k)=sko_uq8(j,k)+ &
  &     (ampl(8,l,j,k)-uq8(j,k))*(ampl(8,l,j,k)-uq8(j,k))/ntime

        sko_uq9(j,k)=sko_uq9(j,k)+ &
  &     (ampl(9,l,j,k)-uq9(j,k))*(ampl(9,l,j,k)-uq9(j,k))/ntime

        sko_uq10(j,k)=sko_uq10(j,k)+ &
  &     (ampl(10,l,j,k)-uq10(j,k))*(ampl(10,l,j,k)-uq10(j,k))/ntime

      end do ! time

      sko_uq1(j,k)=sqrt(sko_uq1(j,k))
      sko_uq2(j,k)=sqrt(sko_uq2(j,k))
      sko_uq3(j,k)=sqrt(sko_uq3(j,k))
      sko_uq4(j,k)=sqrt(sko_uq4(j,k))
      sko_uq5(j,k)=sqrt(sko_uq5(j,k))
      sko_uq6(j,k)=sqrt(sko_uq6(j,k))
      sko_uq7(j,k)=sqrt(sko_uq7(j,k))
      sko_uq8(j,k)=sqrt(sko_uq8(j,k))
      sko_uq9(j,k)=sqrt(sko_uq9(j,k))
      sko_uq10(j,k)=sqrt(sko_uq10(j,k))
     end do
     end do

!---------------------averaging over years

     do k=1,ngit
       do j=1,ilat   

       sko_av(j,k)=(sko_uq1(j,k)+sko_uq2(j,k)+sko_uq3(j,k)+ &
  &                 sko_uq4(j,k)+sko_uq5(j,k))/5.
     
       uq_av(j,k) =(uq1(j,k)+uq2(j,k)+uq3(j,k)+uq4(j,k)+ &
  &                 uq5(j,k))/5.

       sko_av(j,k)=(sko_uq1(j,k)+sko_uq2(j,k)+sko_uq3(j,k)+ &
  &                 sko_uq4(j,k)+sko_uq5(j,k)+sko_uq6(j,k)+ &
  &                 sko_uq7(j,k)+sko_uq8(j,k)+sko_uq9(j,k)+ &
  &                 sko_uq10(j,k))/10.
     
       uq_av(j,k) =(uq1(j,k)+uq2(j,k)+uq3(j,k)+uq4(j,k)+ &
  &           uq5(j,k)+uq6(j,k)+uq7(j,k)+uq8(j,k)+uq9(j,k)+ &
  &           uq10(j,k))/10.

       end do
     end do
!-------------------SDV interyear calculation
     do k=1,ngit
       do j=1,ilat
        dsp_uq1(j,k)=(uq1(j,k)-uq_av(j,k))*(uq1(j,k)-uq_av(j,k))
        dsp_uq2(j,k)=(uq2(j,k)-uq_av(j,k))*(uq2(j,k)-uq_av(j,k))
        dsp_uq3(j,k)=(uq3(j,k)-uq_av(j,k))*(uq3(j,k)-uq_av(j,k))
        dsp_uq4(j,k)=(uq4(j,k)-uq_av(j,k))*(uq4(j,k)-uq_av(j,k))
        dsp_uq5(j,k)=(uq5(j,k)-uq_av(j,k))*(uq5(j,k)-uq_av(j,k))
        dsp_uq6(j,k)=(uq6(j,k)-uq_av(j,k))*(uq6(j,k)-uq_av(j,k))
        dsp_uq7(j,k)=(uq7(j,k)-uq_av(j,k))*(uq7(j,k)-uq_av(j,k))
        dsp_uq8(j,k)=(uq8(j,k)-uq_av(j,k))*(uq8(j,k)-uq_av(j,k))
        dsp_uq9(j,k)=(uq9(j,k)-uq_av(j,k))*(uq9(j,k)-uq_av(j,k))
        dsp_uq10(j,k)=(uq10(j,k)-uq_av(j,k))*(uq10(j,k)-uq_av(j,k))
!-----------------------------------------------------------------
        dsp_sk1(j,k)=(sko_uq1(j,k)-sko_av(j,k))*(sko_uq1(j,k)-sko_av(j,k))
        dsp_sk2(j,k)=(sko_uq2(j,k)-sko_av(j,k))*(sko_uq2(j,k)-sko_av(j,k))
        dsp_sk3(j,k)=(sko_uq3(j,k)-sko_av(j,k))*(sko_uq3(j,k)-sko_av(j,k))
        dsp_sk4(j,k)=(sko_uq4(j,k)-sko_av(j,k))*(sko_uq4(j,k)-sko_av(j,k))
        dsp_sk5(j,k)=(sko_uq5(j,k)-sko_av(j,k))*(sko_uq5(j,k)-sko_av(j,k))
        dsp_sk6(j,k)=(sko_uq6(j,k)-sko_av(j,k))*(sko_uq6(j,k)-sko_av(j,k))
        dsp_sk7(j,k)=(sko_uq7(j,k)-sko_av(j,k))*(sko_uq7(j,k)-sko_av(j,k))
        dsp_sk8(j,k)=(sko_uq8(j,k)-sko_av(j,k))*(sko_uq8(j,k)-sko_av(j,k))
        dsp_sk9(j,k)=(sko_uq9(j,k)-sko_av(j,k))*(sko_uq9(j,k)-sko_av(j,k))
        dsp_sk10(j,k)=(sko_uq10(j,k)-sko_av(j,k))*(sko_uq10(j,k)-sko_av(j,k))
       end do
     end do
     do k=1,ngit
       do j=1,ilat
       uq_LA(j,k)=uq_av(j,k)
       sko_LA(j,k)=sko_av(j,k)

        sdv_LA(j,k)=sqrt((dsp_uq1(j,k)+dsp_uq2(j,k)+dsp_uq3(j,k)+ &
  &              dsp_uq4(j,k)+dsp_uq5(j,k)+dsp_uq6(j,k)+dsp_uq7(j,k)+ &
  &              dsp_uq8(j,k)+dsp_uq9(j,k)+dsp_uq10(j,k))/9.)
        sdv_LA_sk(j,k)=sqrt((dsp_sk1(j,k)+dsp_sk2(j,k)+dsp_sk3(j,k)+ &
  &               dsp_sk4(j,k)+dsp_sk5(j,k)+dsp_sk6(j,k)+dsp_sk7(j,k)+ &
  &               dsp_sk8(j,k)+dsp_sk9(j,k)+dsp_sk10(j,k))/9.)

!---------- END OF AVERAGING Easterly QBO
!-------------------------------  t_TEST
     sdv_AB(j,k)= &
  &   SQRT((10.*sdv_EL(j,k)*sdv_EL(j,k)+10.*sdv_LA(j,k)*sdv_LA(j,k))/20.)
     diff_AB(j,k)=uq_EL(j,k)-uq_LA(j,k) 
     t_EXP(j,k) = ABS(diff_AB(j,k))/sdv_AB(j,k)/SQRT(1./10.+1./10.)
     tn_u(j,k)=9.* &
  & (sdv_EL(j,k)*sdv_EL(j,k)+sdv_LA(j,k)*sdv_LA(j,k))* &
  & (sdv_EL(j,k)*sdv_EL(j,k)+sdv_LA(j,k)*sdv_LA(j,k))/ &
  & (sdv_EL(j,k)*sdv_EL(j,k)*sdv_EL(j,k)*sdv_EL(j,k)+  &
  &  sdv_LA(j,k)*sdv_LA(j,k)*sdv_LA(j,k)*sdv_LA(j,k))
!-----------------------------------------------------------------------
     sdv_sk(j,k)= &
  &   SQRT((10.*sdv_EL_sk(j,k)*sdv_EL_sk(j,k)+10.*sdv_LA_sk(j,k)*sdv_LA_sk(j,k))/20.)
     diff_sk(j,k)=sko_EL(j,k)-sko_LA(j,k)
     t_SKO(j,k) = ABS(diff_sk(j,k))/sdv_sk(j,k)/SQRT(1./10.+1./10.)
     tn_s(j,k)=9.* &
  & (sdv_EL_sk(j,k)*sdv_EL_sk(j,k)+sdv_LA_sk(j,k)*sdv_LA_sk(j,k))* &
  & (sdv_EL_sk(j,k)*sdv_EL_sk(j,k)+sdv_LA_sk(j,k)*sdv_LA_sk(j,k))/ & 
  & (sdv_EL_sk(j,k)*sdv_EL_sk(j,k)*sdv_EL_sk(j,k)*sdv_EL_sk(j,k)+  &
  &  sdv_LA_sk(j,k)*sdv_LA_sk(j,k)*sdv_LA_sk(j,k)*sdv_LA_sk(j,k))
!---------------------------------------------- Significance
    x=tn_u(j,k)/(tn_u(j,k)+t_EXP(j,k)*t_EXP(j,k))
    a=tn_u(j,k)/2.
    b=0.5
    p_u(j,k)=betai(a,b,x)
    x=tn_s(j,k)/(tn_s(j,k)+t_SKO(j,k)*t_SKO(j,k))
    a=tn_s(j,k)/2.
    b=0.5
    p_s(j,k)=betai(a,b,x)      
       end do
    end do
          write(44,rec=1) uq_LA, sdv_LA, sko_LA
          write(55,rec=1) diff_AB,sdv_AB,t_EXP,tn_u,p_u,diff_sk,sdv_sk,t_SKO,tn_s,p_s
          write(66,rec=1) uq_av, diff_AB