var/da/da_setup_structures/da_setup_be_ncep_gfs.inc

   1 subroutine da_setup_be_ncep_gfs( grid, be )
   2 !------------------------------------------------------------------------------
   3 !  PURPOSE: Define and allocate components of background errors for cv_option 3.
   4 !
   5 !  METHOD:  Allocate components in turn.
   6 !
   7 !  HISTORY: 08/02/2002 - Creation of F90 version.           Wan-Shu Wu
   8 !
   9 !  PARENT_MODULE: DA_Setup_Structures
  10 !
  11 !  Modified by Yong-Run Guo,  03/07/2009, for WRFVar 3.1
  12 !
  13 !------------------------------------------------------------------------------
  14
  15    IMPLICIT NONE
  16
  17    type (domain), intent(in)     :: grid
  18    TYPE (be_type), INTENT(INOUT) :: be                    ! Back. errors structure.
  19
  20    INTEGER                     :: ij,ijk                ! Scalar.
  21    INTEGER                     :: i, j, k, ic, jc, ii   ! Loop counters.
  22    INTEGER                     :: ier, be_unit          ! error index
  23
  24
  25 ! added for AVN
  26    integer                     :: nlath
  27    integer                     :: nsig
  28    integer                     :: m,n,m1,n1,n4            ! loop counter
  29    integer                     :: msig,mlath,nmdszh,kcap  ! dummy variables
  30    REAL, ALLOCATABLE           :: corz_kz(:,:)
  31    REAL, ALLOCATABLE           :: cord_kz(:,:)
  32    REAL, ALLOCATABLE           :: corh_kz(:,:)
  33    REAL, ALLOCATABLE           :: corq_kz(:,:)
  34    REAL, ALLOCATABLE           :: corz_avn(:,:)
  35    REAL, ALLOCATABLE           :: cord_avn(:,:)
  36    REAL, ALLOCATABLE           :: corh_avn(:,:)
  37    REAL, ALLOCATABLE           :: corq_avn(:,:)
  38    REAL, ALLOCATABLE           :: corp_avn(:)
  39    REAL, ALLOCATABLE           :: clat_avn(:),sigma_avn(:)
  40    REAL, ALLOCATABLE           :: hwll_avn(:,:,:),hwllp_avn(:),hwll_kz(:,:,:)
  41    REAL, ALLOCATABLE           :: vztdq_avn(:,:,:),vztdq_kz(:,:,:)
  42
  43    REAL, ALLOCATABLE           :: agv_avn(:,:,:),agv_kz(:,:,:)
  44    REAL, ALLOCATABLE           :: bv_avn(:,:),wgv_avn(:,:),bv_kz(:,:),wgv_kz(:,:)
  45    REAL, ALLOCATABLE           :: dsh(:),turn(:,:)
  46
  47    REAL, DIMENSION( kts: kte,  kts: kte) :: vv
  48    REAL, DIMENSION( ids: ide,  jds: jde,  kds: kde, 1:4) :: hwll
  49    REAL, DIMENSION( ids: ide,  jds: jde)                 :: hwllp, &
  50                                                             coef1, coef2, &
  51                                                             global_lat, global_fac
  52    INTEGER, DIMENSION( ids: ide,  jds: jde)              :: mlat
  53
  54    INTEGER :: nta,ndeg,ku,kz
  55    real    :: samp,s2u,tin,as(5),slim
  56    character(len=256) :: mesg
  57
  58    integer                     :: be_rf_unit, be_print_unit, it
  59    real                        :: xsum
  60    real, dimension(1)          :: xxsum
  61    real, dimension(ids: ide,  jds: jde)                   :: global_2d
  62    real, dimension(ids: ide,  jds: jde,  kds: kde)        :: global_3d
  63    real, dimension(ims: ime,  jms: jme,  kms: kme)        :: corz_3d, vz_3d
  64    real, dimension(ims: ime,  jms: jme)                   :: corp_2d
  65    character(len=6)            :: vname
  66
  67    write (6,'(A)') ' ----------------------------------------------------------'
  68    write (6,'(A,I3)') ' [3.0] Set up background errors (be) for cv_option:', cv_options
  69    write (6,'(A)') ' ----------------------------------------------------------'
  70    write (6,*)
  71
  72    if (cv_options /= 3) then
  73       write(unit=message(1),fmt='(A)') 'Something has gone horribly wrong here'
  74       write(unit=message(2),fmt='(A,I4)') &
  75           'This subroutine is for cv_options = 3, yet cv_options = ', cv_options
  76       call da_error(__FILE__,__LINE__,message(1:2))
  77    endif
  78
  79 !!!!!!!!! cv_options=3
  80    be % v1 % name = 'psi  '           ! Streamfunction
  81    be % v2 % name = 'chi_u'           ! Uncorrelated velocity potential.
  82    be % v3 % name = 't_u'             ! Unbalanced temperature.
  83    be % v4 % name = 'q/qsg'
  84    be % v5 % name = 'psfc'            ! surface pressure
  85    write(6,'(3x,A)')' DA_Setup_Background_Errors: 3DVAR dry control variables are:'
  86    write(6,'(4x,7A)')TRIM(be % v1 % name), ', ', TRIM(be % v2 % name), ', ', &
  87                   TRIM(be % v3 % name), ' and ', TRIM(be % v5 % name)
  88
  89    write(6,'(3x,A,A)')' DA_Setup_Background_Errors: 3DVAR humidity control variable is ',&
  90                      TRIM(be % v4 % name)
  91
  92    write(6,*)
  93
  94    be % mix =  ide -  ids + 1
  95    be % mjy =  jde -  jds + 1
  96
  97    ij = ( ite- its+1)*( jte- jts+1)
  98    ijk = ij * ( kte- kts+1)
  99
 100    be % v1 % mz =  kde- kds+1
 101    be % v2 % mz =  kde- kds+1
 102    be % v3 % mz =  kde- kds+1
 103    be % v4 % mz =  kde- kds+1
 104    be % v5 % mz = 1
 105
 106    be%cv_mz(1) = be%v1%mz
 107    be%cv_mz(2) = be%v2%mz
 108    be%cv_mz(3) = be%v3%mz
 109    be%cv_mz(4) = be%v4%mz
 110    be%cv_mz(5) = be%v5%mz
 111
 112    be % cv % size1  = ijk
 113    be % cv % size2  = ijk
 114    be % cv % size3  = ijk
 115    be % cv % size4  = ijk
 116    be % cv % size5  = ij
 117
 118    be % cv % size = be % cv % size1 + be % cv % size2 + be % cv % size3 + &
 119                     be % cv % size4 + be % cv % size5
 120
 121    cv_size = be % cv % size
 122
 123
 124    call da_get_unit(be_unit)
 125    open (unit=be_unit,file="be.dat",status="old",form="unformatted")
 126
 127    rewind(be_unit)
 128
 129    be_rf_unit    = unit_end + 3
 130
 131    read(be_unit, iostat= ier) nsig,nlath
 132
 133    if (ier /= 0) then
 134       write(unit=message(1), fmt='(A,I4,A,I4)') &
 135            'cv_options:', cv_options,' Reading error in unit=',be_unit
 136       write(unit=message(2), fmt='(A)') 'There is likely a problem with your be.dat file'
 137       call da_error(__FILE__,__LINE__,message(1:2))
 138    endif
 139
 140    write(unit=message(1),fmt='(A,I10)') 'Number of vertical level for stats = ', nsig
 141    write(unit=message(2),fmt='(A,I10)') 'Number of latitude           nlath = ', nlath
 142    call da_message(message(1:2))
 143
 144    kz =  kde- kds+1
 145
 146    if(nsig.ne.kz)then
 147       write(unit=message(1),fmt='(A,I6)') 'Number of vertical level for WRFVar=', kz
 148       call da_message(message(1:1))
 149    end if
 150
 151 ! 1, Allocate the arrays
 152
 153 !   1.1 for BK STATS at WRF eta-levels:
 154
 155   ! Variances for psi(corz), chi_u(cord), t_u(corh), and p-rh(corq)
 156    ALLOCATE ( corz_kz(1:2*nlath+1,1:kz),cord_kz(1:2*nlath+1,1:kz) )
 157    ALLOCATE ( corh_kz(1:2*nlath+1,1:kz),corq_kz(1:2*nlath+1,1:kz) )
 158   ! Scale lengths: horizontal (hwll), vertical (vztdq)
 159    ALLOCATE ( hwll_kz(0:nlath*2+1,1:kz,1:4)                      )
 160    ALLOCATE ( vztdq_kz(1:kz,0:nlath*2+1,1:4)                     )
 161   ! Regression coefficients: t(agv), chi(bv), psfc(wgv)
 162    ALLOCATE ( agv_kz(0:nlath*2+1,1:kz,1:kz)                      )
 163    ALLOCATE ( bv_kz(0:nlath*2+1,1:kz),wgv_kz(0:nlath*2+1,1:kz)   )
 164
 165 !   1.2 for BK STATS inputed from NCEP GFS sigma levels:
 166
 167   ! Variances for psi(corz), chi_u(cord), t_u(corh), and p-rh(corq)
 168    ALLOCATE ( corz_avn(1:2*nlath+1,1:nsig),cord_avn(1:2*nlath+1,1:nsig) )
 169    ALLOCATE ( corh_avn(1:2*nlath+1,1:nsig),corq_avn(1:2*nlath+1,1:nsig) )
 170   ! Variance for Psfc:
 171    ALLOCATE ( corp_avn(1:2*nlath+1),clat_avn(1:2*nlath),sigma_avn(1:nsig) )
 172   ! Scale lengths: horizontal (hwll), vertical (vztdq)
 173    ALLOCATE ( hwll_avn(0:nlath*2+1,1:nsig,1:4),hwllp_avn(0:nlath*2+1) )
 174    ALLOCATE ( vztdq_avn(1:nsig,0:nlath*2+1,1:4)                     )
 175   ! Regression coefficients: t(agv), chi(bv), psfc(wgv)
 176    ALLOCATE ( agv_avn(0:nlath*2+1,1:nsig,1:nsig)                    )
 177    ALLOCATE ( bv_avn(0:nlath*2+1,1:nsig),wgv_avn(0:nlath*2+1,1:nsig) )
 178
 179 !   1.3 for BK STATS at the WRF model grids:
 180
 181    ALLOCATE ( be % corz(its:ite,jts:jte,kts:kte,1:4) )
 182    ALLOCATE ( be % corp(its:ite,jts:jte) )
 183
 184    ALLOCATE ( be % vz(kts:kte,its:ite,jts:jte,1:4) )
 185
 186    ALLOCATE ( be % agvz(its:ite,jts:jte,kts:kte,kts:kte) )
 187    ALLOCATE ( be % bvz(its:ite,jts:jte,kts:kte) )
 188    ALLOCATE ( be % wgvz(its:ite,jts:jte,kts:kte) )
 189 !
 190
 191 ! 2, load the WRF model latitude and map factor:
 192
 193 #ifdef DM_PARALLEL
 194    ij = ( ide- ids+1)*( jde- jds+1)
 195 !  Collect xb component of lat into global buffer.
 196    call da_patch_to_global( grid, grid%xb%lat, global_lat )
 197    call wrf_dm_bcast_real( global_lat, ij )
 198 !  Collect xb component of fac into global buffer.
 199    call da_patch_to_global( grid, grid%xb%map_factor, global_fac )
 200    call wrf_dm_bcast_real( global_fac, ij )
 201 #else
 202    do i = ids, ide
 203       do j = jds, jde
 204          global_lat(i,j) = grid%xb%lat(i,j)
 205          global_fac(i,j) = grid%xb%map_factor(i,j)
 206       enddo
 207    enddo
 208 #endif
 209
 210 ! 3, Read in the NCEP GFS BK STATS:
 211
 212 !   3.1 Latitude and sigma values:
 213
 214    read(be_unit, iostat=ier) clat_avn,(sigma_avn(k),k=1,nsig)
 215    if (ier /= 0) then
 216       write(unit=message(1), fmt='(A,I4,A,I4)') &
 217            'cv_options:', cv_options,' Reading error in unit=',be_unit
 218       write(unit=message(2), fmt='(A)') 'There is likely a problem with your be.dat file'
 219       call da_error(__FILE__,__LINE__,message(1:2))
 220    endif
 221
 222 !   3.2 Variances:
 223
 224    m=2*nlath+1
 225    read (be_unit) &
 226                    ((corz_avn(i,k),i=1,m),k=1,nsig),   &
 227                    ((cord_avn(i,k),i=1,m),k=1,nsig),   &
 228                    ((corh_avn(i,k),i=1,m),k=1,nsig),   &
 229                    ((corq_avn(i,k),i=1,m),k=1,nsig),corp_avn
 230
 231 !   3.3 Scale lengths
 232
 233    ! horizontal
 234    read(be_unit) (((hwll_avn(i,k,m),i=0,nlath*2+1),k=1,nsig),m=1,4),   &
 235                      hwllp_avn
 236    ! vertical:
 237    read(be_unit) (((vztdq_avn(k,i,m),k=1,nsig),i=0,nlath*2+1),m=1,4)
 238
 239 !   3.4 Regression coefficients:
 240
 241    read(be_unit) (((agv_avn(i,k,m),i=0,nlath*2+1),k=1,nsig),m=1,nsig), &
 242                          ((bv_avn(i,k),i=0,nlath*2+1),k=1,nsig), &
 243                          ((wgv_avn(i,k),i=0,nlath*2+1),k=1,nsig)
 244
 245 ! 4, incorporate the tuning factors:
 246
 247 !   4.1 Horizontal scales:
 248
 249    as(1)=as1(2)
 250    as(2)=as2(2)
 251    as(3)=as3(2)
 252    as(4)=as4(2)
 253    do m=1,4
 254       do k=1,nsig
 255          do i=0,nlath*2+1
 256             hwll_avn(i,k,m)=hwll_avn(i,k,m)*as(m)
 257          enddo
 258       enddo
 259    enddo
 260    do i=0,nlath*2+1
 261       hwllp_avn(i)=hwllp_avn(i)*as5(2)
 262    enddo
 263
 264 !   4.2 Vertical scales:
 265    as(1)=as1(3)
 266    as(2)=as2(3)
 267    as(3)=as3(3)
 268    as(4)=as4(3)
 269    do m=1,4
 270       do i=0,nlath*2+1
 271          do k=1,nsig
 272             vztdq_avn(k,i,m)=vztdq_avn(k,i,m)*as(m)
 273          enddo
 274        enddo
 275    enddo
 276
 277 ! 5, determine the level ku, which locates above and nearest sigma=0.15:
 278
 279    ku= kde+1
 280    k_loop: do k= kds+1, kde
 281       if (grid%xb%sigmah(k-1)>0.15 .and. grid%xb%sigmah(k)<=0.15) then
 282          ku=k
 283          exit k_loop
 284       endif
 285    end do k_loop
 286    write(mesg,'(a,i5,a)') "==> level ku =", ku ," locates above and nearest sigma=0.15"
 287    call wrf_debug ( 1 , mesg )
 288
 289 ! 6, Vertical interpolation of BK STATS:
 290 !    to convert the NCEP GFS sigma levels to WRF model eta levels
 291
 292    call da_chgvres(nlath,nsig,kz,grid%xb%sigmah,sigma_avn,&
 293                    corz_avn,cord_avn,corh_avn,corq_avn,hwll_avn,vztdq_avn,agv_avn,bv_avn,wgv_avn,&
 294                    corz_kz, cord_kz, corh_kz, corq_kz, hwll_kz, vztdq_kz, agv_kz, bv_kz, wgv_kz)
 295
 296 ! 7, Horizontal interpolation
 297
 298 !   7.1 Calculate the interpolation coefficients:
 299
 300    !$OMP PARALLEL DO &
 301    !$OMP PRIVATE ( j, i, m, m1 )
 302    do j= jds, jde
 303       do i= ids,  ide
 304
 305          if (global_lat(i,j).ge.clat_avn(2*nlath)) then
 306             ! 7.1.1 Model lat >= max AVN lat:
 307             mlat(i,j)=nlath*2-1
 308             coef1(i,j)=0.
 309             coef2(i,j)=1.
 310          else
 311             ! 7.1.2 Model lat < max AVN lat:
 312             do m=1,2*nlath-1
 313                m1=m+1
 314                if ((global_lat(i,j).ge.clat_avn(m)).and.  &
 315                   (global_lat(i,j).lt.clat_avn(m1))) then
 316                   mlat(i,j)=m
 317                   exit
 318                end if
 319             end do
 320
 321             coef2(i,j)=(global_lat(i,j)-clat_avn(m))/(clat_avn(m1)-clat_avn(m))
 322             coef1(i,j)=1.-coef2(i,j)
 323          endif
 324
 325       end do
 326    end do
 327    !$OMP END PARALLEL DO
 328
 329 !   7.2 interpolation of the covariance
 330
 331    ! Psfc:
 332    !$OMP PARALLEL DO &
 333    !$OMP PRIVATE ( ij, j, i, m, m1 )
 334    do ij = 1, grid%num_tiles
 335
 336    do j=grid%j_start(ij), grid%j_end(ij)
 337       do i=its,ite
 338          m=mlat(i,j)
 339          m1=m+1
 340          be%corp(i,j)=corp_avn(m)*coef1(i,j)+corp_avn(m1)*coef2(i,j)
 341       enddo
 342    enddo
 343
 344    enddo
 345    !$OMP END PARALLEL DO
 346
 347    ! psi, chi, t, and rh:
 348    !$OMP PARALLEL DO &
 349    !$OMP PRIVATE ( ij, j, i, k, m, m1 )
 350    do ij = 1, grid%num_tiles
 351
 352    do k=kts,kte
 353       do j=grid%j_start(ij), grid%j_end(ij)
 354          do i=its,ite
 355             m=mlat(i,j)
 356             m1=m+1
 357             be%corz(i,j,k,1)=corz_kz(m,k)*coef1(i,j)+corz_kz(m1,k)*coef2(i,j)
 358             be%corz(i,j,k,2)=cord_kz(m,k)*coef1(i,j)+cord_kz(m1,k)*coef2(i,j)
 359             be%corz(i,j,k,3)=corh_kz(m,k)*coef1(i,j)+corh_kz(m1,k)*coef2(i,j)
 360             be%corz(i,j,k,4)=corq_kz(m,k)*coef1(i,j)+corq_kz(m1,k)*coef2(i,j)
 361          end do
 362       end do
 363    end do
 364
 365    end do
 366    !$OMP END PARALLEL DO
 367
 368 !   7.3 interpolation of the horizontal scale lengths
 369
 370    ic = 1 + ( ide- ids)/2
 371    jc = 1 + ( jde- jds)/2
 372
 373    ! Psfc:
 374    do j= jds,  jde
 375       do i= ids,  ide
 376          m=mlat(i,j)
 377          m1=m+1
 378          hwllp(i,j)=hwllp_avn(m)*coef1(i,j)+hwllp_avn(m1)*coef2(i,j)
 379       end do
 380    end do
 381    write(mesg,'(a,2i5,a,f20.10)') 'Horizontal scale length (m) for Psfc at (',ic/2,jc/2,') hwllp=',hwllp(ic/2,jc/2)
 382    call wrf_debug ( 1 , mesg )
 383
 384    ! psi, chi, t, and rh:
 385    do n4=1,4
 386       do k= kds,  kde
 387          do j= jds,  jde
 388             do i= ids,  ide
 389                m=mlat(i,j)
 390                m1=m+1
 391                hwll(i,j,k,n4)=hwll_kz(m,k,n4)*coef1(i,j)+hwll_kz(m1,k,n4)*coef2(i,j)
 392             end do
 393          end do
 394       end do
 395    end do
 396
 397 !   7.4 interpolation of the vertical scale lengths
 398
 399    do n4=1,4
 400       do j=jts,jte
 401          do i=its,ite
 402             m=mlat(i,j)
 403             m1=m+1
 404             do k=kts,kte
 405                be%vz(k,i,j,n4)=vztdq_kz(k,m,n4)*coef1(i,j)+vztdq_kz(k,m1,n4)*coef2(i,j)
 406             end do
 407          end do
 408       end do
 409    end do
 410    write(mesg,'(a,3i5,a,4f20.16)') 'Vertical scale length for Psi, Chi, t, and rh at (k,i,j)=', &
 411                                   10,its+2,jts+5,'  vz=',(be%vz(10,its+2,jts+5,n4),n4=1,4)
 412    call wrf_debug ( 1 , mesg )
 413
 414 !   7.5 interpolation of the regression coefficients
 415
 416    ! Temperature:
 417    do k=kts,kte
 418       do n=kts,kte
 419          do j=jts,jte
 420             do i=its,ite
 421                m=mlat(i,j)
 422                m1=m+1
 423                be%agvz(i,j,n,k)=agv_kz(m,n,k)*coef1(i,j)+agv_kz(m1,n,k)*coef2(i,j)
 424             end do
 425          end do
 426       end do
 427    end do
 428
 429    ! potential velocity:
 430    do k=kts,kte
 431       do j=jts,jte
 432          do i=its,ite
 433             m=mlat(i,j)
 434             m1=m+1
 435             be%bvz(i,j,k)=bv_kz(m,k)*coef1(i,j)+bv_kz(m1,k)*coef2(i,j)
 436          end do
 437       end do
 438    end do
 439
 440    ! Surface pressure:
 441    do k=kts,kte
 442       do j=jts,jte
 443          do i=its,ite
 444             m=mlat(i,j)
 445             m1=m+1
 446             be%wgvz(i,j,k)=wgv_kz(m,k)*coef1(i,j)+wgv_kz(m1,k)*coef2(i,j)
 447          end do
 448       end do
 449    end do
 450
 451 !   7.6 Deallocate the arrays:
 452
 453    ! For NCEP GFS BK STATS:
 454    DEALLOCATE ( corz_avn,cord_avn )
 455    DEALLOCATE ( corh_avn,corq_avn )
 456    DEALLOCATE ( corp_avn,clat_avn,sigma_avn )
 457    DEALLOCATE ( hwll_avn,hwllp_avn )
 458    DEALLOCATE ( vztdq_avn )
 459    DEALLOCATE ( agv_avn )
 460    DEALLOCATE ( bv_avn,wgv_avn )
 461
 462    ! For WRF model levels:
 463    DEALLOCATE ( corz_kz,cord_kz )
 464    DEALLOCATE ( corh_kz,corq_kz )
 465    DEALLOCATE ( hwll_kz )
 466    DEALLOCATE ( vztdq_kz )
 467    DEALLOCATE ( agv_kz )
 468    DEALLOCATE ( bv_kz,wgv_kz )
 469
 470 !
 471 ! 8, Create the parameters for recursive filter
 472
 473 !  call da_prerf(grid%xb,be)
 474 !
 475 !   8.1 Set the constants and generate be%be, samp, be%table:
 476
 477    ! Constant-1: ndeg:
 478    ndeg=4
 479
 480    be%ndeg=ndeg
 481
 482    ALLOCATE ( turn (1:ndeg,1:ndeg) )
 483    ALLOCATE ( be % be (1:ndeg) )
 484    ALLOCATE ( be % rate (1:ndeg) )
 485
 486    CALL RFDPAR1(be%BE,be%RATE,ndeg)
 487    CALL RFDPAR2(be%BE,be%RATE,TURn,SAMP,ndeg)
 488
 489     ! Constant-2: nta:
 490    nta=5600
 491
 492    be%nta=nta
 493    allocate (dsh(1:nta)        )
 494    ALLOCATE ( be % table (1:nta,1:ndeg) )
 495
 496     ! Constant-3: be%swidth:
 497    be%swidth=10.
 498
 499    tin=be%swidth/dble(nta)
 500    do i=1,nta
 501       dsh(i)=dble(i-1)*tin
 502    enddo
 503
 504    call  RFDPARV(DSH,be%RATE,be%table,nta,ndeg )
 505
 506 !   8.2 Deallocate the working arrays:
 507
 508    deallocate (dsh )
 509    deallocate (turn )
 510    be_print_unit = unit_end + 4
 511
 512 !  8.2.1 Save the CV3 BE at model grids before tuned:
 513
 514 if ( max_ext_its > 1 ) then
 515 #ifdef DM_PARALLEL
 516    write (*,'(/a,i3)') 'mpp code ==> Write CV3 BE to fort.',be_rf_unit
 517
 518    if (rootproc) then
 519    write (be_rf_unit) ids, ide, jds, jde, kds, kde, ku, samp
 520    write (*,'(a,7i5,f15.5)') "ids, ide, jds, jde, kds, kde, ku, samp:", &
 521                         ids, ide, jds, jde, kds, kde, ku, samp
 522    endif
 523
 524    ij = ( ide - ids + 1)*( jde - jds + 1)
 525    ijk = ij * ( kde - kds + 1)
 526 !  Collect the variance and write out:
 527    do ii = 1, 4
 528 ! tile --> patch:
 529      corz_3d(its:ite,jts:jte,kts:kte) = be%corz(its:ite,jts:jte,kts:kte,ii)
 530 ! patch -->global:
 531      call da_patch_to_global(grid, corz_3d, global_3d)
 532 ! broadcast:
 533      call wrf_dm_bcast_real(global_3d, ijk)
 534
 535      if (ii == 1) vname = 'PSI   '
 536      if (ii == 2) vname = 'CHI_u '
 537      if (ii == 3) vname = 'TMP_u '
 538      if (ii == 4) vname = 'PSD_RH'
 539
 540      xxsum(1) = sum (be%corz(:,:,:,ii)*be%corz(:,:,:,ii))
 541      call da_proc_sum_real(xxsum)
 542      if (rootproc) then
 543        write (be_rf_unit) 'VARIANCE:', vname, ii, global_3d
 544 !       write (*,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xxsum(1)
 545      endif
 546    enddo
 547 ! tile --> patch:
 548    corp_2d(its:ite,jts:jte) = be%corp(its:ite,jts:jte)
 549 ! patch --> global:
 550    call da_patch_to_global(grid, corp_2d, global_2d)
 551 ! broadcast:
 552    call wrf_dm_bcast_real(global_2d, ij)
 553    xxsum(1) = sum (be%corp*be%corp)
 554    call da_proc_sum_real(xxsum)
 555    if (rootproc) then
 556      write (be_rf_unit) 'VARIANCE:', 'PSFC_u',  1,  global_2d
 557 !     write (*,'(9x,a,2x,"sum^2=",e20.12)') 'PSFC_u', xxsum(1)
 558    endif
 559 !
 560    do ii = 1, 4
 561 ! tile --> patch:
 562      do i = its, ite
 563      do j = jts, jte
 564      do k = kts,kte
 565         vz_3d(i,j,k) = be%vz(k,i,j,ii)
 566      enddo
 567      enddo
 568      enddo
 569 ! patch -->global:
 570      call da_patch_to_global(grid, vz_3d, global_3d)
 571 ! broadcast:
 572      call wrf_dm_bcast_real(global_3d, ijk)
 573
 574      if (ii == 1) vname = 'PSI   '
 575      if (ii == 2) vname = 'CHI_u '
 576      if (ii == 3) vname = 'TMP_u '
 577      if (ii == 4) vname = 'PSD_RH'
 578
 579      xxsum(1) = sum (be%vz(:,:,:,ii)*be%vz(:,:,:,ii))
 580      call da_proc_sum_real(xxsum)
 581      if (rootproc) then
 582        write (be_rf_unit) 'VZ-SCALE:', vname, ii, global_3d
 583 !       write (*,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xxsum(1)
 584      endif
 585    enddo
 586 !
 587 !  Horizontal scales:
 588     if (rootproc) then
 589       write (be_rf_unit) hwll
 590       xsum = sum (hwll*hwll)
 591 !      write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PSI, CHI_u, T_u, RH_u SCALES:', xsum
 592       write (be_rf_unit) hwllp
 593       xsum = sum (hwllp*hwllp)
 594 !      write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xsum
 595     endif
 596
 597 #else
 598    write (*,'(/a,i3)') 'Serial code ==> Write CV3 BE to fort.',be_rf_unit
 599
 600    write (be_rf_unit) ids, ide, jds, jde, kds, kde, ku, samp
 601    write (*,'(a,7i5,f15.5)') "ids, ide, jds, jde, kds, kde, ku, samp:", &
 602                         ids, ide, jds, jde, kds, kde, ku, samp
 603
 604 ! Save the variances before the normalized:
 605 !   write (*,'(3x,a)')  'VARIANCE:'
 606    do ii = 1, 4
 607      if (ii == 1) vname = 'PSI   '
 608      if (ii == 2) vname = 'CHI_u '
 609      if (ii == 3) vname = 'TMP_u '
 610      if (ii == 4) vname = 'PSD_RH'
 611      global_3d = be%corz(:,:,:,ii)
 612      write (be_rf_unit) 'VARIANCE:', vname, ii, global_3d
 613      xsum = sum (be%corz(:,:,:,ii)*be%corz(:,:,:,ii))
 614 !     write (*,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xsum
 615    enddo
 616    global_2d = be%corp
 617    write (be_rf_unit) 'VARIANCE:', 'PSFC_u',  1,  global_2d
 618    xsum = sum (be%corp*be%corp)
 619 !   write (*,'(9x,a,2x,"sum^2=",e20.12)') 'PSFC_u', xsum
 620 !
 621 !   write (*,'(3x,a)')  'VERTICAL SCALES:'
 622    do ii = 1, 4
 623      if (ii == 1) vname = 'PSI   '
 624      if (ii == 2) vname = 'CHI_u '
 625      if (ii == 3) vname = 'TMP_u '
 626      if (ii == 4) vname = 'PSD_RH'
 627      do i = ids, ide
 628      do j = jds, jde
 629      do k = kds, kde
 630         global_3d(i,j,k) = be%vz(k,i,j,ii)
 631      enddo
 632      enddo
 633      enddo
 634      write (be_rf_unit) 'VZ-SCALE:', vname, ii, global_3d
 635      xsum = sum (be%vz(:,:,:,ii)*be%vz(:,:,:,ii))
 636 !     write (*,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xsum
 637    enddo
 638 !
 639 !  Horizontal scales:
 640 !   write (*,'(3x,a)')  'HORIZONTAL SCALES:'
 641    write (be_rf_unit) hwll
 642    xsum = sum (hwll*hwll)
 643 !   write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PSI, CHI_u, T_u, RH_u SCALES:', xsum
 644
 645    write (be_rf_unit) hwllp
 646    xsum = sum (hwllp*hwllp)
 647 !   write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xsum
 648
 649 #endif
 650
 651 !  8.2.2 Write out the CV3 BE information to "be_print_unit":
 652
 653    if (PRINT_DETAIL_BE) then
 654
 655    write (be_print_unit,'(a,7i5,f15.5)') "ids, ide, jds, jde, kds, kde, ku, samp:", &
 656                                           ids, ide, jds, jde, kds, kde, ku, samp
 657
 658    write (be_print_unit,'(3x,a)')  'VARIANCE:'
 659    do ii = 1, 4
 660      if (ii == 1) vname = 'PSI   '
 661      if (ii == 2) vname = 'CHI_u '
 662      if (ii == 3) vname = 'TMP_u '
 663      if (ii == 4) vname = 'PSD_RH'
 664      xxsum(1) = sum (be%corz(its:ite,jts:jte,kts:kte,ii)*be%corz(its:ite,jts:jte,kts:kte,ii))
 665      call da_proc_sum_real(xxsum)
 666      if (rootproc) &
 667      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xxsum(1)
 668    enddo
 669
 670 !  Pscf Variance before the normalization.
 671    xxsum(1) = sum (be%corp(its:ite,jts:jte)*be%corp(its:ite,jts:jte))
 672    call da_proc_sum_real(xxsum)
 673    if (rootproc) &
 674    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'PSFC_u', xxsum(1)
 675 !
 676 !  Collect the vertical scales and write out:
 677    write (be_print_unit,'(3x,a)')  'VERTICAL SCALES:'
 678    do ii = 1, 4
 679      if (ii == 1) vname = 'PSI   '
 680      if (ii == 2) vname = 'CHI_u '
 681      if (ii == 3) vname = 'TMP_u '
 682      if (ii == 4) vname = 'PSD_RH'
 683      xxsum(1) = sum (be%vz(:,:,:,ii)*be%vz(:,:,:,ii))
 684      call da_proc_sum_real(xxsum)
 685      if (rootproc) &
 686      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xxsum(1)
 687    enddo
 688
 689 !  Horizontal scales:
 690    if (rootproc) &
 691    write (be_print_unit,'(3x,a)')  'READ IN THE HORIZONTAL SCALES:'
 692    if (rootproc) then
 693    write (be_print_unit,'(3x,a)')  'HORIZONTAL SCALES:'
 694    xxsum(1) = sum (hwll*hwll)
 695    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PSI, CHI_u, T_u, RH_u SCALES:', xxsum(1)
 696    endif
 697
 698    if (rootproc) then
 699    xxsum(1) = sum (hwllp*hwllp)
 700    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xxsum(1)
 701    endif
 702
 703 !  Collect the regression coefficients:
 704    write (be_print_unit,'(3x,a)')  'REGRESSION COEFF. T, CHI, PSFC:'
 705    do ii = kts, kte
 706      xxsum(1) = sum (be%agvz(:,:,:,ii)*be%agvz(:,:,:,ii))
 707      call da_proc_sum_real(xxsum)
 708      if (rootproc) &
 709      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, 'TMP-PSI', xxsum(1)
 710    enddo
 711
 712 !  Rg. Coeff. has already stored in "be" data structure, not need to be save.
 713    xxsum(1) = sum (be%bvz*be%bvz)
 714    call da_proc_sum_real(xxsum)
 715    if (rootproc) &
 716    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'CHI-PSI', xxsum(1)
 717
 718 !  Rg. Coeff. has already stored in "be" data structure, not need to be save.
 719    xxsum(1) = sum (be%wgvz*be%wgvz)
 720    call da_proc_sum_real(xxsum)
 721    if (rootproc) &
 722    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'PSF-PSI', xxsum(1)
 723
 724    endif
 725 end if   ! max_ext_its > 1
 726
 727 ! 9, Incorporate the tuning factors for covariance
 728
 729 ! sli in scale  unit (map_factor come with ds )
 730 !           variance* amp for 3d/2d RF
 731
 732    as(1)=sqrt(as1(1))
 733    as(2)=sqrt(as2(1))
 734    as(3)=sqrt(as3(1))
 735    as(4)=sqrt(as4(1))
 736    as(5)=sqrt(as5(1))
 737
 738 !   9.1 Scale the horizontal scale in unit of grid-point:
 739
 740    s2u= 1./grid%xb%ds
 741    hwll=hwll*s2u
 742    hwllp=hwllp*s2u
 743
 744 !   9.2 Re-scale the covariance for psi, chi, t, and rh:
 745
 746    !$OMP PARALLEL DO &
 747    !$OMP PRIVATE (ij, n, j, i, vv, k)
 748    do ij = 1, grid%num_tiles
 749    do n=1,4
 750       do j=grid%j_start(ij), grid%j_end(ij)
 751          do i=its,ite
 752
 753             vv=0.
 754             do k=kts,kte
 755                vv(k,k)=1.
 756             enddo
 757
 758      ! Recursive filter routie applied in vertical with
 759      ! the vertical scale length be%vz:
 760             call da_rfz0(vv,kz,kz,be%ndeg,&
 761                          be%vz(kts:kte,i,j,n),be%be,be%table,be%nta,be%swidth)
 762
 763      ! Re-scale the covariance for psi, chi, t, and rh:
 764             do k=kts,kte
 765                be % corz(i,j,k,n)=be % corz(i,j,k,n)*as(n) &
 766                                   *samp/hwll(i,j,k,n)/vv(k,k)/global_fac(i,j)
 767             enddo
 768
 769          enddo
 770       enddo
 771    enddo
 772    enddo
 773    !$OMP END PARALLEL DO
 774
 775 !   9.3 Re-scale the covariance for Psfc:
 776
 777     be % corp(its:ite,jts:jte)=be % corp(its:ite,jts:jte)*as(5) &
 778          *samp/hwllp(its:ite,jts:jte)/global_fac(its:ite,jts:jte)
 779
 780     write(mesg,*) 'Re-scaled covariance for Psfc: sum(be%corp*be%corp)=', &
 781                   sum(be%corp*be%corp)
 782     call wrf_debug ( 1 , mesg )
 783
 784 !
 785 ! 10, Assign the inverse of scale length fields for recursive filter:
 786 !
 787 !   10.1 allocate the arrays for be scales: y-direction and x-direction:
 788
 789     ALLOCATE ( be % sljpy ( grid%xp%ipsy: grid%xp%ipey, grid%xp%jpsy: grid%xp%jpey) )
 790     ALLOCATE ( be % sljy ( grid%xp%ipsy: grid%xp%ipey, grid%xp%jpsy: grid%xp%jpey, grid%xp%kpsy: grid%xp%kpey,1:4) )
 791
 792     ALLOCATE ( be % slipx ( grid%xp%ipsx: grid%xp%ipex, grid%xp%jpsx: grid%xp%jpex) )
 793     ALLOCATE ( be % slix ( grid%xp%ipsx: grid%xp%ipex, grid%xp%jpsx: grid%xp%jpex, grid%xp%kpsx: grid%xp%kpex,1:4) )
 794
 795 !   10.2 Y-direction:
 796
 797     ! 3-D fields: psi, chi, t, and rh:
 798     do n=1,4
 799        do k= grid%xp%kpsy, grid%xp%kpey
 800           do j= grid%xp%jpsy, grid%xp%jpey
 801              do i= grid%xp%ipsy, grid%xp%ipey
 802                 be%sljy(i,j,k,n)=1./global_fac(i,j)/hwll(i,j,k,n)
 803              enddo
 804           enddo
 805        enddo
 806     enddo
 807
 808     ! Above level ku,the sljy fields are set to a constant
 809     ! for psi and chi, i.e. homogenous:
 810     do n=1,2
 811        do k=max(ku, grid%xp%kpsy), grid%xp%kpey
 812           slim=1./global_fac(ic,jc)/hwll(ic,jc,k,n)
 813           do j= grid%xp%jpsy, grid%xp%jpey
 814              do i= grid%xp%ipsy, grid%xp%ipey
 815                 be%sljy(i,j,k,n)=slim
 816              enddo
 817           enddo
 818        enddo
 819     enddo
 820
 821     ! 2-D field: Psfc:
 822     do j= grid%xp%jpsy, grid%xp%jpey
 823        do i= grid%xp%ipsy, grid%xp%ipey
 824           be%sljpy(i,j)=1./global_fac(i,j)/hwllp(i,j)
 825        enddo
 826     enddo
 827
 828 !   10.3 X-direction:
 829
 830    ! 3-D fields: psi, chi, t, and rh:
 831    do n=1,4
 832       do k= grid%xp%kpsx, grid%xp%kpex
 833          do j= grid%xp%jpsx, grid%xp%jpex
 834             do i= grid%xp%ipsx, grid%xp%ipex
 835                be%slix(i,j,k,n)=1./global_fac(i,j)/hwll(i,j,k,n)
 836             enddo
 837          enddo
 838       enddo
 839    enddo
 840
 841    ! Above level ku,the sljy fields are set to a constant
 842    ! for psi and chi, i.e. homogenous:
 843    do n=1,2
 844       do k=max(ku, grid%xp%kpsx), grid%xp%kpex
 845          slim=1./global_fac(ic,jc)/hwll(ic,jc,k,n)
 846          do j= grid%xp%jpsx, grid%xp%jpex
 847             do i= grid%xp%ipsx, grid%xp%ipex
 848                be%slix(i,j,k,n)=slim
 849             enddo
 850          enddo
 851       enddo
 852    enddo
 853
 854    ! 2-D field: Psfc:
 855    do j= grid%xp%jpsx, grid%xp%jpex
 856       do i= grid%xp%ipsx, grid%xp%ipex
 857          be%slipx(i,j)=1./global_fac(i,j)/hwllp(i,j)
 858       enddo
 859    enddo
 860
 861    close(be_unit)
 862    call da_free_unit(be_unit)
 863
 864 end subroutine da_setup_be_ncep_gfs
 865