var/da/da_vtox_transforms/da_vertical_transform.inc

   1 subroutine da_vertical_transform(grid, string, be, vertical_wgt, vv, vp)
   2
   3    !---------------------------------------------------------------------
   4    ! Purpose: perform vertical transform Uv using eigenvector/eigenvalue
   5    !           of vertical covariance
   6    !
   7    ! Zhiquan (Jake) Liu, NCAR/MMM, 2015-09
   8    !    1. add appropriate comments on transform and variables
   9    !    2. replace inverse transform da_transform_vptovv
  10    !                              by da_transform_vvtovp_inv
  11    !---------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type (domain),    intent(in)    :: grid
  16    character(len=*), intent(in)    :: string      ! Character operation
  17    type (be_type),   intent(in)    :: be          ! Background error structure.
  18    real,             intent(in)    :: vertical_wgt(ims:ime,jms:jme,kms:kme) ! Weighting.
  19    type (vp_type),   intent(inout) :: vv          ! CV in gridpt/EOF space.
  20    type (vp_type),   intent(inout) :: vp          ! CV in gridpt/level space.
  21
  22    integer                         :: j, m, n     ! Loop counters.
  23    real                            :: alpha_stddev_inv ! 1/ sigma_alpha
  24    real                            :: size_inv    ! 1 / size.
  25    real                            :: alpha_me, alpha_ms, alpha_sd ! Alpha statistics.
  26
  27    if (trace_use) call da_trace_entry("da_vertical_transform")
  28
  29    select case(string)
  30
  31    case ('u');
  32
  33       !-------------------------------------------------------------------
  34       ! [1.0] Perform vp(i,j,k) = E L^{1/2} vv(i,j,m) transform:
  35       !-------------------------------------------------------------------
  36
  37       if (be % v1 % mz > 0) then
  38          call da_transform_vvtovp (grid, be % v1 % evec, be % v1 % val, vertical_wgt, &
  39             vv % v1, vp % v1, be % v1 % mz, kte) ! psi (stream function) or u (if cv7)
  40       else
  41          vp % v1(its:ite,jts:jte,kts:kte) = 0.0
  42       end if
  43
  44       if (be % v2 % mz > 0) then
  45          call da_transform_vvtovp (grid, be % v2 % evec, be % v2 % val, vertical_wgt, &
  46             vv % v2, vp % v2, be % v2 % mz, kte) ! chi_u (unbalanced chi) or v (if cv7)
  47       else
  48          vp % v2(its:ite,jts:jte,kts:kte) = 0.0
  49       end if
  50
  51       if (be % v3 % mz > 0) then
  52          call da_transform_vvtovp (grid, be % v3 % evec, be % v3 % val, vertical_wgt, &
  53             vv % v3, vp % v3, be % v3 % mz, kte) ! T_u (unbalanced T) or T (if cv7)
  54       else
  55          vp % v3(its:ite,jts:jte,kts:kte) = 0.0
  56       end if
  57
  58       if (be % v4 % mz > 0) then
  59          call da_transform_vvtovp (grid, be % v4 % evec, be % v4 % val, vertical_wgt, &
  60             vv % v4, vp % v4, be % v4 % mz, kte) ! pseudo rh=q/qs(background)
  61       else
  62          vp % v4(its:ite,jts:jte,kts:kte) = 0.0
  63       end if
  64
  65       if (be % v5 % mz > 0) then
  66          if (global) then
  67             vp % v5(its:ite,jts:jte,1) = vv % v5(its:ite,jts:jte,1)
  68          else
  69          call da_transform_vvtovp (grid, be % v5 % evec, be % v5 % val, vertical_wgt, &
  70             vv % v5, vp % v5, be % v5 % mz, kts) ! Ps_u (unbalanced Ps) or Ps (if cv7)
  71          end if
  72       else
  73          vp % v5(its:ite,jts:jte,kts:kts) = 0.0
  74       end if
  75
  76       ! for cloud_cv_options<=1 and not use_cv_w
  77       vp % v6  = 0.0 ! cloud water qcw
  78       vp % v7  = 0.0 ! rain water qrain
  79       vp % v8  = 0.0 ! cloud ice qice
  80       vp % v9  = 0.0 ! snow qsnow
  81       vp % v10 = 0.0 ! qgraupel
  82       vp % v11 = 0.0 ! vertical velocity w
  83
  84       if ( cloud_cv_options == 2 ) then
  85          if (be % v6 % mz > 0) then
  86             call da_transform_vvtovp (grid, be % v6 % evec, be % v6 % val, vertical_wgt, &
  87                vv % v6, vp % v6, be % v6 % mz, kte)
  88          end if
  89
  90          if (be % v7 % mz > 0) then
  91             call da_transform_vvtovp (grid, be % v7 % evec, be % v7 % val, vertical_wgt, &
  92                vv % v7, vp % v7, be % v7 % mz, kte)
  93          end if
  94
  95          if (be % v8 % mz > 0) then
  96             call da_transform_vvtovp (grid, be % v8 % evec, be % v8 % val, vertical_wgt, &
  97                vv % v8, vp % v8, be % v8 % mz, kte)
  98          end if
  99
 100          if (be % v9 % mz > 0) then
 101             call da_transform_vvtovp (grid, be % v9 % evec, be % v9 % val, vertical_wgt, &
 102                vv % v9, vp % v9, be % v9 % mz, kte)
 103          end if
 104
 105          if (be % v10% mz > 0) then
 106             call da_transform_vvtovp (grid, be % v10 % evec, be % v10 % val, vertical_wgt, &
 107                vv % v10, vp % v10, be % v10 % mz, kte)
 108          end if
 109       else if ( cloud_cv_options == 3 ) then
 110          if ( be % v6  % mz > 0) vp % v6  = vv % v6
 111          if ( be % v7  % mz > 0) vp % v7  = vv % v7
 112          if ( be % v8  % mz > 0) vp % v8  = vv % v8
 113          if ( be % v9  % mz > 0) vp % v9  = vv % v9
 114          if ( be % v10 % mz > 0) vp % v10 = vv % v10
 115       end if
 116
 117       if ( use_cv_w ) then
 118          if (be % v11 % mz > 0) then
 119             if ( cloud_cv_options == 2 ) then
 120                call da_transform_vvtovp (grid, be % v11 % evec, be % v11 % val, vertical_wgt, &
 121                   vv % v11, vp % v11, be % v11 % mz, kte)
 122             else if ( cloud_cv_options == 3 ) then
 123                vp % v11 = vv % v11
 124             end if
 125          end if
 126       end if
 127
 128       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
 129          do n = 1, be % ne
 130             if ( anal_type_hybrid_dual_res ) then
 131                call da_transform_vvtovp_dual_res (grid%intermediate_grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 132                                        vv % alpha(:,:,:,n), vp % alpha(:,:,:,n), be % alpha % mz, kte_int)
 133             else
 134                call da_transform_vvtovp (grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 135                                        vv % alpha(:,:,:,n), vp % alpha(:,:,:,n), be % alpha % mz, kte_int)
 136             endif
 137          end do
 138
 139 !        Calculate alpha standard deviation diagnostic:
 140 !         size_inv = 1.0 / ( (ite-its+1) * ( jte-jts+1) * be % ne * be % alpha % mz )
 141 !         alpha_me = sum(vp % alpha(its:ite,jts:jte,:)) * size_inv
 142 !         alpha_ms = sum(vp % alpha(its:ite,jts:jte,:) * vp % alpha(its:ite,jts:jte,:)) * &
 143 !                    size_inv
 144 !         alpha_sd = sqrt( alpha_ms - alpha_me * alpha_me )
 145 !         write(6,'(a,f15.5)')' Alpha std. dev = ', alpha_sd
 146       end if
 147
 148    case ('u_inv');
 149
 150       !-------------------------------------------------------------------
 151       ! [2.0] Perform inverse transform: vv(i,j,m) = L^{-1/2} E^T vp(i,j,k)
 152       !-------------------------------------------------------------------
 153
 154       if (be % v1 % mz > 0) then
 155          call da_transform_vvtovp_inv (grid, be % v1 % evec, be % v1 % val, vertical_wgt, &
 156             vp % v1, vv % v1, be % v1 % mz, kte)
 157       end if
 158
 159       if (be % v2 % mz > 0) then
 160          call da_transform_vvtovp_inv (grid, be % v2 % evec, be % v2 % val, vertical_wgt, &
 161             vp % v2, vv % v2, be % v2 % mz, kte)
 162       end if
 163
 164       if (be % v3 % mz > 0) then
 165          call da_transform_vvtovp_inv (grid, be % v3 % evec, be % v3 % val, vertical_wgt, &
 166             vp % v3, vv % v3, be % v3 % mz, kte)
 167       end if
 168
 169       if (be % v4 % mz > 0) then
 170          call da_transform_vvtovp_inv (grid, be % v4 % evec, be % v4 % val, vertical_wgt, &
 171             vp % v4, vv % v4, be % v4 % mz, kte)
 172       end if
 173
 174       if (be % v5 % mz > 0) then
 175          if (global) then
 176             vv % v5(its:ite,jts:jte,1) = vp % v5(its:ite,jts:jte,1)
 177          else
 178             call da_transform_vvtovp_inv (grid, be % v5 % evec, be % v5 % val, vertical_wgt, &
 179                vp % v5, vv % v5, be % v5 % mz, kts)
 180          end if
 181       end if
 182
 183       if ( cloud_cv_options == 2 ) then
 184          if (be % v6 % mz > 0) then
 185            call da_transform_vvtovp_inv (grid, be % v6 % evec, be % v6 % val, vertical_wgt, &
 186               vp % v6, vv % v6, be % v6 % mz, kte)
 187          end if
 188
 189          if (be % v7 % mz > 0) then
 190            call da_transform_vvtovp_inv (grid, be % v7 % evec, be % v7 % val, vertical_wgt, &
 191               vp % v7, vv % v7, be % v7 % mz, kte)
 192          end if
 193
 194          if (be % v8 % mz > 0) then
 195            call da_transform_vvtovp_inv (grid, be % v8 % evec, be % v8 % val, vertical_wgt, &
 196               vp % v8, vv % v8, be % v8 % mz, kte)
 197          end if
 198
 199          if (be % v9 % mz > 0) then
 200            call da_transform_vvtovp_inv (grid, be % v9 % evec, be % v9 % val, vertical_wgt, &
 201               vp % v9, vv % v9, be % v9 % mz, kte)
 202          end if
 203
 204          if (be % v10 % mz > 0) then
 205            call da_transform_vvtovp_inv (grid, be % v10 % evec, be % v10 % val, vertical_wgt, &
 206               vp % v10, vv % v10, be % v10 % mz, kte)
 207          end if
 208
 209       else if ( cloud_cv_options == 3 ) then
 210          if (be % v6 % mz > 0) then
 211             vv % v6 = vp % v6
 212          end if
 213
 214          if (be % v7 % mz > 0) then
 215             vv % v7 = vp % v7
 216          end if
 217
 218          if (be % v8 % mz > 0) then
 219             vv % v8 = vp % v8
 220          end if
 221
 222          if (be % v9 % mz > 0) then
 223             vv % v9 = vp % v9
 224          end if
 225
 226          if (be % v10 % mz > 0) then
 227             vv % v10 = vp % v10
 228          end if
 229
 230       end if
 231
 232       if ( use_cv_w ) then
 233          if (be % v11 % mz > 0) then
 234              if ( cloud_cv_options == 2 ) then
 235                 call da_transform_vvtovp_inv (grid, be % v11 % evec, be % v11 % val, vertical_wgt, &
 236                    vp % v11, vv % v11, be % v11 % mz, kte)
 237              else if ( cloud_cv_options == 3 ) then
 238                 vv % v11 = vp % v11
 239              end if
 240          end if
 241       end if
 242
 243       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
 244          do n = 1, be % ne
 245 !           call da_transform_vptovv (be % alpha % evec, be % alpha % val, vertical_wgt, &
 246 !                                     vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kds,kde, &
 247 !                                     ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte)
 248 !            call da_transform_vptovv (be % alpha % evec, be % alpha % val, vertical_wgt, &
 249 !                                      vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kds_int,kde_int, &
 250 !                                      ims_int,ime_int, jms_int,jme_int, kms_int,kme_int, its_int,ite_int, &
 251 !                                      jts_int,jte_int, kts_int,kte_int)
 252
 253             call da_transform_vvtovp_inv (grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 254                                       vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kte)
 255
 256          end do
 257       end if
 258
 259    case ('u_adj');
 260
 261       !-------------------------------------------------------------------
 262       ! [3.0] Perform adjoint transform: vv_adj = L^{1/2} E^T vp_adj
 263       !-------------------------------------------------------------------
 264
 265       if (be % v1 % mz > 0) then
 266          call da_transform_vvtovp_adj (grid, be % v1 % evec, be % v1 % val, vertical_wgt, &
 267             vp % v1, vv % v1, be % v1 % mz, kte)
 268       end if
 269
 270       if (be % v2 % mz > 0) then
 271          call da_transform_vvtovp_adj (grid, be % v2 % evec, be % v2 % val, vertical_wgt, &
 272             vp % v2, vv % v2, be % v2 % mz, kte)
 273       end if
 274
 275       if (be % v3 % mz > 0) then
 276          call da_transform_vvtovp_adj (grid, be % v3 % evec, be % v3 % val, vertical_wgt, &
 277             vp % v3, vv % v3, be % v3 % mz, kte)
 278       end if
 279
 280       if (be % v4 % mz > 0) then
 281          call da_transform_vvtovp_adj (grid, be % v4 % evec, be % v4 % val, vertical_wgt, &
 282             vp % v4, vv % v4, be % v4 % mz, kte)
 283       end if
 284
 285       if (be % v5 % mz > 0) then
 286          if (global) then
 287             vv % v5(its:ite,jts:jte,1) = vp % v5(its:ite,jts:jte,1)
 288          else
 289             call da_transform_vvtovp_adj (grid, be % v5 % evec, be % v5 % val, vertical_wgt, &
 290                vp % v5, vv % v5, be % v5 % mz, kts)
 291          end if
 292       end if
 293
 294       if ( cloud_cv_options == 2 ) then
 295          if (be % v6 % mz > 0) then
 296             call da_transform_vvtovp_adj (grid, be % v6 % evec, be % v6 % val, vertical_wgt, &
 297                vp % v6, vv % v6, be % v6 % mz, kte)
 298          end if
 299
 300          if (be % v7 % mz > 0) then
 301             call da_transform_vvtovp_adj (grid, be % v7 % evec, be % v7 % val, vertical_wgt, &
 302                vp % v7, vv % v7, be % v7 % mz, kte)
 303          end if
 304
 305          if (be % v8 % mz > 0) then
 306             call da_transform_vvtovp_adj (grid, be % v8 % evec, be % v8 % val, vertical_wgt, &
 307                vp % v8, vv % v8, be % v8 % mz, kte)
 308          end if
 309
 310          if (be % v9 % mz > 0) then
 311             call da_transform_vvtovp_adj (grid, be % v9 % evec, be % v9 % val, vertical_wgt, &
 312                vp % v9, vv % v9, be % v9 % mz, kte)
 313          end if
 314
 315          if (be % v10 % mz > 0) then
 316             call da_transform_vvtovp_adj (grid, be % v10 % evec, be % v10 % val, vertical_wgt, &
 317                vp % v10, vv % v10, be % v10 % mz, kte)
 318          end if
 319
 320       else if ( cloud_cv_options == 3 ) then
 321          if (be % v6 % mz > 0) then
 322             vv % v6 = vp % v6
 323          end if
 324
 325          if (be % v7 % mz > 0) then
 326             vv % v7 = vp % v7
 327          end if
 328
 329          if (be % v8 % mz > 0) then
 330             vv % v8 = vp % v8
 331          end if
 332
 333          if (be % v9 % mz > 0) then
 334             vv % v9 = vp % v9
 335          end if
 336
 337          if (be % v10 % mz > 0) then
 338             vv % v10 = vp % v10
 339          end if
 340
 341       end if
 342
 343       if ( use_cv_w ) then
 344          if (be % v11 % mz > 0) then
 345             if ( cloud_cv_options == 2 ) then
 346                call da_transform_vvtovp_adj (grid, be % v11 % evec, be % v11 % val, vertical_wgt, &
 347                   vp % v11, vv % v11, be % v11 % mz, kte)
 348             else if ( cloud_cv_options == 3 ) then
 349                vv % v11 = vp % v11
 350             end if
 351          end if
 352       end if
 353
 354       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
 355          do n = 1, be % ne
 356             if ( anal_type_hybrid_dual_res ) then
 357                call da_transform_vvtovp_adj_dual_res (grid%intermediate_grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 358                                            vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kte_int)
 359             else
 360                call da_transform_vvtovp_adj (grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 361                                            vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kte_int)
 362             endif
 363          end do
 364       end if
 365
 366    case default;
 367
 368       call da_error(__FILE__,__LINE__, &
 369          (/"Invalid da_vertical_transform option "//string/))
 370
 371    end select
 372
 373    if (trace_use) call da_trace_exit("da_vertical_transform")
 374
 375 end subroutine da_vertical_transform
 376
 377