var/da/da_test/da_test_vxtransform.inc

   1 subroutine da_test_vxtransform( grid, xbx, be, vv, vp)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Perform inverse/adjoint tests on control variable transform.
   5    !
   6    ! Method:  1) Test inverse and adjoint of physical variable transform.
   7    !          2) Test inverse and adjoint of vertical transform.
   8    !          3) Perform adjoint test on complete transform: <x,x> = <v_adj,v>.
   9    !
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type(domain),   intent(inout) :: grid
  15    type(xbx_type), intent(in)    :: xbx   ! Header & non-gridded vars.
  16    type(be_type),  intent(in)    :: be    ! background error structure.
  17    type(vp_type),  intent(out)   :: vp    ! Test CV structure.
  18    type(vp_type),  intent(out)   :: vv    ! Test CV structure.
  19
  20    real          :: cv(1:cv_size) ! Test control variable.
  21    integer       :: i
  22    type(vp_type) :: vp1, vv1     ! Test CV structure.
  23    real          :: field(ims:ime, jms:jme)
  24    real          :: xa2_u  (ims:ime, jms:jme, kms:kme)
  25    real          :: xa2_v  (ims:ime, jms:jme, kms:kme)
  26    real          :: xa2_t  (ims:ime, jms:jme, kms:kme)
  27    real          :: xa2_p  (ims:ime, jms:jme, kms:kme)
  28    real          :: xa2_q  (ims:ime, jms:jme, kms:kme)
  29    real          :: xa2_rh (ims:ime, jms:jme, kms:kme)
  30    real          :: xa2_rho(ims:ime, jms:jme, kms:kme)
  31    real          :: xa2_qt (ims:ime, jms:jme, kms:kme)
  32    real          :: xa2_qcw(ims:ime, jms:jme, kms:kme)
  33    real          :: xa2_qrn(ims:ime, jms:jme, kms:kme)
  34    real          :: xa2_w  (ims:ime, jms:jme, kms:kme)
  35
  36    if (trace_use) call da_trace_entry("da_test_vxtransform")
  37
  38    allocate (vp1 % v1(ims:ime,jms:jme,kms:kme))
  39    allocate (vp1 % v2(ims:ime,jms:jme,kms:kme))
  40    allocate (vp1 % v3(ims:ime,jms:jme,kms:kme))
  41    allocate (vp1 % v4(ims:ime,jms:jme,kms:kme))
  42    allocate (vp1 % v5(ims:ime,jms:jme,kms:kme))
  43 #ifdef CLOUD_CV
  44    allocate (vp1 % v6(ims:ime,jms:jme,kms:kme))
  45    allocate (vp1 % v7(ims:ime,jms:jme,kms:kme))
  46    allocate (vp1 % v8(ims:ime,jms:jme,kms:kme))
  47    allocate (vp1 % v9(ims:ime,jms:jme,kms:kme))
  48    allocate (vp1 % v10(ims:ime,jms:jme,kms:kme))
  49    allocate (vp1 % v11(ims:ime,jms:jme,kms:kme))
  50 #endif
  51
  52    allocate (vv1 % v1(ims:ime,jms:jme,1:kz_vv(1)))
  53    allocate (vv1 % v2(ims:ime,jms:jte,1:kz_vv(2)))
  54    allocate (vv1 % v3(ims:ime,jms:jme,1:kz_vv(3)))
  55    allocate (vv1 % v4(ims:ime,jms:jme,1:kz_vv(4)))
  56    allocate (vv1 % v5(ims:ime,jms:jme,1:kz_vv(5)))
  57 #ifdef CLOUD_CV
  58    allocate (vv1 % v6(ims:ime,jms:jme,1:kz_vv(6)))
  59    allocate (vv1 % v7(ims:ime,jms:jme,1:kz_vv(7)))
  60    allocate (vv1 % v8(ims:ime,jms:jme,1:kz_vv(8)))
  61    allocate (vv1 % v9(ims:ime,jms:jme,1:kz_vv(9)))
  62    allocate (vv1 % v10(ims:ime,jms:jme,1:kz_vv(10)))
  63    allocate (vv1 % v11(ims:ime,jms:jme,1:kz_vv(11)))
  64 #endif
  65    call da_zero_vp_type(vp1)
  66    call da_zero_vp_type(vv1)
  67
  68    write(unit=stdout, fmt='(/a/)') 'da_testvxtransform:'
  69
  70    write(unit=stdout, fmt='(/a/)') '---------------------------------------'
  71
  72    ! Make cv all constant value 1.0
  73
  74    call random_number(cv(:))
  75    cv(:) = cv(:) - 0.5
  76    call da_zero_x(grid%xa)
  77
  78    if ( global ) then
  79       write(unit=stdout, fmt='(/a/)') ' Inverse tests not performed for global application'
  80       go to 1111
  81    end if
  82
  83    if ( cv_options /= 5 ) then
  84       write(unit=stdout, fmt='(/a/)') ' Inverse tests not performed for cv_options=',cv_options
  85       go to 1111
  86    end if
  87
  88 #ifdef DM_PARALLEL
  89    write(unit=stdout, fmt='(/a/)') ' Inverse tests will not be done as it is MPP run'
  90 #else
  91    write(unit=stdout, fmt='(/a/)') &
  92          ' Inverse tests follows:'
  93    call da_transform_vtox( grid, xbx, be, cv, vv, vp)
  94    call da_transform_xtoxa( grid )
  95
  96    ! Store grid%xa, Vv & Vp for inverse test
  97
  98    xa2_u(ims:ime,jms:jme,:) = grid%xa % u(ims:ime,jms:jme,:)
  99    xa2_v(ims:ime,jms:jme,:) = grid%xa % v(ims:ime,jms:jme,:)
 100    xa2_w(ims:ime,jms:jme,:) = grid%xa % w(ims:ime,jms:jme,:)
 101    xa2_t(ims:ime,jms:jme,:) = grid%xa % t(ims:ime,jms:jme,:)
 102    xa2_p(ims:ime,jms:jme,:) = grid%xa % p(ims:ime,jms:jme,:)
 103    xa2_q(ims:ime,jms:jme,:) = grid%xa % q(ims:ime,jms:jme,:)
 104    xa2_rho(ims:ime,jms:jme,:) = grid%xa % rho(ims:ime,jms:jme,:)
 105
 106    if ( cv_options_hum == 2 ) then
 107       xa2_rh(ims:ime,jms:jme,:) = grid%xa % rh(ims:ime,jms:jme,:)
 108    end if
 109
 110    vv1 % v1(its:ite,jts:jte,1:kz_vv(1)) = vv % v1(its:ite,jts:jte,1:kz_vv(1))
 111    vv1 % v2(its:ite,jts:jte,1:kz_vv(2)) = vv % v2(its:ite,jts:jte,1:kz_vv(2))
 112    vv1 % v3(its:ite,jts:jte,1:kz_vv(3)) = vv % v3(its:ite,jts:jte,1:kz_vv(3))
 113    vv1 % v4(its:ite,jts:jte,1:kz_vv(4)) = vv % v4(its:ite,jts:jte,1:kz_vv(4))
 114    vv1 % v5(its:ite,jts:jte,1:kz_vv(5)) = vv % v5(its:ite,jts:jte,1:kz_vv(5))
 115 #ifdef CLOUD_CV
 116    vv1 % v6(its:ite,jts:jte,1:kz_vv(6)) = vv % v6(its:ite,jts:jte,1:kz_vv(6))
 117    vv1 % v7(its:ite,jts:jte,1:kz_vv(7)) = vv % v7(its:ite,jts:jte,1:kz_vv(7))
 118    vv1 % v8(its:ite,jts:jte,1:kz_vv(8)) = vv % v8(its:ite,jts:jte,1:kz_vv(8))
 119    vv1 % v9(its:ite,jts:jte,1:kz_vv(9)) = vv % v9(its:ite,jts:jte,1:kz_vv(9))
 120    vv1 % v10(its:ite,jts:jte,1:kz_vv(10)) = vv % v10(its:ite,jts:jte,1:kz_vv(10))
 121    vv1 % v11(its:ite,jts:jte,1:kz_vv(11)) = vv % v11(its:ite,jts:jte,1:kz_vv(11))
 122 #endif
 123
 124    vp1 % v1(its:ite,jts:jte,1:kz_vp(1)) = vp % v1(its:ite,jts:jte,1:kz_vp(1))
 125    vp1 % v2(its:ite,jts:jte,1:kz_vp(2)) = vp % v2(its:ite,jts:jte,1:kz_vp(2))
 126    vp1 % v3(its:ite,jts:jte,1:kz_vp(3)) = vp % v3(its:ite,jts:jte,1:kz_vp(3))
 127    vp1 % v4(its:ite,jts:jte,1:kz_vp(4)) = vp % v4(its:ite,jts:jte,1:kz_vp(4))
 128    vp1 % v5(its:ite,jts:jte,1:kz_vp(5)) = vp % v5(its:ite,jts:jte,1:kz_vp(5))
 129 #ifdef CLOUD_CV
 130    vp1 % v6(its:ite,jts:jte,1:kz_vp(6)) = vp % v6(its:ite,jts:jte,1:kz_vp(6))
 131    vp1 % v7(its:ite,jts:jte,1:kz_vp(7)) = vp % v7(its:ite,jts:jte,1:kz_vp(7))
 132    vp1 % v8(its:ite,jts:jte,1:kz_vp(8)) = vp % v8(its:ite,jts:jte,1:kz_vp(8))
 133    vp1 % v9(its:ite,jts:jte,1:kz_vp(9)) = vp % v9(its:ite,jts:jte,1:kz_vp(9))
 134    vp1 % v10(its:ite,jts:jte,1:kz_vp(10)) = vp % v10(its:ite,jts:jte,1:kz_vp(10))
 135    vp1 % v11(its:ite,jts:jte,1:kz_vp(11)) = vp % v11(its:ite,jts:jte,1:kz_vp(11))
 136 #endif
 137    !----------------------------------------------------------------------
 138    ! [1.0]: Perform VvToVpToVv test:
 139    !----------------------------------------------------------------------
 140    ! call da_transform_xtovp( grid, xbx, vp, be)
 141
 142    if ( vert_corr == vert_corr_2 ) then
 143       ! perform vv = u_v^{-1} vp transform:
 144       call da_vertical_transform (grid, 'u_inv', be, grid%xb % vertical_inner_product, vv, vp)
 145    else
 146       vv % v1(its:ite,jts:jte,1:kz_vv(1)) = vp % v1(its:ite,jts:jte,1:kz_vp(1))
 147       vv % v2(its:ite,jts:jte,1:kz_vv(2)) = vp % v2(its:ite,jts:jte,1:kz_vp(2))
 148       vv % v3(its:ite,jts:jte,1:kz_vv(3)) = vp % v3(its:ite,jts:jte,1:kz_vp(3))
 149       vv % v4(its:ite,jts:jte,1:kz_vv(4)) = vp % v4(its:ite,jts:jte,1:kz_vp(4))
 150       vv % v5(its:ite,jts:jte,1:kz_vv(5)) = vp % v5(its:ite,jts:jte,1:kz_vp(5))
 151 #ifdef CLOUD_CV
 152       vv % v6(its:ite,jts:jte,1:kz_vv(6)) = vp % v6(its:ite,jts:jte,1:kz_vp(6))
 153       vv % v7(its:ite,jts:jte,1:kz_vv(7)) = vp % v7(its:ite,jts:jte,1:kz_vp(7))
 154       vv % v8(its:ite,jts:jte,1:kz_vv(8)) = vp % v8(its:ite,jts:jte,1:kz_vp(8))
 155       vv % v9(its:ite,jts:jte,1:kz_vv(9)) = vp % v9(its:ite,jts:jte,1:kz_vp(9))
 156       vv % v10(its:ite,jts:jte,1:kz_vv(10)) = vp % v10(its:ite,jts:jte,1:kz_vp(10))
 157       vv % v11(its:ite,jts:jte,1:kz_vv(11)) = vp % v11(its:ite,jts:jte,1:kz_vp(11))
 158 #endif
 159    end if
 160
 161    write(unit=stdout, fmt='(/a/)') 'da_check_vvtovptovv_errors'
 162
 163    call da_check_vp_errors (vv1, vv, kz_vv, its,ite, jts,jte, kts,kte)
 164
 165    !----------------------------------------------------------------------
 166    ! [2.0]: Perform VpToVvToVp test:
 167    !----------------------------------------------------------------------
 168    if ( vert_corr == vert_corr_2 ) then
 169       ! perform vp = u_v (vv) transform:
 170       call da_vertical_transform (grid, 'u', be, grid%xb % vertical_inner_product, vv, vp)
 171    else
 172       vp % v1(its:ite,jts:jte,1:kz_vv(1)) = vv % v1(its:ite,jts:jte,1:kz_vv(1))
 173       vp % v2(its:ite,jts:jte,1:kz_vv(2)) = vv % v2(its:ite,jts:jte,1:kz_vv(2))
 174       vp % v3(its:ite,jts:jte,1:kz_vv(3)) = vv % v3(its:ite,jts:jte,1:kz_vv(3))
 175       vp % v4(its:ite,jts:jte,1:kz_vv(4)) = vv % v4(its:ite,jts:jte,1:kz_vv(4))
 176       vp % v5(its:ite,jts:jte,1:kz_vv(5)) = vv % v5(its:ite,jts:jte,1:kz_vv(5))
 177 #ifdef CLOUD_CV
 178       vp % v6(its:ite,jts:jte,1:kz_vv(6)) = vv % v6(its:ite,jts:jte,1:kz_vv(6))
 179       vp % v7(its:ite,jts:jte,1:kz_vv(7)) = vv % v7(its:ite,jts:jte,1:kz_vv(7))
 180       vp % v8(its:ite,jts:jte,1:kz_vv(8)) = vv % v8(its:ite,jts:jte,1:kz_vv(8))
 181       vp % v9(its:ite,jts:jte,1:kz_vv(9)) = vv % v9(its:ite,jts:jte,1:kz_vv(9))
 182       vp % v10(its:ite,jts:jte,1:kz_vv(10)) = vv % v10(its:ite,jts:jte,1:kz_vv(10))
 183       vp % v11(its:ite,jts:jte,1:kz_vv(11)) = vv % v11(its:ite,jts:jte,1:kz_vv(11))
 184 #endif
 185    end if
 186
 187    ! Check inverse errors:
 188
 189    write(unit=stdout, fmt='(/a/)') 'da_check_vptovvtovp_errors'
 190
 191    call da_check_vp_errors( vp1, vp, kz_vv, its,ite, jts,jte, kts,kte )
 192
 193    !----------------------------------------------------------------------
 194    ! [3.0] Check_CvToVv_Adjoint:
 195    !----------------------------------------------------------------------
 196
 197    call da_check_cvtovv_adjoint( grid, xbx, be, cv, vv)
 198
 199    !----------------------------------------------------------------------
 200    ! [4.0] Test inverse physical variable transform:
 201    ! Note: Currently these test are developed for regional only
 202    !----------------------------------------------------------------------
 203
 204    if (.not. global) then
 205       call da_zero_x(grid%xa)
 206       call da_transform_vptox( grid, vp, be)
 207       ! [4.1] Test XToVpToX differences:
 208
 209       write(unit=stdout, fmt='(/a/)') &
 210          'da_check_xtovptox_errors'
 211
 212       call da_check_xtovptox_errors( grid%xa, xa2_u, xa2_v, xa2_w, xa2_t, &
 213          xa2_p, xa2_q,  xa2_rho, xa2_qt, xa2_qcw, xa2_qrn)
 214
 215       ! [4.2] Perform v_{p} = U_{p}^{-1} x transform (again):
 216
 217       call da_transform_xtovp( grid, xbx, vv, be)
 218
 219       ! [2.4] Check inverse errors:
 220
 221       write(unit=stdout, fmt='(/a/)') 'da_check_vptoxtovp_errors'
 222
 223       call da_check_vp_errors( vp, vv, kz_vp, its,ite, jts,jte, kts,kte )
 224    end if
 225
 226    !----------------------------------------------------------------------
 227    ! [5.0] Perform Vv -> Vp adjoint test:
 228    !----------------------------------------------------------------------
 229
 230    call da_check_vvtovp_adjoint( grid, be, vp, vv)
 231
 232    !----------------------------------------------------------------------
 233    ! [6.0] Perform Vp -> X  adjoint tests:
 234    !----------------------------------------------------------------------
 235
 236    call da_check_vptox_adjoint( grid, be, vp)
 237
 238 #endif
 239 1111 continue
 240
 241    !----------------------------------------------------------------------
 242    ! [7.0]: Perform adjoint test on complete transform: <x,x> = <v_adj,v>
 243    !----------------------------------------------------------------------
 244
 245    call da_check_vtox_adjoint( grid, xbx, be, cv, vv, vp)
 246
 247    !----------------------------------------------------------------------
 248    ! [8.0]: Perform Spectral transform tests for Global
 249    !----------------------------------------------------------------------
 250
 251    if (global) then
 252 #ifdef DM_PARALLEL
 253       write (unit=stdout,fmt=*)  &
 254          ' Inverse tests for spectral transforms will not be done as it is MPP run'
 255 #else
 256       call random_number(field(:,:) )
 257       field(:,:) = field(:,:) - 0.5
 258       call da_test_spectral(xbx, field)
 259 #endif
 260    end if
 261
 262    deallocate (vp1 % v1)
 263    deallocate (vp1 % v2)
 264    deallocate (vp1 % v3)
 265    deallocate (vp1 % v4)
 266    deallocate (vp1 % v5)
 267 #ifdef CLOUD_CV
 268    deallocate (vp1 % v6)
 269    deallocate (vp1 % v7)
 270    deallocate (vp1 % v8)
 271    deallocate (vp1 % v9)
 272    deallocate (vp1 % v10)
 273    deallocate (vp1 % v11)
 274 #endif
 275
 276    deallocate (vv1 % v1)
 277    deallocate (vv1 % v2)
 278    deallocate (vv1 % v3)
 279    deallocate (vv1 % v4)
 280    deallocate (vv1 % v5)
 281 #ifdef CLOUD_CV
 282    deallocate (vv1 % v6)
 283    deallocate (vv1 % v7)
 284    deallocate (vv1 % v8)
 285    deallocate (vv1 % v9)
 286    deallocate (vv1 % v10)
 287    deallocate (vv1 % v11)
 288 #endif
 289    if (trace_use) call da_trace_exit("da_test_vxtransform")
 290
 291 end subroutine da_test_vxtransform
 292
 293