! test ! modified by Elena Savenkova 12.04 !----------------------------------- program ave_QBO ! integer,PARAMETER :: ILAT=36,ngit=23,nvar=10 integer,PARAMETER :: ntime=248 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), & & uq_LA(ILAT,ngit),sdv_LA(ILAT,ngit) real :: sdv_EL_sk(ILAT,ngit),sdv_LA_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_LA(ILAT,ngit),sko_av(ILAT,ngit) real :: diff_AB(ILAT,ngit),sdv_AB(ILAT,ngit),t_EXP(ILAT,ngit) real :: tn_u(ILAT,NGIT),tn_s(ILAT,NGIT),p_u(ILAT,NGIT),P_s(ILAT,NGIT) real :: diff_sk(ILAT,ngit),sdv_sk(ILAT,ngit),t_SKO(ILAT,ngit) external betai !-----------------------input files Wly Open (1, file='tp_a0_1980.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (2, file='tp_a0_1983.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (3, file='tp_a0_1993.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (4, file='tp_a0_1995.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (5, file='tp_a0_1997.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (6, file='tp_a0_1999.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (7, file='tp_a0_2002.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (8, file='tp_a0_2004.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (9, file='tp_a0_2008.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (10, file='tp_a0_2013.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) !-----------------------input files Ely Open (11, file='tp_a0_1982.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (12, file='tp_a0_1984.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (13, file='tp_a0_1994.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (14, file='tp_a0_1998.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (15, file='tp_a0_2000.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (16, file='tp_a0_2003.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (17, file='tp_a0_2005.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (18, file='tp_a0_2007.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (19, file='tp_a0_2010.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) Open (20, file='tp_a0_2012.dx', & & access='direct',status='unknown',recl=4*ilat*ngit) !------------------------files for grads Open (33, file='tp0_Jan_Wly_ansmbl.dx', & & access='direct',status='unknown',recl=4*ilat*ngit*3) Open (44, file='tp0_Jan_Ely_ansmbl.dx', & & access='direct',status='unknown',recl=4*ilat*ngit*3) Open (55, file='diffT_Jan_QBO.dx', & & access='direct',status='unknown',recl=4*ilat*ngit*10) Open (66, file='tp0_Jan_QBO_ave.dx', & & access='direct',status='unknown',recl=4*ilat*ngit*2) ! ntime=372. !----------------------read the data do l=1,ntime 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(10,10) !---------------------sko calculation for every ans member 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 sko_uq1(j,k)=sko_uq1(j,k)+ & & (ampl(1,l,j,k)-uq1(j,k))*(ampl(1,l,j,k)-uq1(j,k))/ntime sko_uq2(j,k)=sko_uq2(j,k)+ & & (ampl(2,l,j,k)-uq2(j,k))*(ampl(2,l,j,k)-uq2(j,k))/ntime sko_uq3(j,k)=sko_uq3(j,k)+ & & (ampl(3,l,j,k)-uq3(j,k))*(ampl(3,l,j,k)-uq3(j,k))/ntime sko_uq4(j,k)=sko_uq4(j,k)+ & & (ampl(4,l,j,k)-uq4(j,k))*(ampl(4,l,j,k)-uq4(j,k))/ntime sko_uq5(j,k)=sko_uq5(j,k)+ & & (ampl(5,l,j,k)-uq5(j,k))*(ampl(5,l,j,k)-uq5(j,k))/ntime 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 10 ans members 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)+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 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(33,rec=1) uq_EL, sdv_EL, sko_EL !###################################################################### !---------- END OF AVERAGING westterly 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 ans member 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 sko_uq1(j,k)=sko_uq1(j,k)+ & & (ampl(1,l,j,k)-uq1(j,k))*(ampl(1,l,j,k)-uq1(j,k))/ntime sko_uq2(j,k)=sko_uq2(j,k)+ & & (ampl(2,l,j,k)-uq2(j,k))*(ampl(2,l,j,k)-uq2(j,k))/ntime sko_uq3(j,k)=sko_uq3(j,k)+ & & (ampl(3,l,j,k)-uq3(j,k))*(ampl(3,l,j,k)-uq3(j,k))/ntime sko_uq4(j,k)=sko_uq4(j,k)+ & & (ampl(4,l,j,k)-uq4(j,k))*(ampl(4,l,j,k)-uq4(j,k))/ntime sko_uq5(j,k)=sko_uq5(j,k)+ & & (ampl(5,l,j,k)-uq5(j,k))*(ampl(5,l,j,k)-uq5(j,k))/ntime 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 10 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)+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 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 stop end program FUNCTION betai(a,b,x) REAL betai,a,b,x !CU USES betacf,gammln REAL bt,betacf,gammln if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai' if(x.eq.0..or.x.eq.1.)then bt=0. else bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x)) endif if(x.lt.(a+1.)/(a+b+2.))then betai=bt*betacf(a,b,x)/a return else betai=1.-bt*betacf(b,a,1.-x)/b return endif END FUNCTION betacf(a,b,x) INTEGER MAXIT REAL betacf,a,b,x,EPS,FPMIN PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) INTEGER m,m2 REAL aa,c,d,del,h,qab,qam,qap qab=a+b qap=a+1. qam=a-1. c=1. d=1.-qab*x/qap if(abs(d).lt.FPMIN)d=FPMIN d=1./d h=d do 11 m=1,MAXIT m2=2*m aa=m*(b-m)*x/((qam+m2)*(a+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d h=h*d*c aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue pause 'a or b too big, or MAXIT too small in betacf' 1 betacf=h return END FUNCTION gammln(xx) REAL gammln,xx INTEGER j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & & -.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END