var/da/da_setup_structures/da_setup_flow_predictors.inc

   1 subroutine da_setup_flow_predictors( ix, jy, kz, ne, ep, its, ite, jts, jte, kts, kte )
   2
   3    !------------------------------------------------------------------------------
   4    ! Purpose: Setup structures for flow-dependent information and read it in.
   5    !------------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer, intent(in)         :: ix, jy, kz            ! EP grid dimensions.
  10    integer, intent(in)         :: its, jts, kts         ! Tile start.
  11    integer, intent(in)         :: ite, jte, kte         ! Tile end.
  12    integer, intent(in)         :: ne                    ! Ensemble size.
  13    type (ep_type), intent(inout):: ep                   ! Flow-dependent info.
  14
  15    character(len=4), allocatable :: ep_dir_name(:)      ! Input ep dir name
  16    character*10                :: ce                    ! Member index -> character.
  17    character(len=filename_len) :: filename              ! Input filename.
  18    character*10                :: var(1:5)              ! Variable name.
  19    integer                     :: ni, nj, nk            ! Grid dimensions.
  20    integer                     :: e                     ! Loop counter
  21    logical                     :: ldum1, ldum2,nkdum    ! Dummy.
  22    real*4, allocatable         :: temp3d_r4(:,:,:)      ! Temporary, real*4 array
  23    real*4, allocatable         :: temp2d_r4(:,:)        ! Temporary, real*4 array
  24    real,   allocatable         :: temp3d(:,:,:)         ! Temporary array
  25    real,   allocatable         :: temp2d(:,:)           ! Temporary array
  26
  27    real                        :: ens_scaling_inv       ! Ensemble scaling of perturbations.
  28
  29    integer                     :: ep_unit,t,te,nte,it,ie
  30    integer                     :: ijk
  31    integer, parameter          :: is_r8 = 1
  32    integer, parameter          :: is_r4 = 11
  33
  34    if (trace_use) call da_trace_entry("da_setup_flow_predictors")
  35
  36    call da_get_unit(ep_unit)
  37
  38    call da_message((/"Set up flow-dependent information"/))
  39
  40    ep % ne = ne
  41    ijk = ix * jy * kz
  42
  43    ens_scaling_inv = 1.0
  44    if (ne > 1) ens_scaling_inv = 1.0 / sqrt(real(ne-1))
  45
  46    ! Decide which space we are introducing flow-dependence:
  47    if (alphacv_method == alphacv_method_xa) then    ! xa space (default)
  48       var(1) = 'u'
  49       var(2) = 'v'
  50       var(3) = 't'
  51       var(4) = 'q'
  52       var(5) = 'ps'
  53       if ( ep_format == is_r4 ) then
  54          var(1) = 'U'
  55          var(2) = 'V'
  56          var(3) = 'T'
  57          var(4) = 'QVAPOR'
  58          var(5) = 'PSFC'
  59       end if
  60    else                               ! vp space
  61       var(1) = 'psi'
  62       var(2) = 'chi_u'
  63       var(3) = 't_u'
  64       var(4) = 'rh'
  65       var(5) = 'ps_u'
  66    end if
  67
  68    if ( num_fgat_time > 1 .and. use_4denvar ) then
  69       nte = num_fgat_time * ne
  70       allocate(ep_dir_name(num_fgat_time))
  71       do it = 1, num_fgat_time
  72          write(ep_dir_name(it), fmt='(a, i2.2)') 'ep', it    !ep01,ep02,...
  73       end do
  74    else
  75       nte = ne
  76       allocate(ep_dir_name(1))
  77       ep_dir_name(1) = 'ep'
  78    end if
  79
  80    if ( rootproc .and. ep_format == is_r4 ) then
  81       allocate(temp3d_r4(1:ix,1:jy,1:kz))
  82       allocate(temp2d_r4(1:ix,1:jy))
  83    end if
  84    allocate(temp3d(1:ix,1:jy,1:kz))
  85    allocate(temp2d(1:ix,1:jy))
  86
  87    !---------------------------------------------------------------------------
  88    ! Input ensemble perturbations
  89    !---------------------------------------------------------------------------
  90    do te = 1, nte !loop over num_fgat_time*ne
  91
  92       it = te/ne + 1       ! retrieve time window index from nte
  93       ie = mod(te, ne)     ! retrieve ensemble memeber index from nte
  94       if ( ie == 0 ) then
  95          it = it - 1
  96          ie = ne
  97       end if
  98       !write(6,*) 'it, ie = ', it, ie
  99
 100       write(unit=ce,fmt='(i3.3)') ie
 101
 102       ! v1:
 103       filename = trim(ep_dir_name(it))//'/'//trim(var(1))//'.e'//trim(ce)
 104       if ( rootproc ) then
 105          open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 106          read(unit=ep_unit) ni, nj, nk
 107
 108          if (ni /= ix .or. nj /= jy .or. nk /= kz) then
 109             write(unit=message(1),fmt='(a)') &
 110                'Inconsistent grid dimensions'
 111             write(unit=message(2),fmt='(a,3i6)') &
 112                ' Grid dims for analysis grid: ', ix, jy
 113             write(unit=message(3),fmt='(a,3i6)') &
 114                ' Grid dims for flow predictors: ', ni, nj
 115             call da_warning(__FILE__,__LINE__,message(1:3))
 116          end if
 117
 118          if ( ep_format == is_r8 ) then
 119             read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 120          else if ( ep_format == is_r4 ) then
 121             read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 122             temp3d = temp3d_r4
 123          end if
 124          close(unit=ep_unit)
 125       end if
 126       call wrf_dm_bcast_real(temp3d, ijk)
 127       ep % v1(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * &
 128                                             temp3d(its:ite,jts:jte,kts:kte)
 129
 130       ! v2:
 131       filename = trim(ep_dir_name(it))//'/'//trim(var(2))//'.e'//trim(ce)
 132       if ( rootproc ) then
 133          open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 134          read(unit=ep_unit) ni, nj, nk
 135          if ( ep_format == is_r8 ) then
 136             read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 137          else if ( ep_format == is_r4 ) then
 138             read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 139             temp3d = temp3d_r4
 140          end if
 141          close(unit=ep_unit)
 142       end if
 143       call wrf_dm_bcast_real(temp3d, ijk)
 144       ep % v2(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * &
 145                                             temp3d(its:ite,jts:jte,kts:kte)
 146
 147       ! v3:
 148       filename = trim(ep_dir_name(it))//'/'//trim(var(3))//'.e'//trim(ce)
 149       if ( rootproc ) then
 150          open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 151          read(unit=ep_unit) ni, nj, nk
 152          if ( ep_format == is_r8 ) then
 153             read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 154          else if ( ep_format == is_r4 ) then
 155             read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 156             temp3d = temp3d_r4
 157          end if
 158          close(unit=ep_unit)
 159       end if
 160       call wrf_dm_bcast_real(temp3d, ijk)
 161       ep % v3(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * &
 162                                             temp3d(its:ite,jts:jte,kts:kte)
 163
 164       ! v4:
 165       filename = trim(ep_dir_name(it))//'/'//trim(var(4))//'.e'//trim(ce)
 166       if ( rootproc ) then
 167          open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 168          read(unit=ep_unit) ni, nj, nk
 169          if ( ep_format == is_r8 ) then
 170             read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 171          else if ( ep_format == is_r4 ) then
 172             read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 173             temp3d = temp3d_r4
 174          end if
 175          close(unit=ep_unit)
 176       end if
 177       call wrf_dm_bcast_real(temp3d, ijk)
 178       ep % v4(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * &
 179                                             temp3d(its:ite,jts:jte,kts:kte)
 180
 181       ! v5:
 182       filename = trim(ep_dir_name(it))//'/'//trim(var(5))//'.e'//trim(ce)
 183       if ( rootproc ) then
 184          open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 185          read(unit=ep_unit) ni, nj, nkdum
 186          if ( ep_format == is_r8 ) then
 187             read(unit=ep_unit) temp2d(1:ix,1:jy)
 188          else if ( ep_format == is_r4 ) then
 189             read(unit=ep_unit) temp2d_r4(1:ix,1:jy)
 190             temp2d = temp2d_r4
 191          end if
 192          close(unit=ep_unit)
 193       end if
 194       call wrf_dm_bcast_real(temp2d, ix*jy)
 195       ep % v5(its:ite,jts:jte,1,te) = ens_scaling_inv * temp2d(its:ite,jts:jte)
 196
 197    end do ! num_fgat_time*ne loop
 198
 199   !  Optional include hydrometeors:
 200
 201    if ( alphacv_method == alphacv_method_xa .and. alpha_hydrometeors ) then  ! xa space
 202
 203       do te = 1, nte
 204
 205          it = te/ne + 1
 206          ie = mod(te, ne)
 207          if ( ie == 0 ) then
 208             it = it - 1
 209             ie = ne
 210          end if
 211          !write(6,*) 'it, ie = ', it, ie
 212
 213          write(unit=ce,fmt='(i3.3)') ie
 214
 215          if ( ep_format == is_r8 ) then
 216             filename = trim(ep_dir_name(it))//'/'//'qcloud'//'.e'//trim(ce)
 217          else if ( ep_format == is_r4 ) then
 218             filename = trim(ep_dir_name(it))//'/'//'QCLOUD'//'.e'//trim(ce)
 219          end if
 220          if ( rootproc ) then
 221             open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 222             read(unit=ep_unit) ni, nj, nk
 223             if ( ep_format == is_r8 ) then
 224                read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 225             else if ( ep_format == is_r4 ) then
 226                read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 227                temp3d = temp3d_r4
 228             end if
 229             close(unit=ep_unit)
 230          end if
 231          call wrf_dm_bcast_real(temp3d, ijk)
 232          ep % cw(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * temp3d(its:ite,jts:jte,kts:kte)
 233
 234          if ( ep_format == is_r8 ) then
 235             filename = trim(ep_dir_name(it))//'/'//'qrain'//'.e'//trim(ce)
 236          else if ( ep_format == is_r4 ) then
 237             filename = trim(ep_dir_name(it))//'/'//'QRAIN'//'.e'//trim(ce)
 238          end if
 239          if ( rootproc ) then
 240             open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 241             read(unit=ep_unit) ni, nj, nk
 242             if ( ep_format == is_r8 ) then
 243                read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 244             else if ( ep_format == is_r4 ) then
 245                read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 246                temp3d = temp3d_r4
 247             end if
 248             close(unit=ep_unit)
 249          end if
 250          call wrf_dm_bcast_real(temp3d, ijk)
 251          ep % rn(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * temp3d(its:ite,jts:jte,kts:kte)
 252
 253          if ( ep_format == is_r8 ) then
 254             filename = trim(ep_dir_name(it))//'/'//'qice'//'.e'//trim(ce)
 255          else if ( ep_format == is_r4 ) then
 256             filename = trim(ep_dir_name(it))//'/'//'QICE'//'.e'//trim(ce)
 257          end if
 258          if ( rootproc ) then
 259             open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 260             read(unit=ep_unit) ni, nj, nk
 261             if ( ep_format == is_r8 ) then
 262                read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 263             else if ( ep_format == is_r4 ) then
 264                read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 265                temp3d = temp3d_r4
 266             end if
 267             close(unit=ep_unit)
 268          end if
 269          call wrf_dm_bcast_real(temp3d, ijk)
 270          ep % ci(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * temp3d(its:ite,jts:jte,kts:kte)
 271
 272          if ( ep_format == is_r8 ) then
 273             filename = trim(ep_dir_name(it))//'/'//'qsnow'//'.e'//trim(ce)
 274          else if ( ep_format == is_r4 ) then
 275             filename = trim(ep_dir_name(it))//'/'//'QSNOW'//'.e'//trim(ce)
 276          end if
 277          if ( rootproc ) then
 278             open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 279             read(unit=ep_unit) ni, nj, nk
 280             if ( ep_format == is_r8 ) then
 281                read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 282             else if ( ep_format == is_r4 ) then
 283                read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 284                temp3d = temp3d_r4
 285             end if
 286             close(unit=ep_unit)
 287          end if
 288          call wrf_dm_bcast_real(temp3d, ijk)
 289          ep % sn(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * temp3d(its:ite,jts:jte,kts:kte)
 290
 291          if ( ep_format == is_r8 ) then
 292             filename = trim(ep_dir_name(it))//'/'//'qgraup'//'.e'//trim(ce)
 293          else if ( ep_format == is_r4 ) then
 294             filename = trim(ep_dir_name(it))//'/'//'QGRAUP'//'.e'//trim(ce)
 295          end if
 296          if ( rootproc ) then
 297             open(unit=ep_unit, file = filename, form = 'unformatted', status = 'old')
 298             read(unit=ep_unit) ni, nj, nk
 299             if ( ep_format == is_r8 ) then
 300                read(unit=ep_unit) temp3d(1:ix,1:jy,1:kz)
 301             else if ( ep_format == is_r4 ) then
 302                read(unit=ep_unit) temp3d_r4(1:ix,1:jy,1:kz)
 303                temp3d = temp3d_r4
 304             end if
 305             close(unit=ep_unit)
 306          end if
 307          call wrf_dm_bcast_real(temp3d, ijk)
 308          ep % gr(its:ite,jts:jte,kts:kte,te) = ens_scaling_inv * temp3d(its:ite,jts:jte,kts:kte)
 309
 310       end do ! num_fgat_time*ne loop
 311
 312    end if ! hydrometeors
 313
 314    if ( rootproc .and. ep_format == is_r4 ) then
 315       deallocate(temp3d_r4)
 316       deallocate(temp2d_r4)
 317    end if
 318    deallocate(temp3d)
 319    deallocate(temp2d)
 320
 321    write(unit=message(1),fmt='(a,e12.5)') &
 322       ' Max U    ep: ', maxval(ep % v1(its:ite,jts:jte,:,:))
 323    write(unit=message(2),fmt='(a,e12.5)') &
 324       ' Max V    ep: ', maxval(ep % v2(its:ite,jts:jte,:,:))
 325    write(unit=message(3),fmt='(a,e12.5)') &
 326       ' Max T    ep: ', maxval(ep % v3(its:ite,jts:jte,:,:))
 327    write(unit=message(4),fmt='(a,e12.5)') &
 328       ' Max Q    ep: ', maxval(ep % v4(its:ite,jts:jte,:,:))
 329    write(unit=message(5),fmt='(a,e12.5)') &
 330       ' Max Psfc ep: ', maxval(ep % v5(its:ite,jts:jte,1,:))
 331    call da_message(message(1:5))
 332
 333    if ( alphacv_method == alphacv_method_xa .and. alpha_hydrometeors ) then
 334       write(unit=message(1),fmt='(a,e12.5)') &
 335          ' Max QCLOUD ep: ', maxval(ep % cw(its:ite,jts:jte,:,:))
 336       write(unit=message(2),fmt='(a,e12.5)') &
 337          ' Max QRAIN  ep: ', maxval(ep % rn(its:ite,jts:jte,:,:))
 338       write(unit=message(3),fmt='(a,e12.5)') &
 339          ' Max QICE   ep: ', maxval(ep % ci(its:ite,jts:jte,:,:))
 340       write(unit=message(4),fmt='(a,e12.5)') &
 341          ' Max QSNOW  ep: ', maxval(ep % sn(its:ite,jts:jte,:,:))
 342       write(unit=message(5),fmt='(a,e12.5)') &
 343          ' Max QGRAUP ep: ', maxval(ep % gr(its:ite,jts:jte,:,:))
 344       call da_message(message(1:5))
 345    end if
 346
 347    deallocate(ep_dir_name)
 348    call da_free_unit(ep_unit)
 349
 350    if (trace_use) call da_trace_exit("da_setup_flow_predictors")
 351
 352 end subroutine da_setup_flow_predictors
 353
 354