updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / var / da / da_test / da_check_balance.inc
blob3d6619a0cb5e2ed3b5400ed0ec473d7d495b5266
1 subroutine da_check_balance(phi, phi_u)
3    !---------------------------------------------------------------------------
4    ! Purpose: Compare balanced mass (phi_b - function of wind) and actual phi.
5    !
6    ! Method:  Calculate correlation between balanced and actual phi.
7    !---------------------------------------------------------------------------
9    implicit none
10       
11    real, intent(in)             :: phi(:,:,:)      ! Total phi.
12    real, intent(in)             :: phi_u(:,:,:)    ! Unbalanced phi.
14    integer                      :: iy              ! Size of 1st dimension.
15    integer                      :: jx              ! Size of 2nd dimension.
16    integer                      :: kz              ! Size of 3rd dimension.
17    integer                      :: i, k            ! Loop counters
18    real                         :: corr_coef       ! Correlation coefficient.
19    real                         :: accurac         ! Accuracy.
20    real, allocatable            :: phi_b1(:)       ! 1D balanced phi.
21    real, allocatable            :: phi_b2(:,:)     ! 2D balanced phi.
22    real, allocatable            :: corr_coeff(:,:) ! Correlation coefficient.
23    real, allocatable            :: accuracy(:,:)   ! Accuracy.
25    if (trace_use) call da_trace_entry("da_check_balance")
26           
27    if (balance_type == balance_geo) then
28       write(unit=stdout, fmt='(a)') ' da_check_balance: Balance is geostrophic.'
29    else if (balance_type == balance_cyc) then
30       write(unit=stdout, fmt='(a)') &
31          ' da_check_balance: Balance is cyclostrophic.'
32    else if (balance_type == balance_geocyc) then
33       write(unit=stdout, fmt='(a)') &
34          ' da_check_balance: Balance is geo/cyclostrophic.'
35    end if
36       
37    write(unit=stdout, fmt='(a)') ' da_check_balance: Correlation/accuracy: '
38       
39    !-------------------------------------------------------------------------
40    ! [1.0]: Initialise:
41    !-------------------------------------------------------------------------  
43    iy = size(phi_u, DIM=1)
44    jx = size(phi_u, DIM=2)
45    kz = size(phi_u, DIM=3)
46       
47    allocate(phi_b1(1:jx))
48    allocate(phi_b2(1:iy,1:jx))
50    allocate(corr_coeff(1:kz,1:iy))
51    corr_coeff(1:kz,1:iy) = 0.0
53    allocate(accuracy(1:kz,1:iy))
54    accuracy(1:kz,1:iy) = 0.0
55       
56    !-------------------------------------------------------------------------
57    ! [2.0]: Calculate correlations/accuracy:
58    !-------------------------------------------------------------------------  
60    do k = 1, kz
61       do i = 1, iy
63          phi_b1(2:jx-1) = phi(i,2:jx-1,k) - phi_u(i,2:jx-1,k)
64             
65          call da_correlation_coeff1d(phi_b1(2:jx-1), phi(i,2:jx-1,k), &
66                                       corr_coeff(k,i), accuracy(k,i))
67      
68          ! write(58,*) corr_coeff(k,i), accuracy(k,i)
69       end do
70          
71       phi_b2(2:iy-1,2:jx-1) = phi(2:iy-1,2:jx-1,k) - phi_u(2:iy-1,2:jx-1,k)
72       call da_correlation_coeff2d(phi_b2(2:iy-1,2:jx-1), &
73                                    phi(2:iy-1,2:jx-1,k), &
74                                    corr_coef, accurac)
76       write(unit=stdout, fmt='(i6,1pe9.2,1pe9.2)') &
77             k, corr_coef, accurac
79    end do
81    !-------------------------------------------------------------------------
82    ! [3.0]: Tidy up:
83    !-------------------------------------------------------------------------  
85    deallocate(phi_b1)
86    deallocate(phi_b2)
87    deallocate(corr_coeff)
88    deallocate(accuracy)
90    if (trace_use) call da_trace_exit("da_check_balance")
92 end subroutine da_check_balance