da_transform_vtox.inc

References to this file elsewhere.
1 subroutine da_transform_vtox(grid, cv_size, xbx, be, ep, cv, vv, vp)
2 
3    !--------------------------------------------------------------------------
4    ! Purpose: Control variable transform x' = Uv. 
5    !--------------------------------------------------------------------------
6 
7    implicit none
8 
9    type(domain),   intent(inout) :: grid
10    integer,        intent(in)    :: cv_size ! Size of cv array.
11    type(xbx_type), intent(in)    :: xbx  ! For header & non-grid arrays.
12    type(be_type),  intent(in)    :: be   ! background errors.
13    type(ep_type),  intent(in)    :: ep   ! Ensemble perturbations.
14    real,           intent(in)    :: cv(1:cv_size)   ! control variables.
15    type(vp_type),  intent(out)   :: vv   ! grdipt/eof cv (local).
16    type(vp_type),  intent(inout) :: vp   ! grdipt/level cv (local).
17 
18    integer :: i, j, k
19    real    :: sdmd, s1md, mu
20    real    :: p(kms:kme)
21    real    :: mr_a(kms:kme)
22    real    :: mr_b(kms:kme)
23    real    :: PU, PD, coeff
24 
25    if (trace_use) call da_trace_entry("da_transform_vtox")
26 
27    call da_zero_x (grid%xa)
28 
29    !----------------------------------------------------------------------
30    ! [1.0]: Perform vv = u_h cv transform:
31    !----------------------------------------------------------------------
32 
33    if (global) then
34       call da_transform_vtovv_global(cv_size, xbx, be, cv, vv) 
35    else
36       call da_transform_vtovv(grid, cv_size, be, cv, vv)
37    end if
38 
39    !----------------------------------------------------------------------
40    ! [2.0]: Perform vp = u_v vv transform:
41    !----------------------------------------------------------------------
42 
43    if (vert_corr == vert_corr_2) then      
44       call da_vertical_transform('u', be, grid%xb % vertical_inner_product, vv, vp)
45    else
46       vp % v1(its:ite,jts:jte,kts:kte) = vv % v1(its:ite,jts:jte,kts:kte)
47       vp % v2(its:ite,jts:jte,kts:kte) = vv % v2(its:ite,jts:jte,kts:kte)
48       vp % v3(its:ite,jts:jte,kts:kte) = vv % v3(its:ite,jts:jte,kts:kte)
49       vp % v4(its:ite,jts:jte,kts:kte) = vv % v4(its:ite,jts:jte,kts:kte)
50       vp % v5(its:ite,jts:jte,kts:kte) = vv % v5(its:ite,jts:jte,kts:kte)
51    end if
52 
53    ! WHY?
54    ! if (be % ne > 0) then
55    !    vp % alpha(its:ite,jts:jte,1:be%ne) = vv%alpha(its:ite,jts:jte,1:be%ne)
56    ! end if
57 
58    !----------------------------------------------------------------------  
59    ! [3.0]: Perform x = u_p vp transform::
60    !----------------------------------------------------------------------
61 
62    call da_transform_vptox(grid, vp, be, ep)
63 
64    !----------------------------------------------------------------------
65    ! [4.0]: Move the following:
66    !----------------------------------------------------------------------
67 
68    do j=jts,jte
69       do i=its,ite
70          if (fg_format == fg_format_wrf) then
71             sdmd=0.0
72             s1md=0.0
73             do k=kts,kte
74                mr_a(k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
75                mr_b(k) = grid%xb%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))
76 
77                sdmd=sdmd+mr_a(k)*grid%xb%dnw(k)
78                s1md=s1md+(1.0+mr_b(k))*grid%xb%dnw(k)
79             end do
80 
81             mu=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md
82 
83             p(kte+1)=0.0
84 
85             do k=kte,kts,-1
86                p(k)=p(k+1)-(mu*(1.0+mr_b(k)) &
87                        + grid%xb%psac(i,j)*mr_a(k))*grid%xb%dnw(k)
88 
89                grid%xa%p(i,j,k)=0.5*(p(k)+p(k+1))
90             end do
91          else if (fg_format == fg_format_kma_global) then
92             do k=kts,kte
93                if (k == kte) then
94                   coeff=grid%xb%KMA_B(K)/(grid%xb%KMA_A(K)+grid%xb%KMA_B(K)*grid%xb%psfc(I,J)/100.0)
95                else
96                   PU = grid%xb%KMA_A(K+1) + grid%xb%KMA_B(K+1)*grid%xb%psfc(I,J)/100.0
97                   PD = grid%xb%KMA_A(K ) + grid%xb%KMA_B(K )*grid%xb%psfc(I,J)/100.0
98                   coeff=grid%xb%KMA_B(K)  *1.0/(PD-PU)**2*(-PU*(LOG(PD)-LOG(PU)) &
99                     + PD-PU)&
100                     + grid%xb%KMA_B(K+1)*1.0/(PD-PU)**2*(PD*(LOG(PD)-LOG(PU))-PD+PU)
101                end if
102                ! Here since grid%xa%psfc holds value in Pa. dlnp -> dp
103                grid%xa%p(i,j,k) =  grid%xb%p(i,j,k) * grid%xa%psfc(I,J)/100.0 * coeff
104 
105             end do
106          end if
107       end do
108    end do
109 
110    call da_pt_to_rho_lin(grid)
111  
112    ! If test_wrfvar = .true., not "XToY" transform needed to do here:
113 
114    if (.not.test_wrfvar) then
115       ! Exchange grid%xa halo region.
116 #ifdef DM_PARALLEL
117 #include "HALO_XA.inc"
118 #endif
119 
120       if (sfc_assi_options == 2) then
121          call da_transform_xtowtq (grid)
122          ! Exchange grid%xa (surface variable) halo region.
123 #ifdef DM_PARALLEL
124 #include "HALO_SFC_XA.inc"
125 #endif
126       end if
127 
128       if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
129           use_ssmitbobs .or. use_ssmiretrievalobs .or. use_gpsrefobs) then
130 
131          ! Now do something for PW
132          call da_transform_xtotpw(grid)
133 
134          ! GPS Refractivity:
135          if (use_gpsrefobs) &
136             call da_transform_xtogpsref_lin(grid)
137 
138          if (use_ssmt1obs .or. use_ssmt2obs .or. &
139               use_ssmitbobs .or. use_ssmiretrievalobs) then
140             if (global) then
141                call da_error(__FILE__,__LINE__, &
142                   (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
143             end if
144             call da_transform_xtoseasfcwind_lin(grid)
145          end if
146          if (use_ssmitbobs) call da_transform_xtotb_lin (grid)
147 
148          ! Exchange grid%xa halo region.
149 #ifdef DM_PARALLEL
150 #include "HALO_SSMI_XA.inc"
151 #endif
152       end if
153    end if
154 
155    ! Compute w increments using Richardson's eqn.
156 
157    if (use_radarobs) then
158       call da_uvprho_to_w_lin(grid)
159 
160       do k=kts,kte
161          do j=jts,jte
162             do i=its,ite
163                grid%xa%wh(i,j,k)=0.5*(grid%xa%w(i,j,k)+grid%xa%w(i,j,k+1))
164             end do
165          end do
166       end do
167 
168 #ifdef DM_PARALLEL
169 #include "HALO_RADAR_XA_W.inc"
170 #endif
171  
172       if (use_radarobs .and. use_radar_rf) then
173          ! Partition of hydrometeor increments via warm rain process
174          call da_moist_phys_lin(grid)
175       end if
176    end if
177 
178    !---------------------------------------------------------------
179    ! Polar treatment for Global 
180    !---------------------------------------------------------------
181 
182    if (global)  then   
183       call da_get_vpoles(grid%xa%u,grid%xa%v, &
184          ids, ide, jds, jde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
185       call da_get_spoles(grid%xa%t, &
186          ids, ide, jds, jde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
187       call da_get_spoles(grid%xa%p, &
188          ids, ide, jds, jde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
189       call da_get_spoles(grid%xa%q, &
190          ids, ide, jds, jde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
191       call da_get_spoles(grid%xa%psfc, &
192          ids, ide, jds, jde, ims, ime, jms, jme,   1,   1, its, ite, jts, jte,   1,   1)
193       call da_set_boundary_xa(grid)
194    end if   
195 
196    if (trace_use) call da_trace_exit("da_transform_vtox")
197 
198 end subroutine da_transform_vtox
199 
200