var/da/da_test/da_transform_xtovp.inc

   1 subroutine da_transform_xtovp(grid, xb, xbx, xa, vp, be)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Transforms analysis to control variables (Vp-type)
   5    !    Updated for Analysis on Arakawa-C grid
   6    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   7    !
   8    ! Updates:
   9    !
  10    !       Implementation of multi-variate BE
  11    !       Syed RH Rizvi,  MMM/NESL/NCAR,  Date: 02/01/2010
  12    !---------------------------------------------------------------------------
  13
  14    implicit none
  15
  16    type(domain),            intent(inout) :: grid
  17    type(xb_type),           intent(in)    :: xb         ! First guess structure.
  18    type(xbx_type),          intent(in)    :: xbx        ! Header/non-gridded vars.
  19    type(x_type),            intent(inout) :: xa         ! Analysis increments.
  20    type(vp_type),           intent(out)   :: vp         ! CV on grid structure.
  21    type(be_type), optional, intent(in)    :: be         ! Background errors.
  22
  23    real    :: vor(ims:ime,jms:jme,kms:kme) ! Vorticity.
  24    real    :: div(ims:ime,jms:jme,kms:kme) ! Divergence.
  25
  26    real    :: one_over_m2(ims:ime,jms:jme) !   Multiplicative coeff.
  27
  28    integer :: i, j, k , k1                 ! Loop counters.
  29
  30    if (trace_use) call da_trace_entry("da_transform_xtovp")
  31
  32    if ( (cv_options == 3) .or. (cv_options == 7) ) then
  33       write(unit=message(1),fmt='(A,I3)') 'Cannot perform transform_xtovp for cv_options == ',cv_options
  34       call da_error(__FILE__,__LINE__,message(1:1))
  35    endif
  36
  37    !----------------------------------------------------------------
  38    ! [1.0] Perform transform v = U^{-1} x
  39    !----------------------------------------------------------------
  40
  41    call da_zero_vp_type (vp)
  42
  43    ! [2.2] Transform u, v to streamfunction via vorticity:
  44
  45 #ifdef A2C
  46   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  47        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  48      ipe = ipe + 1
  49      ide = ide + 1
  50   end if
  51
  52   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  53        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  54      jpe = jpe + 1
  55      jde = jde + 1
  56   end if
  57 #endif
  58
  59 #ifdef DM_PARALLEL
  60 #include "HALO_PSICHI_UV_ADJ.inc"
  61 #endif
  62
  63
  64 #ifdef A2C
  65   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  66        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  67      ipe = ipe - 1
  68      ide = ide - 1
  69   end if
  70
  71   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  72        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  73      jpe = jpe - 1
  74      jde = jde - 1
  75   end if
  76 #endif
  77
  78    call da_uv_to_vorticity(xb, xa % u, xa % v, vor)
  79
  80    ! Convert vorticity to Del**2 psi:
  81    if (.not. global) then
  82       if (fg_format == fg_format_wrf_arw_regional) then
  83          one_over_m2(its:ite,jts:jte) = 1.0 / (xb % map_factor(its:ite,jts:jte) * &
  84                                         xb % map_factor(its:ite,jts:jte))
  85          do k = kts, kte
  86             vor(its:ite,jts:jte,k) = &
  87               one_over_m2(its:ite,jts:jte)*vor(its:ite,jts:jte,k)
  88          end do
  89       else if (fg_format == fg_format_wrf_nmm_regional) then
  90          write(unit=message(1),fmt='(A,I5)') &
  91          "Needs to be developed for fg_format_nmm_regional = ",fg_format
  92          call da_error(__FILE__,__LINE__,message(1:1))
  93       else
  94          write(unit=message(1),fmt='(A,I5,A,L10)') &
  95             ' Wrong choice of fg_format= ',fg_format,' with global = ',global
  96          call da_error(__FILE__,__LINE__,message(1:1))
  97       end if
  98    end if
  99
 100    ! Calculate psi:
 101    write (unit=stdout,fmt=*) ' calling Solve_PoissonEquation for Psi'
 102    call da_solve_poissoneqn_fct(grid,xbx, vor, vp%v1)
 103
 104    ! [2.3] Transform u, v to velocity potential via divergence:
 105
 106    call da_message((/'calling UV_To_Divergence'/))
 107    call da_uv_to_divergence(xb, xa % u, xa % v, div)
 108
 109    ! Convert divergence to Del**2 chi:
 110    if (.not. global)  then
 111       if (fg_format == fg_format_wrf_arw_regional) then
 112          do k = kts, kte
 113             div(its:ite,jts:jte,k) = &
 114                one_over_m2(its:ite,jts:jte) * div(its:ite,jts:jte,k)
 115          end do
 116       else if (fg_format == fg_format_wrf_nmm_regional) then
 117          write(unit=message(1),fmt='(A,I5)') &
 118          "Needs to be developed for fg_format_nmm_regional = ",fg_format
 119          call da_error(__FILE__,__LINE__,message(1:1))
 120       else
 121          write(unit=message(1),fmt='(A,I5,A,L10)') &
 122             ' Wrong choice of fg_format= ',fg_format,' with global = ',global
 123          call da_error(__FILE__,__LINE__,message(1:1))
 124       end if
 125    end if
 126
 127    ! Calculate chi:
 128
 129    call da_message((/' calling Solve_PoissonEquation for Chi'/))
 130    call da_solve_poissoneqn_fct(grid,xbx, div, vp%v2)
 131
 132    ! [2.4] Transform chi to chi_u:
 133    call da_message((/' calculating chi_u'/))
 134    do k=kts,kte
 135       do j=jts,jte
 136          vp%v2(its:ite,j,k) = vp%v2(its:ite,j,k) - &
 137             be%reg_psi_chi(j,k)*vp%v1(its:ite,j,k)
 138       end do
 139    end do
 140
 141    ! [2.5] Compute t_u:
 142    call da_message((/' calculating t_u'/))
 143    do k1=kts,kte
 144     do k=kts,kte
 145       do j=jts,jte
 146             vp%v3(its:ite,j,k) = vp%v3(its:ite,j,k) - be%reg_psi_t(j,k,k1)*vp%v1(its:ite,j,k1)
 147       end do
 148     end do
 149    end do
 150    if ( cv_options == 6 ) then
 151       do k1=kts,kte
 152          do k=kts,kte
 153             do j=jts,jte
 154                vp%v3(its:ite,j,k) = vp%v3(its:ite,j,k) + be%reg_chi_u_t(j,k,k1)*vp%v2(its:ite,j,k1)
 155             end do
 156          end do
 157       end do
 158    end if
 159
 160    ! [2.6] Choice of moisture control variable:
 161    call da_message((/' calculating psudo rh'/))
 162    vp % v4(its:ite,jts:jte,kts:kte) = xa % q  (its:ite,jts:jte,kts:kte) /   &
 163       xb % qs (its:ite,jts:jte,kts:kte)
 164
 165    if ( cv_options == 6 ) then
 166       do k1 = kts, kte
 167          do k = kts, kte
 168             do j = jts, jte
 169                vp%v4(its:ite,j,k) = vp%v4(its:ite,j,k)    - &
 170                                     be%reg_psi_rh(j,k,k1)*vp%v1(its:ite,j,k1)    + &
 171                                     be%reg_chi_u_rh(j,k,k1)*vp%v2(its:ite,j,k1)  + &
 172                                     be%reg_t_u_rh(j,k,k1)*vp%v3(its:ite,j,k1) + &
 173                                     be%reg_ps_u_rh(j,k1)*vp%v5(its:ite,j,1)
 174             end do
 175          end do
 176       end do
 177    end if
 178
 179    ! [2.7] compute psfc_u:
 180    call da_message((/' calculating psfc_u '/))
 181    do j=jts,jte
 182       do i=its,ite
 183          vp % v5(i,j,1) = xa%psfc(i,j) - be%reg_psi_ps(j,k)*vp%v1(i,j,k)
 184       end do
 185    end do
 186    if ( cv_options == 6 ) then
 187       do j=jts,jte
 188          do i=its,ite
 189             vp % v5(i,j,1) = xa%psfc(i,j) + be%reg_chi_u_ps(j,k)*vp%v2(i,j,k)
 190          end do
 191       end do
 192    end if
 193
 194    if (trace_use) call da_trace_exit("da_transform_xtovp")
 195
 196 end subroutine da_transform_xtovp
 197
 198