wrfv2_fire/var/gen_be/gen_be_ep2.f90

   1 program gen_be_ep2
   2 !
   3 !----------------------------------------------------------------------
   4 !  Purpose : To convert WRF ensemble to format required for use as
   5 !  flow-dependent perturbations in WRF-Var (alpha control variable,
   6 !  alphacv_method = 2).
   7 !
   8 !  Owner: Dale Barker (NCAR/MMM)
   9 !  Please acknowledge author/institute in work that uses this code.
  10 !
  11 !----------------------------------------------------------------------
  12
  13 #ifdef crayx1
  14 #define iargc ipxfargc
  15 #endif
  16
  17    use da_control, only : stderr, stdout, filename_len
  18    use da_tools_serial, only : da_get_unit, da_free_unit
  19    use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh
  20
  21    implicit none
  22
  23    character (len=filename_len)   :: directory                 ! General filename stub.
  24    character (len=filename_len)   :: filename                  ! General filename stub.
  25    character (len=filename_len)   :: input_file                ! Input file.
  26    character (len=filename_len)   :: output_file               ! Output file.
  27    character (len=10)    :: date                      ! Character date.
  28    character (len=10)    :: var                       ! Variable to search for.
  29    character (len=3)     :: cne                       ! Ensemble size.
  30    character (len=3)     :: ce                        ! Member index -> character.
  31
  32    integer, external     :: iargc
  33    integer               :: numarg
  34    integer               :: ne                        ! Ensemble size.
  35    integer               :: i, j, k, member           ! Loop counters.
  36    integer               :: dim1                      ! Dimensions of grid (T points).
  37    integer               :: dim1s                     ! Dimensions of grid (vor/psi pts).
  38    integer               :: dim2                      ! Dimensions of grid (T points).
  39    integer               :: dim2s                     ! Dimensions of grid (vor/psi pts).
  40    integer               :: dim3                      ! Dimensions of grid (T points).
  41    real                  :: member_inv                ! 1 / member.
  42    real                  :: ds                        ! Grid resolution.
  43    logical               :: remove_mean               ! Remove mean from standard fields.
  44
  45    real, allocatable     :: u(:,:,:)                  ! u-wind.
  46    real, allocatable     :: v(:,:,:)                  ! v-wind.
  47    real, allocatable     :: temp(:,:,:)               ! Temperature.
  48    real, allocatable     :: q(:,:,:)                  ! Specific humidity.
  49    real, allocatable     :: qcloud(:,:,:)             ! Cloud.
  50    real, allocatable     :: qrain(:,:,:)              ! Rain.
  51    real, allocatable     :: psfc(:,:)                 ! Surface pressure.
  52    real, allocatable     :: u_mean(:,:,:)             ! u-wind.
  53    real, allocatable     :: v_mean(:,:,:)             ! v-wind.
  54    real, allocatable     :: temp_mean(:,:,:)          ! Temperature.
  55    real, allocatable     :: q_mean(:,:,:)             ! Specific humidity.
  56    real, allocatable     :: qcloud_mean(:,:,:)        ! Cloud.
  57    real, allocatable     :: qrain_mean(:,:,:)         ! Rain.
  58    real, allocatable     :: psfc_mean(:,:)            ! Surface pressure.
  59    real, allocatable     :: u_mnsq(:,:,:)             ! u-wind.
  60    real, allocatable     :: v_mnsq(:,:,:)             ! v-wind.
  61    real, allocatable     :: temp_mnsq(:,:,:)          ! Temperature.
  62    real, allocatable     :: q_mnsq(:,:,:)             ! Specific humidity.
  63    real, allocatable     :: qcloud_mnsq(:,:,:)        ! Cloud.
  64    real, allocatable     :: qrain_mnsq(:,:,:)         ! Rain.
  65    real, allocatable     :: psfc_mnsq(:,:)            ! Surface pressure.
  66
  67    real, allocatable     :: utmp(:,:)                 ! u-wind.
  68    real, allocatable     :: vtmp(:,:)                 ! v-wind.
  69    real, allocatable     :: ttmp(:,:)                 ! temperature.
  70    real, allocatable     :: dummy(:,:)                ! dummy.
  71
  72    integer :: gen_be_iunit, gen_be_ounit
  73
  74    stderr = 0
  75    stdout = 6
  76
  77 !---------------------------------------------------------------------------------------------
  78    write(6,'(/a)')' [1] Initialize information.'
  79 !---------------------------------------------------------------------------------------------
  80
  81    call da_get_unit(gen_be_iunit)
  82    call da_get_unit(gen_be_ounit)
  83
  84    remove_mean = .true.
  85
  86    numarg = iargc()
  87    if ( numarg /= 4 )then
  88       write(UNIT=6,FMT='(a)') &
  89         "Usage: gen_be_ep2 date ne <directory> <filename> Stop"
  90       stop
  91    end if
  92
  93    ! Initialse to stop Cray compiler complaining
  94    date=""
  95    cne=""
  96    directory=""
  97    filename=""
  98
  99    call getarg( 1, date )
 100    call getarg( 2, cne )
 101    read(cne,'(i3)')ne
 102    call getarg( 3, directory )
 103    call getarg( 4, filename )
 104
 105    if ( remove_mean ) then
 106       write(6,'(a,a)')' Computing gen_be ensemble perturbation files for date ', date
 107    else
 108       write(6,'(a,a)')' Computing gen_be ensemble forecast files for date ', date
 109    end if
 110    write(6,'(a)')' Perturbations are in MODEL SPACE (u, v, t, q, qcloud, qrain, ps)'
 111    write(6,'(a,i4)')' Ensemble Size = ', ne
 112    write(6,'(a,a)')' Directory = ', trim(directory)
 113    write(6,'(a,a)')' Filename = ', trim(filename)
 114
 115 !---------------------------------------------------------------------------------------------
 116    write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
 117 !---------------------------------------------------------------------------------------------
 118
 119 !  Get grid dimensions from first T field:
 120    var = "T"
 121    input_file = trim(directory)//'.e001/'//trim(filename)
 122    call da_stage0_initialize( input_file, var, dim1, dim2, dim3, ds )
 123    dim1s = dim1+1 ! u i dimension is 1 larger.
 124    dim2s = dim2+1 ! v j dimension is 1 larger.
 125
 126 !  Allocate arrays in output fields:
 127    allocate( u(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
 128    allocate( v(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
 129    allocate( temp(1:dim1,1:dim2,1:dim3) )
 130    allocate( q(1:dim1,1:dim2,1:dim3) )
 131    allocate( qcloud(1:dim1,1:dim2,1:dim3) )
 132    allocate( qrain(1:dim1,1:dim2,1:dim3) )
 133    allocate( psfc(1:dim1,1:dim2) )
 134    allocate( u_mean(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
 135    allocate( v_mean(1:dim1,1:dim2,1:dim3) )
 136    allocate( temp_mean(1:dim1,1:dim2,1:dim3) )
 137    allocate( q_mean(1:dim1,1:dim2,1:dim3) )
 138    allocate( qcloud_mean(1:dim1,1:dim2,1:dim3) )
 139    allocate( qrain_mean(1:dim1,1:dim2,1:dim3) )
 140    allocate( psfc_mean(1:dim1,1:dim2) )
 141    allocate( u_mnsq(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
 142    allocate( v_mnsq(1:dim1,1:dim2,1:dim3) )
 143    allocate( temp_mnsq(1:dim1,1:dim2,1:dim3) )
 144    allocate( q_mnsq(1:dim1,1:dim2,1:dim3) )
 145    allocate( qcloud_mnsq(1:dim1,1:dim2,1:dim3) )
 146    allocate( qrain_mnsq(1:dim1,1:dim2,1:dim3) )
 147    allocate( psfc_mnsq(1:dim1,1:dim2) )
 148    u_mean = 0.0
 149    v_mean = 0.0
 150    temp_mean = 0.0
 151    q_mean = 0.0
 152    qcloud_mean = 0.0
 153    qrain_mean = 0.0
 154    psfc_mean = 0.0
 155    u_mnsq = 0.0
 156    v_mnsq = 0.0
 157    temp_mnsq = 0.0
 158    q_mnsq = 0.0
 159    qcloud_mnsq = 0.0
 160    qrain_mnsq = 0.0
 161    psfc_mnsq = 0.0
 162
 163 !  Temporary arrays:
 164    allocate( utmp(1:dim1s,1:dim2) ) ! u on Arakawa C-grid.
 165    allocate( vtmp(1:dim1,1:dim2s) ) ! v on Arakawa C-grid.
 166    allocate( ttmp(1:dim1,1:dim2) )
 167    allocate( dummy(1:dim1,1:dim2) )
 168
 169 !---------------------------------------------------------------------------------------------
 170    write(6,'(/a)')' [3] Extract necessary fields from input NETCDF files and output.'
 171 !---------------------------------------------------------------------------------------------
 172
 173    do member = 1, ne
 174
 175       write(UNIT=ce,FMT='(i3.3)')member
 176       input_file = trim(directory)//'.e'//trim(ce)//'/'//trim(filename)
 177
 178       do k = 1, dim3
 179
 180          ! Read u, v:
 181          var = "U"
 182          call da_get_field( input_file, var, 3, dim1s, dim2, dim3, k, utmp )
 183          var = "V"
 184          call da_get_field( input_file, var, 3, dim1, dim2s, dim3, k, vtmp )
 185
 186 !        Interpolate u to mass pts:
 187          do j = 1, dim2
 188             do i = 1, dim1
 189                u(i,j,k) = 0.5 * ( utmp(i,j) + utmp(i+1,j) )
 190                v(i,j,k) = 0.5 * ( vtmp(i,j) + vtmp(i,j+1) )
 191             end do
 192          end do
 193
 194 !        Read theta, and convert to temperature:
 195          call da_get_trh( input_file, dim1, dim2, dim3, k, ttmp, dummy )
 196          temp(:,:,k) = ttmp(:,:)
 197
 198 !        Read mixing ratio, and convert to specific humidity:
 199          var = "QVAPOR"
 200          call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, dummy )
 201          q(:,:,k) = dummy(:,:) / ( 1.0 + dummy(:,:) )
 202
 203 !        Read hydrometeors (need better method to read all e.g. qsn automatically):
 204          var = "QCLOUD"
 205          call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, dummy )
 206          qcloud(:,:,k) = dummy(:,:)
 207          var = "QRAIN"
 208          call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, dummy )
 209          qrain(:,:,k) = dummy(:,:)
 210
 211       end do
 212
 213 !     Finally, extract surface pressure:
 214       var = "PSFC"
 215       call da_get_field( input_file, var, 2, dim1, dim2, dim3, 1, psfc )
 216
 217 !     Write out ensemble forecasts for this member:
 218       output_file = 'tmp.e'//ce
 219       open (gen_be_ounit, file = output_file, form='unformatted')
 220       write(gen_be_ounit)date, dim1, dim2, dim3
 221       write(gen_be_ounit)u
 222       write(gen_be_ounit)v
 223       write(gen_be_ounit)temp
 224       write(gen_be_ounit)q
 225       write(gen_be_ounit)qcloud
 226       write(gen_be_ounit)qrain
 227       write(gen_be_ounit)psfc
 228       close(gen_be_ounit)
 229
 230 !     Calculate accumulating mean and mean square:
 231       member_inv = 1.0 / real(member)
 232       u_mean = ( real( member-1 ) * u_mean + u ) * member_inv
 233       v_mean = ( real( member-1 ) * v_mean + v ) * member_inv
 234       temp_mean = ( real( member-1 ) * temp_mean + temp ) * member_inv
 235       q_mean = ( real( member-1 ) * q_mean + q ) * member_inv
 236       qcloud_mean = ( real( member-1 ) * qcloud_mean + qcloud ) * member_inv
 237       qrain_mean = ( real( member-1 ) * qrain_mean + qrain ) * member_inv
 238       psfc_mean = ( real( member-1 ) * psfc_mean + psfc ) * member_inv
 239       u_mnsq = ( real( member-1 ) * u_mnsq + u * u ) * member_inv
 240       v_mnsq = ( real( member-1 ) * v_mnsq + v * v ) * member_inv
 241       temp_mnsq = ( real( member-1 ) * temp_mnsq + temp * temp ) * member_inv
 242       q_mnsq = ( real( member-1 ) * q_mnsq + q * q ) * member_inv
 243       qcloud_mnsq = ( real( member-1 ) * qcloud_mnsq + qcloud * qcloud ) * member_inv
 244       qrain_mnsq = ( real( member-1 ) * qrain_mnsq + qrain * qrain ) * member_inv
 245       psfc_mnsq = ( real( member-1 ) * psfc_mnsq + psfc * psfc ) * member_inv
 246
 247    end do
 248
 249    deallocate( utmp )
 250    deallocate( vtmp )
 251    deallocate( ttmp )
 252    deallocate( dummy )
 253
 254 !---------------------------------------------------------------------------------------------
 255    write(6,'(/a)')' [4] Compute perturbations and output'
 256 !---------------------------------------------------------------------------------------------
 257
 258    if ( remove_mean ) then
 259       write(6,'(a)') "    Calculate ensemble perturbations"
 260    else
 261       write(6,'(a)') "    WARNING: Not removing ensemble mean (outputs are full-fields)"
 262    end if
 263
 264    do member = 1, ne
 265       write(UNIT=ce,FMT='(i3.3)')member
 266
 267 !     Re-read ensemble member standard fields:
 268       input_file = 'tmp.e'//ce
 269       open (gen_be_iunit, file = input_file, form='unformatted')
 270       read(gen_be_iunit)date, dim1, dim2, dim3
 271       read(gen_be_iunit)u
 272       read(gen_be_iunit)v
 273       read(gen_be_iunit)temp
 274       read(gen_be_iunit)q
 275       read(gen_be_iunit)qcloud
 276       read(gen_be_iunit)qrain
 277       read(gen_be_iunit)psfc
 278       close(gen_be_iunit)
 279
 280       if ( remove_mean ) then
 281          u = u - u_mean
 282          v = v - v_mean
 283          temp = temp - temp_mean
 284          q = q - q_mean
 285          qcloud = qcloud - qcloud_mean
 286          qrain = qrain - qrain_mean
 287          psfc = psfc - psfc_mean
 288       end if
 289
 290 !     Write out perturbations for this member:
 291
 292       output_file = 'u.e'//trim(ce) ! Output u.
 293       open (gen_be_ounit, file = output_file, form='unformatted')
 294       write(gen_be_ounit)dim1, dim2, dim3
 295       write(gen_be_ounit)u
 296       close(gen_be_ounit)
 297
 298       output_file = 'v.e'//trim(ce) ! Output v.
 299       open (gen_be_ounit, file = output_file, form='unformatted')
 300       write(gen_be_ounit)dim1, dim2, dim3
 301       write(gen_be_ounit)v
 302       close(gen_be_ounit)
 303
 304       output_file = 't.e'//trim(ce) ! Output t.
 305       open (gen_be_ounit, file = output_file, form='unformatted')
 306       write(gen_be_ounit)dim1, dim2, dim3
 307       write(gen_be_ounit)temp
 308       close(gen_be_ounit)
 309
 310       output_file = 'q.e'//trim(ce) ! Output q.
 311       open (gen_be_ounit, file = output_file, form='unformatted')
 312       write(gen_be_ounit)dim1, dim2, dim3
 313       write(gen_be_ounit)q
 314       close(gen_be_ounit)
 315
 316       output_file = 'qcloud.e'//trim(ce) ! Output qcloud.
 317       open (gen_be_ounit, file = output_file, form='unformatted')
 318       write(gen_be_ounit)dim1, dim2, dim3
 319       write(gen_be_ounit)qcloud
 320       close(gen_be_ounit)
 321
 322       output_file = 'qrain.e'//trim(ce) ! Output qrain.
 323       open (gen_be_ounit, file = output_file, form='unformatted')
 324       write(gen_be_ounit)dim1, dim2, dim3
 325       write(gen_be_ounit)qrain
 326       close(gen_be_ounit)
 327
 328       output_file = 'ps.e'//trim(ce) ! Output ps.
 329       open (gen_be_ounit, file = output_file, form='unformatted')
 330       write(gen_be_ounit)dim1, dim2, dim3
 331       write(gen_be_ounit)psfc
 332       close(gen_be_ounit)
 333
 334    end do
 335
 336 !  Write out mean/stdv fields (stdv stored in mnsq arrays):
 337    u_mnsq = sqrt( u_mnsq - u_mean * u_mean )
 338    v_mnsq = sqrt( v_mnsq - v_mean * v_mean )
 339    temp_mnsq = sqrt( temp_mnsq - temp_mean * temp_mean )
 340    q_mnsq = sqrt( q_mnsq - q_mean * q_mean )
 341    qcloud_mnsq = sqrt( qcloud_mnsq - qcloud_mean * qcloud_mean )
 342    qrain_mnsq = sqrt( qrain_mnsq - qrain_mean * qrain_mean )
 343    psfc_mnsq = sqrt( psfc_mnsq - psfc_mean * psfc_mean )
 344
 345    output_file = 'u.mean' ! Output u.
 346    open (gen_be_ounit, file = output_file, form='unformatted')
 347    write(gen_be_ounit)dim1, dim2, dim3
 348    write(gen_be_ounit)u_mean
 349    close(gen_be_ounit)
 350
 351    output_file = 'u.stdv' ! Output u.
 352    open (gen_be_ounit, file = output_file, form='unformatted')
 353    write(gen_be_ounit)dim1, dim2, dim3
 354    write(gen_be_ounit)u_mnsq
 355    close(gen_be_ounit)
 356
 357    output_file = 'v.mean' ! Output v.
 358    open (gen_be_ounit, file = output_file, form='unformatted')
 359    write(gen_be_ounit)dim1, dim2, dim3
 360    write(gen_be_ounit)v_mean
 361    close(gen_be_ounit)
 362
 363    output_file = 'v.stdv' ! Output v.
 364    open (gen_be_ounit, file = output_file, form='unformatted')
 365    write(gen_be_ounit)dim1, dim2, dim3
 366    write(gen_be_ounit)v_mnsq
 367    close(gen_be_ounit)
 368
 369    output_file = 't.mean' ! Output t.
 370    open (gen_be_ounit, file = output_file, form='unformatted')
 371    write(gen_be_ounit)dim1, dim2, dim3
 372    write(gen_be_ounit)temp_mean
 373    close(gen_be_ounit)
 374
 375    output_file = 't.stdv' ! Output t.
 376    open (gen_be_ounit, file = output_file, form='unformatted')
 377    write(gen_be_ounit)dim1, dim2, dim3
 378    write(gen_be_ounit)temp_mnsq
 379    close(gen_be_ounit)
 380
 381    output_file = 'q.mean' ! Output q.
 382    open (gen_be_ounit, file = output_file, form='unformatted')
 383    write(gen_be_ounit)dim1, dim2, dim3
 384    write(gen_be_ounit)q_mean
 385    close(gen_be_ounit)
 386
 387    output_file = 'q.stdv' ! Output q.
 388    open (gen_be_ounit, file = output_file, form='unformatted')
 389    write(gen_be_ounit)dim1, dim2, dim3
 390    write(gen_be_ounit)q_mnsq
 391    close(gen_be_ounit)
 392
 393    output_file = 'qcloud.mean' ! Output qcloud.
 394    open (gen_be_ounit, file = output_file, form='unformatted')
 395    write(gen_be_ounit)dim1, dim2, dim3
 396    write(gen_be_ounit)qcloud_mean
 397    close(gen_be_ounit)
 398
 399    output_file = 'qcloud.stdv' ! Output qcloud.
 400    open (gen_be_ounit, file = output_file, form='unformatted')
 401    write(gen_be_ounit)dim1, dim2, dim3
 402    write(gen_be_ounit)qcloud_mnsq
 403    close(gen_be_ounit)
 404
 405    output_file = 'qrain.mean' ! Output qrain.
 406    open (gen_be_ounit, file = output_file, form='unformatted')
 407    write(gen_be_ounit)dim1, dim2, dim3
 408    write(gen_be_ounit)qrain_mean
 409    close(gen_be_ounit)
 410
 411    output_file = 'qrain.stdv' ! Output qrain.
 412    open (gen_be_ounit, file = output_file, form='unformatted')
 413    write(gen_be_ounit)dim1, dim2, dim3
 414    write(gen_be_ounit)qrain_mnsq
 415    close(gen_be_ounit)
 416
 417    output_file = 'ps.mean' ! Output ps.
 418    open (gen_be_ounit, file = output_file, form='unformatted')
 419    write(gen_be_ounit)dim1, dim2, dim3
 420    write(gen_be_ounit)psfc_mean
 421    close(gen_be_ounit)
 422
 423    output_file = 'ps.stdv' ! Output ps.
 424    open (gen_be_ounit, file = output_file, form='unformatted')
 425    write(gen_be_ounit)dim1, dim2, dim3
 426    write(gen_be_ounit)psfc_mnsq
 427    close(gen_be_ounit)
 428
 429    call da_free_unit(gen_be_iunit)
 430    call da_free_unit(gen_be_ounit)
 431
 432 #ifdef crayx1
 433 contains
 434
 435    subroutine getarg(i, harg)
 436      implicit none
 437      character(len=*) :: harg
 438      integer :: ierr, ilen, i
 439
 440      call pxfgetarg(i, harg, ilen, ierr)
 441      return
 442    end subroutine getarg
 443 #endif
 444
 445 end program gen_be_ep2
 446