Merge remote-tracking branch 'origin/release-v4.5.2'
[WRF.git] / var / da / da_setup_structures / da_interpolate_regcoeff.inc
blob5da1057630c8b291098db6048a7969db8ef388a2
1 subroutine da_interpolate_regcoeff (iy, iys, kz, kzs, meanl_stats, meanl_xb, meanp_stats, meanp_xb, &
2    pb_vert_reg_stats, pb_vert_reg)
4    !---------------------------------------------------------------------------
5    ! Purpose: Interpolate statistical regression coefficient to new domain.
6    !
7    ! Method:  i,k,k Interpolation.
8    !---------------------------------------------------------------------------
10    implicit none
12    integer, intent(in)  :: iy                       ! Number of rows in  xb.
13    integer, intent(in)  :: iys                      ! Number of rows in stats.
14    integer, intent(in)  :: kz                       ! Number of levels in xb.
15    integer, intent(in)  :: kzs                      ! Number of levels in stats.
16    real,    intent(in)  :: meanl_stats(:)           ! Mean latitude on stats rows.
17    real,    intent(in)  :: meanl_xb(:)              ! Mean latitude on xb rows.
18    real,    intent(in)  :: meanp_stats(:)           ! Mean pressure on stats levs.
19    real,    intent(in)  :: meanp_xb(:)              ! Mean pressure on xb levs.
20    real*8,  intent(in)  :: pb_vert_reg_stats(:,:,:) ! Coefficient on stats grid.
21    real*8,  intent(out) :: pb_vert_reg(:,:,:)       ! Coefficient on xb grid.
22      
23    integer :: i, is, i_south           ! Loop counters.
24    integer :: k1, k2, k, ks            ! Loop counters.
25    integer :: k1s, k2s
26    real    :: lat_wgt
28    integer :: k_above(1:kz)
29    real    :: pb_vert_reg_temp(1:iys,1:kz,1:kz)
30    real    :: weight(1:kz)
32    if (trace_use) call da_trace_entry("da_interpolate_regcoeff")
34    pb_vert_reg = 0.0
36    !------------------------------------------------------------------------
37    ! [1.0] Find xb levels/rows bounded by stats domain:
38    !------------------------------------------------------------------------
40    do k = 1, kz
41       if (meanp_xb(k) <= meanp_stats(1)) then
42          weight(k) = 1.0e-6
43          k_above(k) = 1
44       else if (meanp_xb(k) >= meanp_stats(kzs)) then
45          weight(k) = 1.0-1.0e-6
46          k_above(k) = kzs-1
47       else
48          do ks = 1, kzs-1
49             if (meanp_xb(k) >= meanp_stats(ks) .AND. meanp_xb(k) <= meanp_stats(ks+1)) then
50                weight(k) = (meanp_xb(k) - meanp_stats(ks)) / (meanp_stats(ks+1) - meanp_stats(ks))
51                k_above(k) = ks
52                exit
53             end if
54          end do
55       end if
56    end do
58    !------------------------------------------------------------------------   
59    ! [3.0] Interpolate regression coefficient from stats to xb levels:
60    !------------------------------------------------------------------------
62    pb_vert_reg_temp(1:iys,1:kz,1:kz) = 0.0
64    do is = 1, iys
65       do k1 = 1, kz
66          k1s = k_above(k1)
67          do k2 = 1, kz
68             k2s = k_above(k2)
70             pb_vert_reg_temp(is,k1,k2) = (1.0-weight(k1)) * (1.0-weight(k2)) * &
71                                          pb_vert_reg_stats(is,k1s,k2s) + &
72                                          weight(k1) * (1.0-weight(k2)) * &
73                                          pb_vert_reg_stats(is,k1s+1,k2s) + &
74                                          weight(k2) * (1.0-weight(k1)) * &
75                                          pb_vert_reg_stats(is,k1s,k2s+1) + &
76                                          weight(k2) * weight(k1) * &
77                                          pb_vert_reg_stats(is,k1s+1,k2s+1)
78          end do
79       end do     
80    end do
81          
82    !------------------------------------------------------------------------
83    ! [4.0] Interpolate to from statistics latitudes to xb latitudes:
84    !------------------------------------------------------------------------
86    i_south = 2
88    do i = 1, iy
89    
90       ! Find position of xb latitude in statistics rows:
92       if (meanl_xb(i) <= meanl_stats(2)) then
93          i_south = 2
94          lat_wgt = 0.5
95       else if (meanl_xb(i) >= meanl_stats(iys-1)) then
96          i_south = iys-2
97          lat_wgt = 0.5
98       else
99          do is = 1, iys-1
100             if (meanl_xb(i) >= meanl_stats(is) .AND. meanl_xb(i) <= meanl_stats(is+1)) then
102                lat_wgt = (meanl_xb(i) - meanl_stats(is)) / (meanl_stats(is+1) - meanl_stats(is))
103                i_south = is
104                exit
105             end if
106          end do
107       end if
108    
109       do k1 = 1, kz
110          do k2 = 1, kz
111             pb_vert_reg(i,k1,k2) = lat_wgt * pb_vert_reg_temp(i_south+1,k1,k2) + &
112                (1.0 - lat_wgt) * pb_vert_reg_temp(i_south,k1,k2)
113          end do
114       end do     
115    end do
117    if (print_detail_regression) then
118       call da_array_print(1, pb_vert_reg_stats(1,:,:), 'pb_vert_reg_stats(1,:,:)')
119       call da_array_print(1, pb_vert_reg(1,:,:),       'pb_vert_reg(1,:,:)')
120       call da_array_print(1, pb_vert_reg_stats(:,1,:), 'pb_vert_reg_stats(:,1,:)')
121       call da_array_print(1, pb_vert_reg(:,1,:),       'pb_vert_reg(:,1,:)')
122    end if
124    if (trace_use) call da_trace_exit("da_interpolate_regcoeff")
125    
126 end subroutine da_interpolate_regcoeff