var/da/da_setup_structures/da_scale_background_errors_cv3.inc

   1 subroutine da_scale_background_errors_cv3 ( grid, be, it )
   2 !
   3    type (domain), intent(in)     :: grid
   4    TYPE (be_type), INTENT(INOUT) :: be     ! Back. errors structure
   5    INTEGER,        INTENT(IN)    :: it     ! outer-loop index
   6 !
   7    REAL, DIMENSION( ids: ide,  jds: jde,  kds: kde, 1:4) :: hwll
   8    REAL, DIMENSION( ids: ide,  jds: jde)                 :: hwllp, global_fac
   9    REAL, DIMENSION( kts: kte,  kts: kte)                 :: vv
  10    integer                     :: n, i, j, k, ic, jc, ii, ij, ijk
  11    integer, dimension(1)       :: iis, iie, jjs, jje, kks, kke, ivar
  12    real, dimension(1)          :: xsum
  13    real, dimension(ids: ide,  jds: jde)                   :: global_2d
  14    real, dimension(ids: ide,  jds: jde,  kds: kde)        :: global_3d
  15    character(len=6)            :: vname
  16    character(len=9)            :: chr
  17    character(len=8)            :: i_char
  18
  19    INTEGER :: nta, ndeg, ku(1), kz
  20    real    :: samp(1),s2u,tin,as(5),as0(5),slim
  21    character(len=256) :: mesg
  22
  23    integer                     :: ier, be_rf_unit, be_print_unit
  24 !
  25 ! 1, get global_fac from grid:
  26
  27 #ifdef DM_PARALLEL
  28    ij = ( ide- ids+1)*( jde- jds+1)
  29    ijk = ij * ( kde - kds + 1)
  30 !  Collect xb component of fac into global buffer.
  31    call da_patch_to_global( grid, grid%xb%map_factor, global_fac )
  32    call wrf_dm_bcast_real( global_fac, ij )
  33 #else
  34    do i = ids, ide
  35       do j = jds, jde
  36          global_fac(i,j) = grid%xb%map_factor(i,j)
  37       enddo
  38    enddo
  39 #endif
  40
  41 ! 2, Read in the variance and horizontal scales
  42 !
  43 ! 2.1 assign the units for reading and printing:
  44
  45 ! Rewind the unit:
  46     be_rf_unit    = unit_end + 3
  47     be_print_unit = unit_end + 5
  48     if (rootproc) rewind (be_rf_unit)
  49
  50     if (rootproc .and. print_detail_be) &
  51     write(be_print_unit,'(//10x,a,i2,a)') '====== SCALE CV3 BE for IT =', it, ' ======'
  52 !
  53 ! 2.2 Read the dimensions
  54 !
  55     if (rootproc) &
  56        read(be_rf_unit) iis, iie, jjs, jje, kks, kke, ku, samp
  57     call wrf_dm_bcast_integer( iis , 1 )
  58     call wrf_dm_bcast_integer( iie , 1 )
  59     call wrf_dm_bcast_integer( jjs , 1 )
  60     call wrf_dm_bcast_integer( jje , 1 )
  61     call wrf_dm_bcast_integer( kks , 1 )
  62     call wrf_dm_bcast_integer( kke , 1 )
  63     call wrf_dm_bcast_integer( ku , 1 )
  64     call wrf_dm_bcast_real( samp , 1 )
  65     if (rootproc .and. print_detail_be) &
  66     write (be_print_unit,'(/a,7i5,f15.5)') &
  67               "Read in: ids, ide, jds, jde, kds, kde, ku, samp:", &
  68                         iis, iie, jjs, jje, kks, kke, ku, samp
  69 !
  70 ! 2.3 Read in the variances before normalization and write to "be_print_unit":
  71
  72 !    write (*,'(3x,a)')  'READ IN VARIANCE BEFORE NORMALIZATION:'
  73 ! Print out the sum(x*X) for verification:
  74 !
  75    if (rootproc .and. print_detail_be) &
  76    write (be_print_unit,'(3x,a)')  'VARIANCE:'
  77    do ii = 1, 4
  78      if (rootproc) &
  79         read (be_rf_unit) chr, vname, ivar, global_3d
  80      call wrf_dm_bcast_string( chr, 9 )
  81      call wrf_dm_bcast_string( vname, 6 )
  82      call wrf_dm_bcast_integer( ivar, 1 )
  83      call wrf_dm_bcast_real( global_3d, ijk )
  84      be%corz(its:ite,jts:jte,kts:kte,ivar(1)) = global_3d(its:ite,jts:jte,kts:kte)
  85      xsum(1) = sum (be%corz(its:ite,jts:jte,kts:kte,ivar(1))*be%corz(its:ite,jts:jte,kts:kte,ivar(1)))
  86      call da_proc_sum_real(xsum)
  87      if (rootproc .and. print_detail_be) &
  88      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xsum(1)
  89 !     write (*,'(3x,a9,2x,a6,2x,i3)') chr, vname, ivar
  90    enddo
  91
  92 !  Psfc Variance before the normalization.
  93    if (rootproc) &
  94       read (be_rf_unit) chr, vname, ivar, global_2d
  95    call wrf_dm_bcast_string( chr, 9 )
  96    call wrf_dm_bcast_string( vname, 6 )
  97    call wrf_dm_bcast_integer( ivar, 1 )
  98    call wrf_dm_bcast_real( global_2d, ij )
  99    be%corp(its:ite,jts:jte) = global_2d(its:ite,jts:jte)
 100    xsum(1) = sum (be%corp(its:ite,jts:jte)*be%corp(its:ite,jts:jte))
 101    call da_proc_sum_real(xsum)
 102    if (rootproc .and. print_detail_be) &
 103    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'PSFC_u', xsum(1)
 104 !   write (*,'(3x,a9,2x,a6,2x,i3)') chr, vname, ivar
 105 !
 106 ! 2.4 Read in the vertical scales to "be_print_unit":
 107 !
 108 !   write (*,'(3x,a)')  'READ IN VERTICAL SCALE BEFORE TUNING:'
 109 !
 110    if (rootproc .and. print_detail_be) &
 111    write (be_print_unit,'(3x,a)')  'VERTICAL SCALE:'
 112    do ii = 1, 4
 113      if (rootproc) &
 114          read (be_rf_unit) chr, vname, ivar, global_3d
 115      call wrf_dm_bcast_string( chr, 9 )
 116      call wrf_dm_bcast_string( vname, 6 )
 117      call wrf_dm_bcast_integer( ivar, 1 )
 118      call wrf_dm_bcast_real( global_3d, ijk )
 119      do i = its, ite
 120      do j = jts, jte
 121      do k = kts, kte
 122         be%vz(k,i,j,ivar(1)) = global_3d(i,j,k)
 123      enddo
 124      enddo
 125      enddo
 126      xsum(1) = sum (be%vz(kts:kte,its:ite,jts:jte,ivar(1))*be%vz(kts:kte,its:ite,jts:jte,ivar(1)))
 127      call da_proc_sum_real(xsum)
 128      if (rootproc .and. print_detail_be) &
 129      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, vname, xsum(1)
 130 !     write (*,'(3x,a9,2x,a6,2x,i3)') chr, vname, ivar
 131    enddo
 132 !
 133 ! 2.5 Read in the Horizontal scales, and write out vertical scales to "be_print_unit":
 134
 135    if (rootproc .and. print_detail_be) &
 136    write (be_print_unit,'(3x,a)')  'READ IN THE HORIZONTAL SCALES:'
 137    if (rootproc) &
 138       read (be_rf_unit) hwll
 139    call wrf_dm_bcast_real( hwll, 4*ijk )
 140    if (rootproc .and. print_detail_be) then
 141    write (be_print_unit,'(3x,a)')  'HORIZONTAL SCALES:'
 142    xsum(1) = sum (hwll*hwll)
 143    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PSI, CHI_u, T_u, RH_u SCALES:', xsum(1)
 144    endif
 145
 146    if (rootproc) &
 147       read (be_rf_unit) hwllp
 148    call wrf_dm_bcast_real( hwllp, ij )
 149    if (rootproc .and. print_detail_be) then
 150    xsum(1) = sum (hwllp*hwllp)
 151    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xsum(1)
 152 !   write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xsum(1)
 153    endif
 154
 155 ! 2.6 Write out the regression coefficients to "be_print_unit":
 156
 157    if (print_detail_be) then
 158
 159    if (rootproc) &
 160    write (be_print_unit,'(3x,a)')  'REGRESSION COEFF. T, CHI, PSFC:'
 161    do ii = kts, kte
 162      xsum(1) = sum (be%agvz(:,:,:,ii)*be%agvz(:,:,:,ii))
 163      call da_proc_sum_real(xsum)
 164      if (rootproc) &
 165      write (be_print_unit,'(5x,i3,1x,a,2x,"sum^2=",e20.12)') ii, 'TMP-PSI', xsum(1)
 166    enddo
 167
 168 !  Rg. Coeff. already stored in "be" data structure, no need to save.
 169    xsum(1) = sum (be%bvz*be%bvz)
 170    call da_proc_sum_real(xsum)
 171    if (rootproc) &
 172    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'CHI-PSI', xsum(1)
 173
 174 !  Rg. Coeff. already stored in "be" data structure, no need to save.
 175    xsum(1) = sum (be%wgvz*be%wgvz)
 176    call da_proc_sum_real(xsum)
 177    if (rootproc) &
 178    write (be_print_unit,'(9x,a,2x,"sum^2=",e20.12)') 'PSF-PSI', xsum(1)
 179 !
 180    endif
 181
 182 ! 3.0 Re-scale the scale-lengths
 183 !
 184 ! 3.1 horizontal scales
 185
 186    as0(1)=as1(2)
 187    as0(2)=as2(2)
 188    as0(3)=as3(2)
 189    as0(4)=as4(2)
 190
 191    as(1)=as1(2+(it-1)*3)
 192    as(2)=as2(2+(it-1)*3)
 193    as(3)=as3(2+(it-1)*3)
 194    as(4)=as4(2+(it-1)*3)
 195
 196    do n=1, 4
 197       hwll(:,:,:,n) = hwll(:,:,:,n) * as(n) / as0(n)
 198    enddo
 199
 200    hwllp = hwllp * as5(2+(it-1)*3) / as5(2)
 201
 202 ! 3.2 vertical scales
 203
 204    as0(1)=as1(3)
 205    as0(2)=as2(3)
 206    as0(3)=as3(3)
 207    as0(4)=as4(3)
 208
 209    as(1)=as1(3+(it-1)*3)
 210    as(2)=as2(3+(it-1)*3)
 211    as(3)=as3(3+(it-1)*3)
 212    as(4)=as4(3+(it-1)*3)
 213
 214    do n=1, 4
 215       be%vz(kts:kte,its:ite,jts:jte,n) = be%vz(kts:kte,its:ite,jts:jte,n) * as(n) / as0(n)
 216    enddo
 217
 218 ! 4, Normalize the variances
 219
 220    kz =  kde - kds + 1
 221 !
 222 ! 4.1 Get the square-root of the variance tuning factor:
 223
 224    n = 1+(it-1)*3
 225    if ( as1(n) <= 1.0E-7 .or. as2(n) <= 1.0E-7 .or. as3(n) <= 1.0E-7 .or. as4(n) <= 1.0E-7 .or. as5(n) <= 1.0E-7 ) then
 226      Write (i_char, '(i8)') it
 227      call da_error(__FILE__,__LINE__, &
 228                     (/"Missing or invalid AS1 ~ AS5 value for IT = "//i_char/))
 229    endif
 230
 231    as(1)=sqrt(as1(1+(it-1)*3))
 232    as(2)=sqrt(as2(1+(it-1)*3))
 233    as(3)=sqrt(as3(1+(it-1)*3))
 234    as(4)=sqrt(as4(1+(it-1)*3))
 235    as(5)=sqrt(as5(1+(it-1)*3))
 236
 237 ! 4.2 Scale the horizontal scale in unit of grid-point:
 238
 239    s2u= 1./grid%xb%ds
 240    hwll=hwll*s2u
 241    hwllp=hwllp*s2u
 242
 243 ! 4.3 Normalize the covariance for psi, chi, t, and rh:
 244
 245    !$OMP PARALLEL DO &
 246    !$OMP PRIVATE (ij, n, j, i, vv, k)
 247    do ij = 1, grid%num_tiles
 248    do n=1,4
 249
 250       do j=grid%j_start(ij), grid%j_end(ij)
 251          do i=its,ite
 252      ! Initialize the matrix vv:
 253             vv=0.
 254             do k=kts,kte
 255                vv(k,k)=1.
 256             enddo
 257
 258      ! Apply pseudo-Gaussian "right-left" smoother in vertical with
 259      ! the vertical scale length be%vz:
 260             call da_rfz0(vv,kz,kz,be%ndeg,&
 261                          be%vz(kts:kte,i,j,n),be%be,be%table,be%nta,be%swidth)
 262
 263      ! Normalize the covariance for psi, chi, t, and rh:
 264             do k=kts,kte
 265                be % corz(i,j,k,n)=be % corz(i,j,k,n)*as(n) &
 266                                   *samp(1)/hwll(i,j,k,n)/vv(k,k)/global_fac(i,j)
 267             enddo
 268
 269          enddo
 270       enddo
 271    enddo
 272    enddo
 273    !$OMP END PARALLEL DO
 274
 275 ! 4.4 Normalize the covariance for Psfc:
 276
 277 !    xsum(1) = sum (be%corp(its:ite,jts:jte)*be%corp(its:ite,jts:jte))
 278 !    call da_proc_sum_real(xsum)
 279 !    if (rootproc) &
 280 !    write (*,'("BEFORE: ",a,2x,"sum^2=",e20.12)') 'PSFC_u', xsum(1)
 281
 282 !   if (rootproc) then
 283 !   xsum(1) = sum (hwllp*hwllp)
 284 !   write (*,'(9x,a,2x,"sum^2=",e20.12)') 'HORIZONTAL PS_u SCALES:', xsum(1)
 285 !   endif
 286 !   write (*,'("as(5)=",f15.5,2x,"samp=",f15.5)') as(5), samp
 287
 288     be % corp(its:ite,jts:jte)=be % corp(its:ite,jts:jte)*as(5) &
 289          *samp(1)/hwllp(its:ite,jts:jte)/global_fac(its:ite,jts:jte)
 290
 291 !    xsum(1) = sum (be%corp(its:ite,jts:jte)*be%corp(its:ite,jts:jte))
 292 !    call da_proc_sum_real(xsum)
 293 !    if (rootproc) &
 294 !    write (*,'(" AFTER: ",a,2x,"sum^2=",e20.12)') 'PSFC_u', xsum(1)
 295
 296     write(mesg,*) 'Normalized covariance for Psfc: sum(be%corp*be%corp)=', &
 297                   sum(be%corp*be%corp)
 298     call wrf_debug ( 1 , mesg )
 299 !
 300 ! 5, Assign the inverse of scale length fields for recursive filter:
 301
 302    ic = 1 + ( ide- ids)/2
 303    jc = 1 + ( jde- jds)/2
 304
 305 !   if (rootproc) &
 306 !   write (*,'("i=",i4,2x,"jc=",i4)') ic, jc
 307 !
 308 ! 5.1 allocate the arrays for be scales: y-direction and x-direction:
 309
 310     ALLOCATE ( be % sljpy ( grid%xp%ipsy: grid%xp%ipey, grid%xp%jpsy: grid%xp%jpey) )
 311     ALLOCATE ( be % sljy ( grid%xp%ipsy: grid%xp%ipey, grid%xp%jpsy: grid%xp%jpey, grid%xp%kpsy: grid%xp%kpey,1:4) )
 312
 313     ALLOCATE ( be % slipx ( grid%xp%ipsx: grid%xp%ipex, grid%xp%jpsx: grid%xp%jpex) )
 314     ALLOCATE ( be % slix ( grid%xp%ipsx: grid%xp%ipex, grid%xp%jpsx: grid%xp%jpex, grid%xp%kpsx: grid%xp%kpex,1:4) )
 315
 316 ! 5.2 Y-direction:
 317
 318 ! 5.2.1 3-D fields: psi, chi, t, and rh:
 319
 320     do n=1,4
 321        do k= grid%xp%kpsy, grid%xp%kpey
 322           do j= grid%xp%jpsy, grid%xp%jpey
 323              do i= grid%xp%ipsy, grid%xp%ipey
 324                 be%sljy(i,j,k,n)=1./global_fac(i,j)/hwll(i,j,k,n)
 325              enddo
 326           enddo
 327        enddo
 328     enddo
 329 !    write (*,'("a, be%sljy:sum^2=",e20.12)') sum(be%sljy*be%sljy)
 330
 331     ! Above level ku, the sljy fields are set to a constant
 332     ! for psi and chi, i.e. homogenous:
 333     do n=1,2
 334        do k=max(ku(1), grid%xp%kpsy), grid%xp%kpey
 335           slim=1./global_fac(ic,jc)/hwll(ic,jc,k,n)
 336           do j= grid%xp%jpsy, grid%xp%jpey
 337              do i= grid%xp%ipsy, grid%xp%ipey
 338                 be%sljy(i,j,k,n)=slim
 339              enddo
 340           enddo
 341        enddo
 342     enddo
 343 !    write (*,'("b, be%sljy:sum^2=",e20.12)') sum(be%sljy*be%sljy)
 344
 345 ! 5.2.2 2-D field: Psfc:
 346     do j= grid%xp%jpsy, grid%xp%jpey
 347        do i= grid%xp%ipsy, grid%xp%ipey
 348           be%sljpy(i,j)=1./global_fac(i,j)/hwllp(i,j)
 349        enddo
 350     enddo
 351 !    write (*,'("   be%sljpy:sum^2=",e20.12)') sum(be%sljpy*be%sljpy)
 352
 353 ! 5.3 X-direction:
 354
 355 ! 5.3.1 3-D fields: psi, chi, t, and rh:
 356
 357     do n=1,4
 358       do k= grid%xp%kpsx, grid%xp%kpex
 359          do j= grid%xp%jpsx, grid%xp%jpex
 360             do i= grid%xp%ipsx, grid%xp%ipex
 361                be%slix(i,j,k,n)=1./global_fac(i,j)/hwll(i,j,k,n)
 362             enddo
 363          enddo
 364       enddo
 365     enddo
 366 !    write (*,'("a, be%slix:sum^2=",e20.12)') sum(be%slix*be%slix)
 367
 368    ! Above level ku, the sljy fields are set to a constant
 369    ! for psi and chi, i.e. homogenous:
 370     do n=1,2
 371       do k=max(ku(1), grid%xp%kpsx), grid%xp%kpex
 372          slim=1./global_fac(ic,jc)/hwll(ic,jc,k,n)
 373          do j= grid%xp%jpsx, grid%xp%jpex
 374             do i= grid%xp%ipsx, grid%xp%ipex
 375                be%slix(i,j,k,n)=slim
 376             enddo
 377          enddo
 378       enddo
 379     enddo
 380 !    write (*,'("b, be%slix:sum^2=",e20.12)') sum(be%slix*be%slix)
 381
 382 ! 5.3.2 2-D field: Psfc:
 383     do j= grid%xp%jpsx, grid%xp%jpex
 384       do i= grid%xp%ipsx, grid%xp%ipex
 385          be%slipx(i,j)=1./global_fac(i,j)/hwllp(i,j)
 386       enddo
 387     enddo
 388 !    write (*,'("   be%slipx:sum^2=",e20.12)') sum(be%slipx*be%slipx)
 389
 390 end subroutine da_scale_background_errors_cv3