var/da/da_update_bc/da_update_bc.f90

   1 program da_update_bc
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: update BC file from wrfvar output.
   5    ! current version reads only wrf-netcdf file format
   6    !
   7    ! Y.-R. Guo, 03/18/2008:
   8    !   1) Fixed the bug for low_bdy_only;
   9    !   2) Introducing another namelist variable: update_lsm
  10    !      update_lsm = .true. --- The LSM predicted variables:
  11    !                         TSLB, SMOIS, SNOW, SH2O, RHOSN, CANWAT, SNOWH
  12    !                              will be updated based on wrf_input file
  13    !                 = .false. -- no updated, default.
  14    !
  15    !-----------------------------------------------------------------------
  16
  17    use da_netcdf_interface, only : da_get_var_3d_real_cdf, &
  18       da_put_var_3d_real_cdf, da_get_dims_cdf, da_put_var_2d_real_cdf, &
  19       da_get_var_2d_real_cdf, da_get_var_2d_int_cdf, da_get_bdytimestr_cdf, &
  20       da_get_times_cdf, da_get_bdyfrq, stderr, stdout, da_put_var_2d_int_cdf, &
  21       da_get_var_1d_real_cdf, da_get_gl_att_int_cdf
  22
  23    use da_module_couple_uv, only : da_couple_uv
  24
  25    implicit none
  26
  27    include 'netcdf.inc'
  28
  29    integer, parameter :: max_3d_variables = 20, &
  30                          max_2d_variables = 25
  31
  32    character(len=512) :: da_file,      &
  33                          da_file_02,   &
  34                          wrf_bdy_file, &
  35                          wrf_input
  36
  37    character(len=20) :: var_pref, var_name, vbt_name
  38
  39    character(len=20) :: var3d(max_3d_variables), &
  40                         varsf(max_2d_variables)
  41
  42    character(len=10), dimension(4) :: bdyname, tenname
  43
  44    ! for WRF hybrid coordinate
  45    integer           :: nlevf, nlevh
  46    integer           :: hybrid_opt, use_theta_m
  47    real, allocatable :: c1f(:), c2f(:), c1h(:), c2h(:)
  48
  49    integer           :: ids, ide, jds, jde, kds, kde
  50    integer           :: num3d, num2d, ndims
  51    integer           :: time_level
  52    integer           :: i,j,k,l,m,n
  53
  54    integer, dimension(4) :: dims
  55
  56    real, allocatable, dimension(:,:,:) :: tend3d, scnd3d, frst3d, full3d, full3d2
  57
  58    real, allocatable, dimension(:,:,:) :: u, v, u2, v2
  59
  60    real, allocatable, dimension(:,  :) :: mu, mub, msfu, msfv, msfm, &
  61                                           mu2, tend2d, scnd2d, frst2d, full2d
  62
  63    real, allocatable, dimension(:,  :) :: tsk, tsk_wrfvar
  64    real, allocatable, dimension(:,:)   :: snow, snowc, snowh
  65
  66    integer, allocatable, dimension(:,:) :: ivgtyp, full2dint
  67
  68    character(len=80), allocatable, dimension(:) :: times, &
  69                                                    thisbdytime, nextbdytime
  70
  71    integer :: east_end, north_end, io_status, cdfid, varid, domain_id, iswater
  72    integer :: iostatus(4)
  73
  74    logical :: debug, update_lateral_bdy, update_low_bdy, update_lsm, keep_tsk_wrf
  75    logical :: keep_snow_wrf, var4d_lbc
  76
  77    real :: bdyfrq
  78
  79    character(len=512) :: wrfvar_output_file    ! obsolete. Kept for backward compatibility
  80    logical            :: cycling, low_bdy_only ! obsolete. Kept for backward compatibility
  81
  82    integer, parameter :: namelist_unit = 7, &
  83                          ori_unit = 11, &
  84                          new_unit = 12
  85
  86    namelist /control_param/ da_file,      &
  87                             da_file_02,   &
  88                             wrf_bdy_file, &
  89                             wrf_input, domain_id, var4d_lbc, &
  90                             debug, update_lateral_bdy, update_low_bdy, update_lsm, &
  91                             keep_tsk_wrf, keep_snow_wrf, iswater, &
  92                             wrfvar_output_file, cycling, low_bdy_only
  93
  94    da_file            = 'wrfvar_output'
  95    da_file_02         = 'ana02'
  96    wrf_bdy_file       = 'wrfbdy_d01'
  97    wrf_input          = 'wrfinput_d01'
  98    domain_id          = 1
  99
 100    var4d_lbc          = .false.
 101    debug              = .false.
 102    update_lateral_bdy = .true.
 103    update_low_bdy     = .true.
 104    update_lsm         = .false.
 105    keep_tsk_wrf       = .true.
 106    keep_snow_wrf      = .true.
 107    iswater            = 16      ! USGS water index: 16, MODIS water index: 17
 108
 109    wrfvar_output_file = 'OBSOLETE'
 110    cycling            = .false.
 111    low_bdy_only       = .false.
 112
 113    !---------------------------------------------------------------------
 114    ! Read namelist
 115    !---------------------------------------------------------------------
 116    io_status = 0
 117
 118    open(unit = namelist_unit, file = 'parame.in', &
 119           status = 'old' , access = 'sequential', &
 120           form   = 'formatted', action = 'read', &
 121           iostat = io_status)
 122
 123    if (io_status /= 0) then
 124       write(unit=stdout,fmt=*) 'Error to open namelist file: parame.in.'
 125       write(unit=stdout,fmt=*) 'Will work for updating lateral boundary only.'
 126    else
 127       read(unit=namelist_unit, nml = control_param , iostat = io_status)
 128
 129       if (io_status /= 0) then
 130          write(unit=stdout,fmt=*) 'Error to read control_param. Stopped.'
 131          stop
 132       end if
 133
 134       ! deal with the old namelist
 135       if ( index(wrfvar_output_file, 'OBSOLETE') <= 0 ) then
 136          ! wrfvar_output_file is set in the user's parame.in
 137          ! reset the settings
 138          da_file = wrfvar_output_file
 139          if ( domain_id > 1 ) then
 140             low_bdy_only = .true.
 141          end if
 142          if ( cycling .and. domain_id == 1 ) then
 143             update_lateral_bdy = .true.
 144             update_low_bdy     = .true.
 145          else
 146             if ( low_bdy_only ) then
 147                update_lateral_bdy = .false.
 148                update_low_bdy     = .true.
 149             else
 150                update_lateral_bdy = .true.
 151                update_low_bdy     = .false.
 152             end if
 153          end if
 154       end if
 155
 156       WRITE(unit=stdout, fmt='(2a)') &
 157            'da_file       = ', trim(da_file), &
 158            'da_file_02    = ', trim(da_file_02), &
 159            'wrf_bdy_file  = ', trim(wrf_bdy_file), &
 160            'wrf_input     = ', trim(wrf_input)
 161
 162       WRITE(unit=stdout, fmt='(2(a, L10))')             &
 163            'update_lateral_bdy = ', update_lateral_bdy, &
 164            'update_low_bdy     = ', update_low_bdy
 165
 166       if ( update_lsm ) keep_snow_wrf = .false.
 167
 168       close(unit=namelist_unit)
 169    end if
 170
 171    ! 3D need update
 172    num3d=5
 173    var3d(1)='U'
 174    var3d(2)='V'
 175    var3d(3)='T'
 176    var3d(4)='PH'
 177    var3d(5)='QVAPOR'
 178 !  var3d(6)='W'
 179
 180    ! 2D need update
 181    num2d=23
 182    varsf(1)='MUB'
 183    varsf(2)='MU'
 184    varsf(3)='MAPFAC_U'
 185    varsf(4)='MAPFAC_V'
 186    varsf(5)='MAPFAC_M'
 187    varsf(6)='TMN'
 188    varsf(7)='SST'
 189    varsf(8)='TSK'
 190    varsf(9)='VEGFRA'
 191    varsf(10)='ALBBCK'
 192    varsf(11)='TSLB'
 193    varsf(12)='SMOIS'
 194    varsf(13)='SNOW'
 195    varsf(14)='SEAICE'
 196    varsf(15)='SH2O'
 197    varsf(16)='CANWAT'
 198    varsf(17)='RHOSN'
 199    varsf(18)='SNOWH'
 200    varsf(19)='LANDMASK'
 201    varsf(20)='IVGTYP'
 202    varsf(21)='ISLTYP'
 203    varsf(22)='SNOWC'
 204    varsf(23)='XLAND'
 205
 206    if ( domain_id > 1 ) then
 207       write(unit=stdout, fmt='(a,i2)') 'Nested domain ID=',domain_id
 208       write(unit=stdout, fmt='(a)') &
 209         'No wrfbdy file needed, only low boundary need to be updated.'
 210       if ( update_lateral_bdy ) then
 211          write(unit=stdout, fmt='(a)') &
 212             'Re-setting update_lateral_bdy to be false for nested domain.'
 213          update_lateral_bdy = .false.
 214       end if
 215       update_low_bdy     = .true.
 216    end if
 217
 218    if ( update_lateral_bdy ) then
 219    ! First, the boundary times
 220    call da_get_dims_cdf(wrf_bdy_file, 'Times', dims, ndims, debug)
 221
 222    if (debug) then
 223       write(unit=stdout, fmt='(a,i2,2x,a,4i6)') &
 224            'Times: ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 225    end if
 226
 227    time_level = dims(2)
 228
 229    if (time_level < 1) then
 230       write(unit=stdout, fmt='(a,i2/a)') &
 231            'time_level = ', time_level, &
 232            'We need at least one time-level BDY.'
 233       stop 'Wrong BDY file.'
 234    end if
 235
 236    allocate(times(dims(2)))
 237    allocate(thisbdytime(dims(2)))
 238    allocate(nextbdytime(dims(2)))
 239
 240    call da_get_times_cdf(wrf_bdy_file, times, dims(2), dims(2), debug)
 241
 242    call da_get_bdytimestr_cdf(wrf_bdy_file, 'thisbdytime', thisbdytime, dims(2), debug)
 243    call da_get_bdytimestr_cdf(wrf_bdy_file, 'nextbdytime', nextbdytime, dims(2), debug)
 244
 245    call da_get_bdyfrq(thisbdytime(1), nextbdytime(1), bdyfrq, debug)
 246
 247    if (debug) then
 248       do n=1, dims(2)
 249          write(unit=stdout, fmt='(3(a, i2, 2a,2x))') &
 250            '       times(', n, ')=', trim(times(n)), &
 251            'thisbdytime (', n, ')=', trim(thisbdytime(n)), &
 252            'nextbdytime (', n, ')=', trim(nextbdytime(n))
 253       end do
 254    end if
 255
 256    end if
 257
 258    east_end=0
 259    north_end=0
 260
 261    cdfid = ncopn(da_file, NCWRITE, io_status )
 262
 263    ! for WRF hybrid coordinate
 264
 265    call da_get_gl_att_int_cdf(da_file, 'BOTTOM-TOP_PATCH_END_STAG',   nlevf, debug, io_status)
 266    call da_get_gl_att_int_cdf(da_file, 'BOTTOM-TOP_PATCH_END_UNSTAG', nlevh, debug, io_status)
 267    if ( io_status /= NF_NOERR .or. nlevf <= 0 .or. nlevh <= 0 ) then
 268       write(unit=stdout,fmt=*) 'Error finding the number of vertical levels'
 269       stop
 270    else
 271       allocate ( c1f(nlevf) )
 272       allocate ( c2f(nlevf) )
 273       allocate ( c1h(nlevh) )
 274       allocate ( c2h(nlevh) )
 275    end if
 276
 277    ! initialize as use_theta_m = 0
 278    use_theta_m = 0
 279    call da_get_gl_att_int_cdf(da_file, 'USE_THETA_M', use_theta_m, debug, io_status)
 280    write(stdout,*) 'use_theta_m = ', use_theta_m
 281
 282    ! initialize as hybrid_opt = 0
 283    hybrid_opt = 0
 284    c1f(:) = 1.0
 285    c2f(:) = 0.0
 286    c1h(:) = 1.0
 287    c2h(:) = 0.0
 288
 289    call da_get_gl_att_int_cdf(da_file, 'HYBRID_OPT', hybrid_opt, debug, io_status)
 290    if ( io_status == NF_NOERR .and. hybrid_opt > 0 ) then
 291       ! when HYBRID_OPT is available from the file (V3.9 and later)
 292       ! use C1F,C2F,C1H,C2H from the file
 293       call da_get_var_1d_real_cdf( da_file, 'C1F', c1f, nlevf, 1, debug)
 294       call da_get_var_1d_real_cdf( da_file, 'C2F', c2f, nlevf, 1, debug)
 295       call da_get_var_1d_real_cdf( da_file, 'C1H', c1h, nlevh, 1, debug)
 296       call da_get_var_1d_real_cdf( da_file, 'C2H', c2h, nlevh, 1, debug)
 297    end if
 298
 299    write(stdout,*) 'hybrid_opt = ', hybrid_opt
 300    if ( debug ) then
 301       write(stdout,*) 'c1f = ', c1f
 302       write(stdout,*) 'c2f = ', c2f
 303       write(stdout,*) 'c1h = ', c1h
 304       write(stdout,*) 'c2h = ', c2f
 305    end if
 306
 307    ! For 2D variables
 308    ! Get mu, mub, msfu, and msfv
 309
 310    do n=1,num2d
 311
 312       io_status = nf_inq_varid(cdfid, trim(varsf(n)), varid)
 313       if (io_status /= 0 ) then
 314          print '(/"N=",i2," io_status=",i5,5x,"VAR=",a,a)', &
 315                    n, io_status, trim(varsf(n)), " does not exist"
 316          cycle
 317       endif
 318
 319       call da_get_dims_cdf( da_file, trim(varsf(n)), dims, &
 320          ndims, debug)
 321
 322       select case(trim(varsf(n)))
 323       case ('MU') ;
 324          if ( .not. update_lateral_bdy ) cycle
 325
 326          allocate(mu(dims(1), dims(2)))
 327
 328          call da_get_var_2d_real_cdf( da_file, &
 329             trim(varsf(n)), mu, dims(1), dims(2), 1, debug)
 330
 331          east_end=dims(1)+1
 332          north_end=dims(2)+1
 333
 334          if ( var4d_lbc ) then
 335             allocate(mu2(dims(1), dims(2)))
 336
 337             call da_get_var_2d_real_cdf( da_file_02, &
 338                trim(varsf(n)), mu2, dims(1), dims(2), 1, debug)
 339          end if
 340       case ('MUB') ;
 341          if ( .not. update_lateral_bdy ) cycle
 342
 343          allocate(mub(dims(1), dims(2)))
 344
 345          call da_get_var_2d_real_cdf( da_file, trim(varsf(n)), mub, &
 346                                    dims(1), dims(2), 1, debug)
 347       case ('MAPFAC_U') ;
 348          if ( .not. update_lateral_bdy ) cycle
 349
 350          allocate(msfu(dims(1), dims(2)))
 351
 352          call da_get_var_2d_real_cdf( da_file, trim(varsf(n)), msfu, &
 353                                    dims(1), dims(2), 1, debug)
 354       case ('MAPFAC_V') ;
 355          if ( .not. update_lateral_bdy ) cycle
 356
 357          allocate(msfv(dims(1), dims(2)))
 358
 359          call da_get_var_2d_real_cdf( da_file, trim(varsf(n)), msfv, &
 360                                    dims(1), dims(2), 1, debug)
 361       case ('MAPFAC_M') ;
 362          if ( .not. update_lateral_bdy ) cycle
 363
 364          allocate(msfm(dims(1), dims(2)))
 365
 366          call da_get_var_2d_real_cdf( da_file, trim(varsf(n)), msfm, &
 367                                    dims(1), dims(2), 1, debug)
 368       case ('TSK') ;
 369          if ( .not. update_low_bdy ) cycle
 370
 371          allocate(tsk(dims(1), dims(2)))
 372          allocate(tsk_wrfvar(dims(1), dims(2)))
 373          allocate(ivgtyp(dims(1), dims(2)))
 374
 375          call da_get_var_2d_real_cdf( wrf_input, trim(varsf(n)), tsk, &
 376                                    dims(1), dims(2), 1, debug)
 377
 378          if ( keep_tsk_wrf ) then
 379             call da_get_var_2d_real_cdf( da_file, trim(varsf(n)), tsk_wrfvar, &
 380                                       dims(1), dims(2), 1, debug)
 381             !hcl call da_get_var_2d_int_cdf( da_file, 'IVGTYP', ivgtyp, &
 382             call da_get_var_2d_int_cdf( wrf_input, 'IVGTYP', ivgtyp, &
 383                                       dims(1), dims(2), 1, debug)
 384             ! update TSK.
 385             do j=1,dims(2)
 386                do i=1,dims(1)
 387                   if (ivgtyp(i,j) /= iswater)  tsk(i,j)=tsk_wrfvar(i,j)
 388                end do
 389             end do
 390          end if
 391
 392             call da_put_var_2d_real_cdf( da_file, trim(varsf(n)), tsk, &
 393                                       dims(1), dims(2), 1, debug)
 394             deallocate(tsk)
 395             deallocate(ivgtyp)
 396             deallocate(tsk_wrfvar)
 397
 398          !hcl case ('TMN', 'SST', 'VEGFRA', 'ALBBCK', 'SEAICE') ;
 399          case ('TMN', 'SST', 'VEGFRA', 'ALBBCK', 'SEAICE', 'LANDMASK', 'XLAND') ;
 400             if ( .not. update_low_bdy ) cycle
 401
 402             allocate(full2d(dims(1), dims(2)))
 403
 404             call da_get_var_2d_real_cdf( wrf_input, trim(varsf(n)), full2d, &
 405                                       dims(1), dims(2), 1, debug)
 406
 407             call da_put_var_2d_real_cdf( da_file, trim(varsf(n)), full2d, &
 408                                       dims(1), dims(2), 1, debug)
 409             deallocate(full2d)
 410
 411          case ('IVGTYP', 'ISLTYP') ;  !hcl add
 412             if ( .not. update_low_bdy ) cycle
 413
 414             allocate(full2dint(dims(1), dims(2)))
 415
 416             call da_get_var_2d_int_cdf( wrf_input, trim(varsf(n)), full2dint, &
 417                                       dims(1), dims(2), 1, debug)
 418
 419             call da_put_var_2d_int_cdf( da_file, trim(varsf(n)), full2dint, &
 420                                       dims(1), dims(2), 1, debug)
 421             deallocate(full2dint)
 422
 423          case ('SNOW', 'RHOSN', 'SNOWH', 'SNOWC') ;
 424             if ( (.not. update_lsm) .and. (.not. update_low_bdy) ) cycle
 425             if ( keep_snow_wrf ) cycle
 426                allocate(full2d(dims(1), dims(2)))
 427
 428                call da_get_var_2d_real_cdf( wrf_input, trim(varsf(n)), full2d, &
 429                                       dims(1), dims(2), 1, debug )
 430
 431                call da_put_var_2d_real_cdf( da_file, trim(varsf(n)), full2d, &
 432                                       dims(1), dims(2), 1, debug )
 433                deallocate(full2d)
 434
 435          case ('CANWAT') ;
 436             if ( .not. update_lsm ) cycle
 437                allocate(full2d(dims(1), dims(2)))
 438
 439                call da_get_var_2d_real_cdf( wrf_input, trim(varsf(n)), full2d, &
 440                                       dims(1), dims(2), 1, debug )
 441 !               print *,"sum(full2d^2)=", sum(full2d*full2d)
 442
 443                call da_put_var_2d_real_cdf( da_file, trim(varsf(n)), full2d, &
 444                                       dims(1), dims(2), 1, debug )
 445                deallocate(full2d)
 446
 447          case ('TSLB', 'SMOIS', 'SH2O') ;
 448             if( .not. update_lsm ) cycle
 449                allocate(full3d(dims(1), dims(2), dims(3)))
 450
 451                call da_get_var_3d_real_cdf( wrf_input, trim(varsf(n)), full3d, &
 452                                       dims(1), dims(2), dims(3), 1, debug )
 453 !               print *,"sum(full3d^2)=", sum(full3d*full3d)
 454
 455                call da_put_var_3d_real_cdf( da_file, trim(varsf(n)), full3d, &
 456                                       dims(1), dims(2), dims(3), 1, debug )
 457                deallocate(full3d)
 458
 459          case default ;
 460             write(unit=stdout,fmt=*) 'It is impossible here. varsf(n)=', trim(varsf(n))
 461       end select
 462    end do
 463
 464    ! check for snow over water
 465    iostatus(1) = nf_inq_varid(cdfid, 'IVGTYP', varid)
 466    iostatus(2) = nf_inq_varid(cdfid, 'SNOW',   varid)
 467    iostatus(3) = nf_inq_varid(cdfid, 'SNOWC',  varid)
 468    iostatus(4) = nf_inq_varid(cdfid, 'SNOWH',  varid)
 469    if ( iostatus(1) == 0 ) then
 470       allocate(snow(dims(1), dims(2)))
 471       allocate(snowc(dims(1), dims(2)))
 472       allocate(snowh(dims(1), dims(2)))
 473       allocate(ivgtyp(dims(1), dims(2)))
 474       if ( iostatus(1) == 0 ) then
 475          call da_get_var_2d_int_cdf( da_file, 'IVGTYP', ivgtyp,    &
 476                                dims(1), dims(2), 1, debug)
 477       end if
 478       if ( iostatus(2) == 0 ) then
 479          call da_get_var_2d_real_cdf( da_file, 'SNOW',    snow,     &
 480                                 dims(1), dims(2), 1, debug)
 481       end if
 482       if ( iostatus(3) == 0 ) then
 483          call da_get_var_2d_real_cdf( da_file, 'SNOWC',  snowc,     &
 484                                 dims(1), dims(2), 1, debug)
 485       end if
 486       if ( iostatus(4) == 0 ) then
 487          call da_get_var_2d_real_cdf( da_file, 'SNOWH',  snowh,     &
 488                                 dims(1), dims(2), 1, debug)
 489       end if
 490       if ( iostatus(2) == 0 ) then
 491          do j = 1, dims(2)
 492             do i = 1, dims(1)
 493                if (ivgtyp(i,j) == iswater)  then
 494                   if ( snow(i,j) > 0.0 ) then
 495                      write(unit=stdout,fmt=*) 'Remove snow over water at i, j = ', i, j
 496                      if ( iostatus(2) == 0 ) snow(i,j)  = 0.0
 497                      if ( iostatus(3) == 0 ) snowc(i,j) = 0.0
 498                      if ( iostatus(4) == 0 ) snowh(i,j) = 0.0
 499                   end if
 500                end if
 501             end do
 502          end do
 503       end if
 504       if ( iostatus(2) == 0 ) then
 505          call da_put_var_2d_real_cdf( da_file, 'SNOW',   snow, &
 506                                 dims(1), dims(2), 1, debug)
 507       end if
 508       if ( iostatus(3) == 0 ) then
 509          call da_put_var_2d_real_cdf( da_file, 'SNOWC',  snowc, &
 510                                 dims(1), dims(2), 1, debug)
 511       end if
 512       if ( iostatus(4) == 0 ) then
 513          call da_put_var_2d_real_cdf( da_file, 'SNOWH',  snowh, &
 514                                 dims(1), dims(2), 1, debug)
 515       end if
 516       deallocate(snow)
 517       deallocate(snowc)
 518       deallocate(snowh)
 519       deallocate(ivgtyp)
 520    end if
 521
 522  if ( update_lateral_bdy ) then
 523
 524    if (east_end < 1 .or. north_end < 1) then
 525       write(unit=stdout, fmt='(a)') 'Wrong data for Boundary.'
 526       stop
 527    end if
 528
 529    write(unit=stdout,fmt='(/a/)') 'Processing the lateral boundary condition:'
 530
 531    ! boundary variables
 532    bdyname(1)='_BXS'
 533    bdyname(2)='_BXE'
 534    bdyname(3)='_BYS'
 535    bdyname(4)='_BYE'
 536
 537    ! boundary tendancy variables
 538    tenname(1)='_BTXS'
 539    tenname(2)='_BTXE'
 540    tenname(3)='_BTYS'
 541    tenname(4)='_BTYE'
 542
 543    do m=1,4
 544       var_name='MU' // trim(bdyname(m))
 545       vbt_name='MU' // trim(tenname(m))
 546
 547       call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
 548
 549       allocate(frst2d(dims(1), dims(2)))
 550       allocate(scnd2d(dims(1), dims(2)))
 551       allocate(tend2d(dims(1), dims(2)))
 552
 553       ! Get variable at second time level
 554       if ( .not. var4d_lbc ) then
 555          if (time_level > 1) then
 556             call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), scnd2d, &
 557                                       dims(1), dims(2), 2, debug)
 558          else
 559             call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
 560                                       dims(1), dims(2), 1, debug)
 561             call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
 562                                       dims(1), dims(2), 1, debug)
 563          end if
 564       end if
 565
 566       if (debug) then
 567          write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
 568               'No.', m, 'Variable: ', trim(vbt_name), &
 569               'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 570
 571          call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
 572                                    dims(1), dims(2), 1, debug)
 573
 574          write(unit=ori_unit, fmt='(a, 10i12)') &
 575               ' old ', (i, i=1,dims(2))
 576          do j=1,dims(1)
 577             write(unit=ori_unit, fmt='(i4, 1x, 10e20.7)') &
 578                   j, (tend2d(j,i), i=1,dims(2))
 579          end do
 580       end if
 581
 582       ! calculate variable at first time level
 583       select case(m)
 584       case (1) ;             ! West boundary
 585          do l=1,dims(2)
 586             do j=1,dims(1)
 587                if (time_level < 2 .and. .not. var4d_lbc) &
 588                   scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
 589                if (var4d_lbc) scnd2d(j,l)=mu2(l,j)
 590                frst2d(j,l)=mu(l,j)
 591             end do
 592          end do
 593       case (2) ;             ! East boundary
 594          do l=1,dims(2)
 595             do j=1,dims(1)
 596                if (time_level < 2 .and. .not. var4d_lbc) &
 597                   scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
 598                if (var4d_lbc) scnd2d(j,l)=mu2(east_end-l,j)
 599                frst2d(j,l)=mu(east_end-l,j)
 600             end do
 601          end do
 602       case (3) ;             ! South boundary
 603          do l=1,dims(2)
 604             do i=1,dims(1)
 605                if (time_level < 2 .and. .not. var4d_lbc) &
 606                   scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
 607                if (var4d_lbc) scnd2d(i,l)=mu2(i,l)
 608                frst2d(i,l)=mu(i,l)
 609             end do
 610          end do
 611       case (4) ;             ! North boundary
 612          do l=1,dims(2)
 613             do i=1,dims(1)
 614                if (time_level < 2 .and. .not. var4d_lbc) &
 615                   scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
 616                if (var4d_lbc) scnd2d(i,l)=mu2(i,north_end-l)
 617                frst2d(i,l)=mu(i,north_end-l)
 618             end do
 619          end do
 620       case default ;
 621          write(unit=stdout,fmt=*) 'It is impossible here. mu, m=', m
 622       end select
 623
 624       ! calculate new tendancy
 625       do l=1,dims(2)
 626          do i=1,dims(1)
 627             tend2d(i,l)=(scnd2d(i,l)-frst2d(i,l))/bdyfrq
 628          end do
 629       end do
 630
 631       if (debug) then
 632          write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
 633               'No.', m, 'Variable: ', trim(vbt_name), &
 634               'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 635
 636          write(unit=new_unit, fmt='(a, 10i12)') &
 637               ' new ', (i, i=1,dims(2))
 638
 639          do j=1,dims(1)
 640             write(unit=new_unit, fmt='(i4, 1x, 10e20.7)') &
 641                   j, (tend2d(j,i), i=1,dims(2))
 642          end do
 643       end if
 644
 645       ! output new variable at first time level
 646       call da_put_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
 647                                 dims(1), dims(2), 1, debug)
 648       ! output new tendancy
 649       call da_put_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
 650                                 dims(1), dims(2), 1, debug)
 651
 652       deallocate(frst2d)
 653       deallocate(scnd2d)
 654       deallocate(tend2d)
 655    end do
 656
 657    !---------------------------------------------------------------------
 658    ! For 3D variables
 659
 660    ! Get U
 661    call da_get_dims_cdf( da_file, 'U', dims, ndims, debug)
 662
 663    ! call da_get_att_cdf( da_file, 'U', debug)
 664
 665    allocate(u(dims(1), dims(2), dims(3)))
 666
 667    ids=1
 668    ide=dims(1)-1
 669    jds=1
 670    jde=dims(2)
 671    kds=1
 672    kde=dims(3)
 673
 674    call da_get_var_3d_real_cdf( da_file, 'U', u, &
 675                              dims(1), dims(2), dims(3), 1, debug)
 676    if ( var4d_lbc ) then
 677       allocate(u2(dims(1), dims(2), dims(3)))
 678       call da_get_var_3d_real_cdf( da_file_02, 'U', u2, &
 679                                 dims(1), dims(2), dims(3), 1, debug)
 680    end if
 681
 682    ! do j=1,dims(2)
 683    !    write(unit=stdout, fmt='(2(a,i5), a, f12.8)') &
 684    !       'u(', dims(1), ',', j, ',1)=', u(dims(1),j,1)
 685    ! end do
 686
 687    ! Get V
 688    call da_get_dims_cdf( da_file, 'V', dims, ndims, debug)
 689
 690    ! call da_get_att_cdf( da_file, 'V', debug)
 691
 692    allocate(v(dims(1), dims(2), dims(3)))
 693
 694    call da_get_var_3d_real_cdf( da_file, 'V', v, &
 695                              dims(1), dims(2), dims(3), 1, debug)
 696    if ( var4d_lbc ) then
 697       allocate(v2(dims(1), dims(2), dims(3)))
 698       call da_get_var_3d_real_cdf( da_file_02, 'V', v2, &
 699                                 dims(1), dims(2), dims(3), 1, debug)
 700    end if
 701
 702    ! do i=1,dims(1)
 703    !    write(unit=stdout, fmt='(2(a,i5), a, f12.8)') &
 704    !       'v(', i, ',', dims(2), ',1)=', v(i,dims(2),1)
 705    ! end do
 706
 707    if (debug) then
 708       write(unit=stdout, fmt='(a,e20.12,4x)') &
 709            'Before couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
 710            'Before couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
 711    end if
 712
 713    !---------------------------------------------------------------------
 714    ! Couple u, v.
 715    call da_couple_uv ( u, v, mu, mub, msfu, msfv, c1h, c2h, ids, ide, jds, jde, kds, kde)
 716    if ( var4d_lbc ) then
 717       call da_couple_uv ( u2, v2, mu2, mub, msfu, msfv, c1h, c2h, ids, ide, jds, jde, kds, kde)
 718    end if
 719
 720    if (debug) then
 721       write(unit=stdout, fmt='(a,e20.12,4x)') &
 722            'After  couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
 723            'After  couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
 724    end if
 725
 726    !---------------------------------------------------------------------
 727    !For 3D variables
 728
 729    do n=1,num3d
 730       write(unit=stdout, fmt='(a, i3, 2a)') 'Processing: var3d(', n, ')=', trim(var3d(n))
 731
 732       call da_get_dims_cdf( da_file, trim(var3d(n)), dims, ndims, debug)
 733
 734       allocate(full3d(dims(1), dims(2), dims(3)))
 735       if ( var4d_lbc ) allocate(full3d2(dims(1), dims(2), dims(3)))
 736
 737       east_end=dims(1)+1
 738       north_end=dims(2)+1
 739
 740       select case(trim(var3d(n)))
 741       case ('U') ;           ! U
 742          ! var_pref='R' // trim(var3d(n))
 743          var_pref=trim(var3d(n))
 744          full3d(:,:,:)=u(:,:,:)
 745          if ( var4d_lbc ) full3d2(:,:,:)=u2(:,:,:)
 746       case ('V') ;           ! V
 747          ! var_pref='R' // trim(var3d(n))
 748          var_pref=trim(var3d(n))
 749          full3d(:,:,:)=v(:,:,:)
 750          if ( var4d_lbc ) full3d2(:,:,:)=v2(:,:,:)
 751       case ('W') ;
 752          ! var_pref = 'R' // trim(var3d(n))
 753          var_pref = trim(var3d(n))
 754
 755          call da_get_var_3d_real_cdf( da_file, trim(var3d(n)), &
 756             full3d, dims(1), dims(2), dims(3), 1, debug)
 757          if ( var4d_lbc ) &
 758             call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
 759                full3d2, dims(1), dims(2), dims(3), 1, debug)
 760
 761          if (debug) then
 762             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 763                  'Before couple Sample ', trim(var3d(n)), &
 764                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 765          end if
 766
 767          do k=1,dims(3)
 768             do j=1,dims(2)
 769                do i=1,dims(1)
 770                   full3d(i,j,k)=full3d(i,j,k)*(c1f(k)*(mu(i,j)+mub(i,j))+c2f(k))/msfm(i,j)
 771                   if ( var4d_lbc ) full3d2(i,j,k)=full3d2(i,j,k)*(c1f(k)*(mu2(i,j)+mub(i,j))+c2f(k))/msfm(i,j)
 772                end do
 773             end do
 774          end do
 775
 776          if (debug) then
 777             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 778                  'After  couple Sample ', trim(var3d(n)), &
 779                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 780          end if
 781       case ('T', 'PH') ;
 782          var_pref=trim(var3d(n))
 783
 784          if ( use_theta_m > 0 .and. trim(var3d(n)) == 'T' ) then
 785            call da_get_var_3d_real_cdf( da_file, 'THM', &
 786               full3d, dims(1), dims(2), dims(3), 1, debug)
 787            if ( var4d_lbc ) &
 788               call da_get_var_3d_real_cdf( da_file_02, 'THM', &
 789                  full3d2, dims(1), dims(2), dims(3), 1, debug)
 790          else
 791            call da_get_var_3d_real_cdf( da_file, trim(var3d(n)), &
 792               full3d, dims(1), dims(2), dims(3), 1, debug)
 793            if ( var4d_lbc ) &
 794               call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
 795                  full3d2, dims(1), dims(2), dims(3), 1, debug)
 796          end if
 797
 798          if (debug) then
 799             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 800                  'Before couple Sample ', trim(var3d(n)), &
 801                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 802          end if
 803
 804          if ( dims(3) == nlevf ) then
 805             ! applying c1f and c2f for mu
 806             do k=1,dims(3)
 807                do j=1,dims(2)
 808                   do i=1,dims(1)
 809                      full3d(i,j,k)=full3d(i,j,k)*(c1f(k)*(mu(i,j)+mub(i,j))+c2f(k))
 810                      if ( var4d_lbc ) full3d2(i,j,k)=full3d2(i,j,k)*(c1f(k)*(mu2(i,j)+mub(i,j))+c2f(k))
 811                   end do
 812                end do
 813             end do
 814          else if ( dims(3) == nlevh ) then
 815             ! applying c1h and c2h for mu
 816             do k=1,dims(3)
 817                do j=1,dims(2)
 818                   do i=1,dims(1)
 819                      full3d(i,j,k)=full3d(i,j,k)*(c1h(k)*(mu(i,j)+mub(i,j))+c2h(k))
 820                      if ( var4d_lbc ) full3d2(i,j,k)=full3d2(i,j,k)*(c1h(k)*(mu2(i,j)+mub(i,j))+c2h(k))
 821                   end do
 822                end do
 823             end do
 824          else
 825             write(unit=stdout,fmt=*) 'Error finding matched dimensions.'
 826             stop
 827          end if
 828
 829          if (debug) then
 830             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 831                  'After  couple Sample ', trim(var3d(n)), &
 832                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 833          end if
 834       case ('QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW', 'QGRAUP') ;
 835          ! var_pref='R' // var3d(n)(1:2)
 836          ! var_pref=var3d(n)(1:2)
 837          var_pref=var3d(n)
 838
 839          call da_get_var_3d_real_cdf( da_file, trim(var3d(n)), &
 840             full3d, dims(1), dims(2), dims(3), 1, debug)
 841          if ( var4d_lbc ) &
 842             call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
 843                full3d2, dims(1), dims(2), dims(3), 1, debug)
 844
 845          if (debug) then
 846             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 847                  'Before couple Sample ', trim(var3d(n)), &
 848                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 849          end if
 850
 851          do k=1,dims(3)
 852             do j=1,dims(2)
 853                do i=1,dims(1)
 854                   full3d(i,j,k)=full3d(i,j,k)*(c1h(k)*(mu(i,j)+mub(i,j))+c2h(k))
 855                   if ( var4d_lbc ) full3d2(i,j,k)=full3d2(i,j,k)*(c1h(k)*(mu2(i,j)+mub(i,j))+c2h(k))
 856                end do
 857             end do
 858          end do
 859
 860          if (debug) then
 861             write(unit=stdout, fmt='(3a,e20.12,4x)') &
 862                  'After  couple Sample ', trim(var3d(n)), &
 863                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
 864          end if
 865       case default ;
 866          write(unit=stdout,fmt=*) 'It is impossible here. var3d(', n, ')=', trim(var3d(n))
 867       end select
 868
 869       do m=1,4
 870          var_name=trim(var_pref) // trim(bdyname(m))
 871          vbt_name=trim(var_pref) // trim(tenname(m))
 872
 873          write(unit=stdout, fmt='(a, i3, 2a)') &
 874             'Processing: bdyname(', m, ')=', trim(var_name)
 875
 876          call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
 877
 878          allocate(frst3d(dims(1), dims(2), dims(3)))
 879          allocate(scnd3d(dims(1), dims(2), dims(3)))
 880          allocate(tend3d(dims(1), dims(2), dims(3)))
 881
 882          ! Get variable at second time level
 883          if ( .not. var4d_lbc ) then
 884             if (time_level > 1) then
 885                call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), scnd3d, &
 886                                          dims(1), dims(2), dims(3), 2, debug)
 887             else
 888                call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
 889                                          dims(1), dims(2), dims(3), 1, debug)
 890                call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
 891                                          dims(1), dims(2), dims(3), 1, debug)
 892             end if
 893          end if
 894
 895          if (debug) then
 896             write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
 897                  'No.', m, 'Variable: ', trim(vbt_name), &
 898                  'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 899
 900             call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
 901                                       dims(1), dims(2), dims(3), 1, debug)
 902
 903             write(unit=ori_unit, fmt='(a, 10i12)') &
 904                  ' old ', (i, i=1,dims(3))
 905             do j=1,dims(1)
 906                write(unit=ori_unit, fmt='(i4, 1x, 10e20.7)') &
 907                      j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
 908             end do
 909          end if
 910
 911          select case(trim(bdyname(m)))
 912          case ('_BXS') ;             ! West boundary
 913             do l=1,dims(3)
 914             do k=1,dims(2)
 915             do j=1,dims(1)
 916                if (time_level < 2 .and. .not. var4d_lbc) &
 917                scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
 918                if ( var4d_lbc ) scnd3d(j,k,l)=full3d2(l,j,k)
 919                frst3d(j,k,l)=full3d(l,j,k)
 920             end do
 921             end do
 922             end do
 923          case ('_BXE') ;             ! East boundary
 924             do l=1,dims(3)
 925             do k=1,dims(2)
 926             do j=1,dims(1)
 927                if (time_level < 2 .and. .not. var4d_lbc) &
 928                scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
 929                if ( var4d_lbc ) scnd3d(j,k,l)=full3d2(east_end-l,j,k)
 930                frst3d(j,k,l)=full3d(east_end-l,j,k)
 931             end do
 932             end do
 933             end do
 934          case ('_BYS') ;             ! South boundary
 935             do l=1,dims(3)
 936             do k=1,dims(2)
 937             do i=1,dims(1)
 938                if (time_level < 2 .and. .not. var4d_lbc) &
 939                scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
 940                if ( var4d_lbc )scnd3d(i,k,l)=full3d2(i,l,k)
 941                frst3d(i,k,l)=full3d(i,l,k)
 942             end do
 943             end do
 944             end do
 945          case ('_BYE') ;             ! North boundary
 946             do l=1,dims(3)
 947             do k=1,dims(2)
 948             do i=1,dims(1)
 949                if (time_level < 2 .and. .not. var4d_lbc) &
 950                scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
 951                if ( var4d_lbc ) scnd3d(i,k,l)=full3d2(i,north_end-l,k)
 952                frst3d(i,k,l)=full3d(i,north_end-l,k)
 953             end do
 954             end do
 955             end do
 956          case default ;
 957             write(unit=stdout,fmt=*) 'It is impossible here.'
 958             write(unit=stdout,fmt=*) 'bdyname(', m, ')=', trim(bdyname(m))
 959             stop
 960          end select
 961
 962          write(unit=stdout, fmt='(a, i3, 2a)') &
 963             'cal. tend: bdyname(', m, ')=', trim(vbt_name)
 964
 965          ! calculate new tendancy
 966          do l=1,dims(3)
 967             do k=1,dims(2)
 968                do i=1,dims(1)
 969                   tend3d(i,k,l)=(scnd3d(i,k,l)-frst3d(i,k,l))/bdyfrq
 970                end do
 971             end do
 972          end do
 973
 974          if (debug) then
 975             write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
 976                  'No.', m, 'Variable: ', trim(vbt_name), &
 977                  'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 978
 979             write(unit=new_unit, fmt='(a, 10i12)') &
 980                  ' new ', (i, i=1,dims(3))
 981
 982             do j=1,dims(1)
 983                write(unit=new_unit, fmt='(i4, 1x, 10e20.7)') &
 984                      j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
 985             end do
 986          end if
 987
 988          ! output new variable at first time level
 989          call da_put_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
 990                                 dims(1), dims(2), dims(3), 1, debug)
 991          call da_put_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
 992                                    dims(1), dims(2), dims(3), 1, debug)
 993
 994          deallocate(frst3d)
 995          deallocate(scnd3d)
 996          deallocate(tend3d)
 997       end do
 998
 999       deallocate(full3d)
1000       if ( var4d_lbc ) deallocate(full3d2)
1001    end do
1002
1003    deallocate(mu)
1004    if ( var4d_lbc ) deallocate(mu2)
1005    deallocate(u)
1006    if ( var4d_lbc ) deallocate(u2)
1007    deallocate(v)
1008    if ( var4d_lbc ) deallocate(v2)
1009
1010 !--------------------- second time level-----------------------------------------
1011 !- for var4d_lbc, we need to update the second time level LBC
1012    if ( var4d_lbc .and. update_lateral_bdy .and. time_level > 1 ) then
1013
1014    east_end=0
1015    north_end=0
1016
1017    cdfid = ncopn(da_file_02, NCNOWRIT, io_status )
1018
1019    ! For 2D variables
1020    ! Get mu, mub, msfu, and msfv
1021
1022
1023    io_status = nf_inq_varid(cdfid, trim(varsf(n)), varid)
1024    if (io_status /= 0 ) then
1025       print '(/"N=",i2," io_status=",i5,5x,"VAR=",a,a)', &
1026                 n, io_status, trim(varsf(n)), " does not exist"
1027       stop
1028    endif
1029
1030    call da_get_dims_cdf( da_file_02, 'MU', dims, &
1031       ndims, debug)
1032
1033    allocate(mu(dims(1), dims(2)))
1034
1035    call da_get_var_2d_real_cdf( da_file_02, &
1036       'MU', mu, dims(1), dims(2), 1, debug)
1037
1038    east_end=dims(1)+1
1039    north_end=dims(2)+1
1040
1041    if (east_end < 1 .or. north_end < 1) then
1042       write(unit=stdout, fmt='(a)') 'Wrong data for Boundary.'
1043       stop
1044    end if
1045
1046    write(unit=stdout,fmt='(/a/)') 'Processing the lateral boundary condition:'
1047
1048    ! boundary variables
1049    bdyname(1)='_BXS'
1050    bdyname(2)='_BXE'
1051    bdyname(3)='_BYS'
1052    bdyname(4)='_BYE'
1053
1054    ! boundary tendancy variables
1055    tenname(1)='_BTXS'
1056    tenname(2)='_BTXE'
1057    tenname(3)='_BTYS'
1058    tenname(4)='_BTYE'
1059
1060    do m=1,4
1061       var_name='MU' // trim(bdyname(m))
1062       vbt_name='MU' // trim(tenname(m))
1063
1064       call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
1065
1066       allocate(frst2d(dims(1), dims(2)))
1067       allocate(scnd2d(dims(1), dims(2)))
1068       allocate(tend2d(dims(1), dims(2)))
1069
1070       ! Get variable at third time level
1071       if (time_level > 2) then
1072          call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), scnd2d, &
1073                                    dims(1), dims(2), 3, debug)
1074       else
1075          call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
1076                                    dims(1), dims(2), 2, debug)
1077          call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
1078                                    dims(1), dims(2), 2, debug)
1079       end if
1080
1081       if (debug) then
1082          write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
1083               'No.', m, 'Variable: ', trim(vbt_name), &
1084               'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
1085
1086          call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
1087                                    dims(1), dims(2), 2, debug)
1088
1089          write(unit=ori_unit, fmt='(a, 10i12)') &
1090               ' old ', (i, i=1,dims(2))
1091          do j=1,dims(1)
1092             write(unit=ori_unit, fmt='(i4, 1x, 10e20.7)') &
1093                   j, (tend2d(j,i), i=1,dims(2))
1094          end do
1095       end if
1096
1097       ! calculate variable at second time level
1098       select case(m)
1099       case (1) ;             ! West boundary
1100          do l=1,dims(2)
1101             do j=1,dims(1)
1102                if (time_level < 3) &
1103                   scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
1104                frst2d(j,l)=mu(l,j)
1105             end do
1106          end do
1107       case (2) ;             ! East boundary
1108          do l=1,dims(2)
1109             do j=1,dims(1)
1110                if (time_level < 3) &
1111                   scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
1112                frst2d(j,l)=mu(east_end-l,j)
1113             end do
1114          end do
1115       case (3) ;             ! South boundary
1116          do l=1,dims(2)
1117             do i=1,dims(1)
1118                if (time_level < 3) &
1119                   scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
1120                frst2d(i,l)=mu(i,l)
1121             end do
1122          end do
1123       case (4) ;             ! North boundary
1124          do l=1,dims(2)
1125             do i=1,dims(1)
1126                if (time_level < 3) &
1127                   scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
1128                frst2d(i,l)=mu(i,north_end-l)
1129             end do
1130          end do
1131       case default ;
1132          write(unit=stdout,fmt=*) 'It is impossible here. mu, m=', m
1133       end select
1134
1135       ! calculate new tendancy
1136       do l=1,dims(2)
1137          do i=1,dims(1)
1138             tend2d(i,l)=(scnd2d(i,l)-frst2d(i,l))/bdyfrq
1139          end do
1140       end do
1141
1142       if (debug) then
1143          write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
1144               'No.', m, 'Variable: ', trim(vbt_name), &
1145               'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
1146
1147          write(unit=new_unit, fmt='(a, 10i12)') &
1148               ' new ', (i, i=1,dims(2))
1149
1150          do j=1,dims(1)
1151             write(unit=new_unit, fmt='(i4, 1x, 10e20.7)') &
1152                   j, (tend2d(j,i), i=1,dims(2))
1153          end do
1154       end if
1155
1156       ! output new variable at first time level
1157       call da_put_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
1158                                 dims(1), dims(2), 2, debug)
1159       ! output new tendancy
1160       call da_put_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
1161                                 dims(1), dims(2), 2, debug)
1162
1163       deallocate(frst2d)
1164       deallocate(scnd2d)
1165       deallocate(tend2d)
1166    end do
1167
1168    !---------------------------------------------------------------------
1169    ! For 3D variables
1170
1171    ! Get U
1172    call da_get_dims_cdf( da_file_02, 'U', dims, ndims, debug)
1173
1174    ! call da_get_att_cdf( da_file_02, 'U', debug)
1175
1176    allocate(u(dims(1), dims(2), dims(3)))
1177
1178    ids=1
1179    ide=dims(1)-1
1180    jds=1
1181    jde=dims(2)
1182    kds=1
1183    kde=dims(3)
1184
1185    call da_get_var_3d_real_cdf( da_file_02, 'U', u, &
1186                              dims(1), dims(2), dims(3), 1, debug)
1187
1188    ! Get V
1189    call da_get_dims_cdf( da_file_02, 'V', dims, ndims, debug)
1190
1191    allocate(v(dims(1), dims(2), dims(3)))
1192
1193    call da_get_var_3d_real_cdf( da_file_02, 'V', v, &
1194                              dims(1), dims(2), dims(3), 1, debug)
1195
1196    if (debug) then
1197       write(unit=stdout, fmt='(a,e20.12,4x)') &
1198            'Before couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
1199            'Before couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
1200    end if
1201
1202    !---------------------------------------------------------------------
1203    ! Couple u, v.
1204    call da_couple_uv ( u, v, mu, mub, msfu, msfv, c1h, c2h, ids, ide, jds, jde, kds, kde)
1205
1206    if (debug) then
1207       write(unit=stdout, fmt='(a,e20.12,4x)') &
1208            'After  couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
1209            'After  couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
1210    end if
1211
1212    !---------------------------------------------------------------------
1213    !For 3D variables
1214
1215    do n=1,num3d
1216       write(unit=stdout, fmt='(a, i3, 2a)') 'Processing: var3d(', n, ')=', trim(var3d(n))
1217
1218       call da_get_dims_cdf( da_file_02, trim(var3d(n)), dims, ndims, debug)
1219
1220       allocate(full3d(dims(1), dims(2), dims(3)))
1221
1222       east_end=dims(1)+1
1223       north_end=dims(2)+1
1224
1225       select case(trim(var3d(n)))
1226       case ('U') ;           ! U
1227          ! var_pref='R' // trim(var3d(n))
1228          var_pref=trim(var3d(n))
1229          full3d(:,:,:)=u(:,:,:)
1230       case ('V') ;           ! V
1231          ! var_pref='R' // trim(var3d(n))
1232          var_pref=trim(var3d(n))
1233          full3d(:,:,:)=v(:,:,:)
1234       case ('W') ;
1235          ! var_pref = 'R' // trim(var3d(n))
1236          var_pref = trim(var3d(n))
1237
1238          call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
1239             full3d, dims(1), dims(2), dims(3), 1, debug)
1240
1241          if (debug) then
1242             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1243                  'Before couple Sample ', trim(var3d(n)), &
1244                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1245          end if
1246
1247          do k=1,dims(3)
1248             do j=1,dims(2)
1249                do i=1,dims(1)
1250                   full3d(i,j,k)=full3d(i,j,k)*(c1f(k)*(mu(i,j)+mub(i,j))+c2f(k))/msfm(i,j)
1251                end do
1252             end do
1253          end do
1254
1255          if (debug) then
1256             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1257                  'After  couple Sample ', trim(var3d(n)), &
1258                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1259          end if
1260       case ('T', 'PH') ;
1261          var_pref=trim(var3d(n))
1262
1263          if ( use_theta_m > 0 .and. trim(var3d(n)) == 'T' ) then
1264            call da_get_var_3d_real_cdf( da_file_02, 'THM', &
1265               full3d, dims(1), dims(2), dims(3), 1, debug)
1266          else
1267            call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
1268               full3d, dims(1), dims(2), dims(3), 1, debug)
1269          end if
1270
1271          if (debug) then
1272             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1273                  'Before couple Sample ', trim(var3d(n)), &
1274                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1275          end if
1276
1277          if ( dims(3) == nlevf ) then
1278             ! applying c1f and c2f for mu
1279             do k=1,dims(3)
1280                do j=1,dims(2)
1281                   do i=1,dims(1)
1282                      full3d(i,j,k)=full3d(i,j,k)*(c1f(k)*(mu(i,j)+mub(i,j))+c2f(k))
1283                   end do
1284                end do
1285             end do
1286          else if ( dims(3) == nlevh ) then
1287             ! applying c1h and c2h for mu
1288             do k=1,dims(3)
1289                do j=1,dims(2)
1290                   do i=1,dims(1)
1291                      full3d(i,j,k)=full3d(i,j,k)*(c1h(k)*(mu(i,j)+mub(i,j))+c2h(k))
1292                   end do
1293                end do
1294             end do
1295          else
1296             write(unit=stdout,fmt=*) 'Error finding matched dimensions.'
1297             stop
1298          end if
1299
1300          if (debug) then
1301             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1302                  'After  couple Sample ', trim(var3d(n)), &
1303                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1304          end if
1305       case ('QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW', 'QGRAUP') ;
1306          ! var_pref='R' // var3d(n)(1:2)
1307          ! var_pref=var3d(n)(1:2)
1308          var_pref=var3d(n)
1309
1310          call da_get_var_3d_real_cdf( da_file_02, trim(var3d(n)), &
1311             full3d, dims(1), dims(2), dims(3), 1, debug)
1312
1313          if (debug) then
1314             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1315                  'Before couple Sample ', trim(var3d(n)), &
1316                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1317          end if
1318
1319          do k=1,dims(3)
1320             do j=1,dims(2)
1321                do i=1,dims(1)
1322                   full3d(i,j,k)=full3d(i,j,k)*(c1h(k)*(mu(i,j)+mub(i,j))+c2h(k))
1323                end do
1324             end do
1325          end do
1326
1327          if (debug) then
1328             write(unit=stdout, fmt='(3a,e20.12,4x)') &
1329                  'After  couple Sample ', trim(var3d(n)), &
1330                  '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
1331          end if
1332       case default ;
1333          write(unit=stdout,fmt=*) 'It is impossible here. var3d(', n, ')=', trim(var3d(n))
1334       end select
1335
1336       do m=1,4
1337          var_name=trim(var_pref) // trim(bdyname(m))
1338          vbt_name=trim(var_pref) // trim(tenname(m))
1339
1340          write(unit=stdout, fmt='(a, i3, 2a)') &
1341             'Processing: bdyname(', m, ')=', trim(var_name)
1342
1343          call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
1344
1345          allocate(frst3d(dims(1), dims(2), dims(3)))
1346          allocate(scnd3d(dims(1), dims(2), dims(3)))
1347          allocate(tend3d(dims(1), dims(2), dims(3)))
1348
1349          ! Get variable at second time level
1350          if (time_level > 2) then
1351             call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), scnd3d, &
1352                                       dims(1), dims(2), dims(3), 3, debug)
1353          else
1354             call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
1355                                       dims(1), dims(2), dims(3), 2, debug)
1356             call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
1357                                       dims(1), dims(2), dims(3), 2, debug)
1358          end if
1359
1360          if (debug) then
1361             write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
1362                  'No.', m, 'Variable: ', trim(vbt_name), &
1363                  'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
1364
1365             call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
1366                                       dims(1), dims(2), dims(3), 2, debug)
1367
1368             write(unit=ori_unit, fmt='(a, 10i12)') &
1369                  ' old ', (i, i=1,dims(3))
1370             do j=1,dims(1)
1371                write(unit=ori_unit, fmt='(i4, 1x, 10e20.7)') &
1372                      j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
1373             end do
1374          end if
1375
1376          select case(trim(bdyname(m)))
1377          case ('_BXS') ;             ! West boundary
1378             do l=1,dims(3)
1379             do k=1,dims(2)
1380             do j=1,dims(1)
1381                if (time_level < 3) &
1382                scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
1383                frst3d(j,k,l)=full3d(l,j,k)
1384             end do
1385             end do
1386             end do
1387          case ('_BXE') ;             ! East boundary
1388             do l=1,dims(3)
1389             do k=1,dims(2)
1390             do j=1,dims(1)
1391                if (time_level < 3) &
1392                scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
1393                frst3d(j,k,l)=full3d(east_end-l,j,k)
1394             end do
1395             end do
1396             end do
1397          case ('_BYS') ;             ! South boundary
1398             do l=1,dims(3)
1399             do k=1,dims(2)
1400             do i=1,dims(1)
1401                if (time_level < 3) &
1402                scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
1403                frst3d(i,k,l)=full3d(i,l,k)
1404             end do
1405             end do
1406             end do
1407          case ('_BYE') ;             ! North boundary
1408             do l=1,dims(3)
1409             do k=1,dims(2)
1410             do i=1,dims(1)
1411                if (time_level < 3) &
1412                scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
1413                frst3d(i,k,l)=full3d(i,north_end-l,k)
1414             end do
1415             end do
1416             end do
1417          case default ;
1418             write(unit=stdout,fmt=*) 'It is impossible here.'
1419             write(unit=stdout,fmt=*) 'bdyname(', m, ')=', trim(bdyname(m))
1420             stop
1421          end select
1422
1423          write(unit=stdout, fmt='(a, i3, 2a)') &
1424             'cal. tend: bdyname(', m, ')=', trim(vbt_name)
1425
1426          ! calculate new tendancy
1427          do l=1,dims(3)
1428             do k=1,dims(2)
1429                do i=1,dims(1)
1430                   tend3d(i,k,l)=(scnd3d(i,k,l)-frst3d(i,k,l))/bdyfrq
1431                end do
1432             end do
1433          end do
1434
1435          if (debug) then
1436             write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
1437                  'No.', m, 'Variable: ', trim(vbt_name), &
1438                  'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
1439
1440             write(unit=new_unit, fmt='(a, 10i12)') &
1441                  ' new ', (i, i=1,dims(3))
1442
1443             do j=1,dims(1)
1444                write(unit=new_unit, fmt='(i4, 1x, 10e20.7)') &
1445                      j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
1446             end do
1447          end if
1448
1449          ! output new variable at first time level
1450          call da_put_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
1451                                 dims(1), dims(2), dims(3), 2, debug)
1452          call da_put_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
1453                                    dims(1), dims(2), dims(3), 2, debug)
1454
1455          deallocate(frst3d)
1456          deallocate(scnd3d)
1457          deallocate(tend3d)
1458       end do
1459
1460       deallocate(full3d)
1461    end do
1462
1463    deallocate(mu)
1464    deallocate(u)
1465    deallocate(v)
1466
1467    end if ! end of update second time level LBC for var4d_lbc
1468
1469    deallocate(mub)
1470    deallocate(msfu)
1471    deallocate(msfv)
1472    deallocate(times)
1473    deallocate(thisbdytime)
1474    deallocate(nextbdytime)
1475
1476  end if ! end if update_lateral_bdy
1477
1478  deallocate ( c1f )
1479  deallocate ( c2f )
1480  deallocate ( c1h )
1481  deallocate ( c2h )
1482
1483  write(unit=stdout,fmt=*) &
1484     '=================================================================='
1485  if ( update_lateral_bdy ) then
1486     write(unit=stdout,fmt=*) 'Lateral boundary tendency updated.'
1487  end if
1488  if ( update_low_bdy ) then
1489     write(unit=stdout,fmt=*) 'Low boundary updated with wrf_input fields.'
1490  end if
1491  if ( update_lsm ) then
1492     write(unit=stdout,fmt=*) 'LSM variables updated with wrf_input fields.'
1493  end if
1494
1495    if (io_status == 0) &
1496       write (unit=stdout,fmt=*) "*** Update_bc completed successfully ***"
1497
1498 end program da_update_bc
1499