var/da/da_update_bc/da_update_bc_ad.f90

   1 program da_update_bc_ad
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: adjoint version of update BC file from wrfvar output.
   5    ! current version reads only wrf-netcdf file format
   6    !
   7    ! Xin Zhang, 12/12/2011:
   8    !            for 3DVAR use only with FSO
   9    !
  10    !-----------------------------------------------------------------------
  11
  12    use da_netcdf_interface, only : da_get_var_3d_real_cdf, &
  13       da_put_var_3d_real_cdf, da_get_dims_cdf, da_put_var_2d_real_cdf, &
  14       da_get_var_2d_real_cdf, da_get_var_2d_int_cdf, da_get_bdytimestr_cdf, &
  15       da_get_times_cdf, da_get_bdyfrq, stderr, stdout, da_put_var_2d_int_cdf
  16
  17    use da_module_couple_uv_ad, only : da_couple_uv_b
  18
  19    implicit none
  20
  21    include 'netcdf.inc'
  22
  23    integer, parameter :: max_3d_variables = 5, &
  24                          max_2d_variables = 4
  25
  26    character(len=512) :: wrfinput,      &
  27                          init_sens,     &
  28                          wrf_bdy_file,  &
  29                          bdy_sens
  30
  31    character(len=20) :: var_pref, var_name, vbt_name
  32
  33    character(len=20) :: var3d(max_3d_variables), &
  34                         varsf(max_2d_variables)
  35
  36    character(len=10), dimension(4) :: bdyname, tenname
  37
  38    integer           :: ids, ide, jds, jde, kds, kde
  39    integer           :: num3d, num2d, ndims
  40    integer           :: time_level
  41    integer           :: i,j,k,l,m,n
  42
  43    integer, dimension(4) :: dims
  44
  45    real, allocatable, dimension(:,:,:) :: tend3d, scnd3d, frst3d, full3d, full3d9
  46
  47    real, allocatable, dimension(:,:,:) :: u, v
  48    real, allocatable, dimension(:,:,:) :: a_u, a_v
  49
  50    real, allocatable, dimension(:,  :) :: mu, a_mu, mub, msfu, msfv, &
  51                                           tend2d, frst2d
  52
  53    character(len=80), allocatable, dimension(:) :: times, &
  54                                                    thisbdytime, nextbdytime
  55
  56    integer :: east_end, north_end, io_status, cdfid, varid
  57    integer :: iostatus(4)
  58
  59    logical :: debug, update_lateral_bdy
  60
  61    real :: bdyfrq
  62
  63    integer, parameter :: namelist_unit = 7
  64
  65    namelist /control_param/ wrfinput,      &
  66                             init_sens,    &
  67                             wrf_bdy_file, &
  68                             bdy_sens, &
  69                             debug, update_lateral_bdy
  70
  71
  72    wrfinput           = 'wrfinput_d01'
  73    init_sens          = 'init_sens_d01'
  74    wrf_bdy_file       = 'wrfbdy_d01'
  75    bdy_sens           = 'gradient_wrfbdy_d01'
  76
  77    debug              = .false.
  78    update_lateral_bdy = .true.
  79
  80    !---------------------------------------------------------------------
  81    ! Read namelist
  82    !---------------------------------------------------------------------
  83    io_status = 0
  84
  85    open(unit = namelist_unit, file = 'parame.in', &
  86           status = 'old' , access = 'sequential', &
  87           form   = 'formatted', action = 'read', &
  88           iostat = io_status)
  89
  90    if (io_status /= 0) then
  91       write(unit=stdout,fmt=*) 'Error to open namelist file: parame.in.'
  92       write(unit=stdout,fmt=*) 'Will work for updating lateral boundary only.'
  93    else
  94       read(unit=namelist_unit, nml = control_param , iostat = io_status)
  95
  96       if (io_status /= 0) then
  97          write(unit=stdout,fmt=*) 'Error to read control_param. Stopped.'
  98          stop
  99       end if
 100
 101       WRITE(unit=stdout, fmt='(2a)') &
 102            'wrfinput       = ', trim(wrfinput), &
 103            'init_sens     = ', trim(init_sens), &
 104            'wrf_bdy_file  = ', trim(wrf_bdy_file), &
 105            'bdy_sens  = ', trim(bdy_sens)
 106
 107       WRITE(unit=stdout, fmt='(2(a, L10))')             &
 108            'update_lateral_bdy = ', update_lateral_bdy
 109
 110       close(unit=namelist_unit)
 111    end if
 112
 113    ! 3D need update
 114    num3d=5
 115    var3d(1)='A_U'
 116    var3d(2)='A_V'
 117    var3d(3)='A_T'
 118    var3d(4)='A_PH'
 119    var3d(5)='A_QVAPOR'
 120
 121    ! 2D need update
 122    num2d=4
 123    varsf(1)='MUB'
 124    varsf(2)='MU'
 125    varsf(3)='MAPFAC_U'
 126    varsf(4)='MAPFAC_V'
 127
 128    if ( update_lateral_bdy ) then
 129    ! First, the boundary times
 130    call da_get_dims_cdf(wrf_bdy_file, 'Times', dims, ndims, debug)
 131
 132    if (debug) then
 133       write(unit=stdout, fmt='(a,i2,2x,a,4i6)') &
 134            'Times: ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
 135    end if
 136
 137    time_level = dims(2)
 138
 139    if (time_level < 1) then
 140       write(unit=stdout, fmt='(a,i2/a)') &
 141            'time_level = ', time_level, &
 142            'We need at least one time-level BDY.'
 143       stop 'Wrong BDY file.'
 144    end if
 145
 146    allocate(times(dims(2)))
 147    allocate(thisbdytime(dims(2)))
 148    allocate(nextbdytime(dims(2)))
 149
 150    call da_get_times_cdf(wrf_bdy_file, times, dims(2), dims(2), debug)
 151
 152    call da_get_bdytimestr_cdf(wrf_bdy_file, 'thisbdytime', thisbdytime, dims(2), debug)
 153    call da_get_bdytimestr_cdf(wrf_bdy_file, 'nextbdytime', nextbdytime, dims(2), debug)
 154
 155    call da_get_bdyfrq(thisbdytime(1), nextbdytime(1), bdyfrq, debug)
 156
 157    if (debug) then
 158       do n=1, dims(2)
 159          write(unit=stdout, fmt='(3(a, i2, 2a,2x))') &
 160            '       times(', n, ')=', trim(times(n)), &
 161            'thisbdytime (', n, ')=', trim(thisbdytime(n)), &
 162            'nextbdytime (', n, ')=', trim(nextbdytime(n))
 163       end do
 164    end if
 165
 166    end if
 167
 168    east_end=0
 169    north_end=0
 170
 171    cdfid = ncopn(wrfinput, NCNOWRIT, io_status )
 172
 173 !---------------------------------------------------------------
 174    ! For 2D variables
 175    ! Get mu, mub, msfu, and msfv
 176
 177    do n=1,num2d
 178
 179       io_status = nf_inq_varid(cdfid, trim(varsf(n)), varid)
 180       if (io_status /= 0 ) then
 181          print '(/"N=",i2," io_status=",i5,5x,"VAR=",a,a)', &
 182                    n, io_status, trim(varsf(n)), " does not exist"
 183          cycle
 184       endif
 185
 186       call da_get_dims_cdf( wrfinput, trim(varsf(n)), dims, &
 187          ndims, debug)
 188
 189       select case(trim(varsf(n)))
 190       case ('MU') ;
 191          if ( .not. update_lateral_bdy ) cycle
 192
 193          allocate(mu(dims(1), dims(2)))
 194          allocate(a_mu(dims(1), dims(2)))
 195          a_mu=0.0
 196
 197          call da_get_var_2d_real_cdf( wrfinput, &
 198             trim(varsf(n)), mu, dims(1), dims(2), 1, debug)
 199
 200       case ('MUB') ;
 201          if ( .not. update_lateral_bdy ) cycle
 202
 203          allocate(mub(dims(1), dims(2)))
 204
 205          call da_get_var_2d_real_cdf( wrfinput, trim(varsf(n)), mub, &
 206                                    dims(1), dims(2), 1, debug)
 207       case ('MAPFAC_U') ;
 208          if ( .not. update_lateral_bdy ) cycle
 209
 210          allocate(msfu(dims(1), dims(2)))
 211
 212          call da_get_var_2d_real_cdf( wrfinput, trim(varsf(n)), msfu, &
 213                                    dims(1), dims(2), 1, debug)
 214       case ('MAPFAC_V') ;
 215          if ( .not. update_lateral_bdy ) cycle
 216
 217          allocate(msfv(dims(1), dims(2)))
 218
 219          call da_get_var_2d_real_cdf( wrfinput, trim(varsf(n)), msfv, &
 220                                    dims(1), dims(2), 1, debug)
 221
 222       case default ;
 223             write(unit=stdout,fmt=*) 'It is impossible here. varsf(n)=', trim(varsf(n))
 224       end select
 225    end do
 226
 227
 228    ! Get U
 229    call da_get_dims_cdf( wrfinput, 'U', dims, ndims, debug)
 230
 231    allocate(u(dims(1), dims(2), dims(3)))
 232    allocate(a_u(dims(1), dims(2), dims(3)))
 233    a_u=0.0
 234
 235    ids=1
 236    ide=dims(1)-1
 237    jds=1
 238    jde=dims(2)
 239    kds=1
 240    kde=dims(3)
 241
 242    call da_get_var_3d_real_cdf( wrfinput, 'U', u, &
 243                              dims(1), dims(2), dims(3), 1, debug)
 244
 245    ! Get V
 246    call da_get_dims_cdf( wrfinput, 'V', dims, ndims, debug)
 247
 248    allocate(v(dims(1), dims(2), dims(3)))
 249    allocate(a_v(dims(1), dims(2), dims(3)))
 250    a_v=0.0
 251
 252    call da_get_var_3d_real_cdf( wrfinput, 'V', v, &
 253                              dims(1), dims(2), dims(3), 1, debug)
 254
 255    ! boundary variables
 256    bdyname(1)='_BXS'
 257    bdyname(2)='_BXE'
 258    bdyname(3)='_BYS'
 259    bdyname(4)='_BYE'
 260
 261    ! boundary tendancy variables
 262    tenname(1)='_BTXS'
 263    tenname(2)='_BTXE'
 264    tenname(3)='_BTYS'
 265    tenname(4)='_BTYE'
 266
 267    !For 3D variables
 268
 269    do n=1,num3d
 270       write(unit=stdout, fmt='(a, i3, 2a)') 'Processing: var3d(', n, ')=', trim(var3d(n))
 271
 272       call da_get_dims_cdf( init_sens, trim(var3d(n)), dims, ndims, debug)
 273
 274       allocate(full3d(dims(1), dims(2), dims(3)))
 275       full3d=0.0
 276
 277       allocate(full3d9(dims(1), dims(2), dims(3)))
 278
 279       east_end=dims(1)+1
 280       north_end=dims(2)+1
 281
 282       var_pref=trim(var3d(n))
 283
 284       do m=1,4
 285          var_name=trim(var_pref) // trim(bdyname(m))
 286          vbt_name=trim(var_pref) // trim(tenname(m))
 287
 288          write(unit=stdout, fmt='(a, i3, 2a)') &
 289             'Processing: bdyname(', m, ')=', trim(var_name)
 290
 291          call da_get_dims_cdf( bdy_sens, trim(var_name), dims, ndims, debug)
 292
 293          allocate(frst3d(dims(1), dims(2), dims(3)))
 294          allocate(tend3d(dims(1), dims(2), dims(3)))
 295
 296          write(unit=stdout, fmt='(a, i3, 2a)') &
 297             'adjoint cal. tend: bdyname(', m, ')=', trim(vbt_name)
 298
 299          ! output new variable at first time level
 300          call da_get_var_3d_real_cdf( bdy_sens, trim(var_name), frst3d, &
 301                                    dims(1), dims(2), dims(3), 1, debug)
 302          call da_get_var_3d_real_cdf( bdy_sens, trim(vbt_name), tend3d, &
 303                                    dims(1), dims(2), dims(3), 1, debug)
 304
 305          ! calculate new tendancy
 306          do l=1,dims(3)
 307             do k=1,dims(2)
 308                do i=1,dims(1)
 309                   frst3d(i,k,l)=frst3d(i,k,l)-tend3d(i,k,l)/bdyfrq
 310                   tend3d(i,k,l)=0.0
 311                end do
 312             end do
 313          end do
 314
 315          select case(trim(bdyname(m)))
 316          case ('_BXS') ;             ! West boundary
 317             do l=1,dims(3)
 318             do k=1,dims(2)
 319             do j=1,dims(1)
 320                full3d(l,j,k)=full3d(l,j,k)+frst3d(j,k,l)
 321                frst3d(j,k,l)=0.0
 322             end do
 323             end do
 324             end do
 325          case ('_BXE') ;             ! East boundary
 326             do l=1,dims(3)
 327             do k=1,dims(2)
 328             do j=1,dims(1)
 329                full3d(east_end-l,j,k)=full3d(east_end-l,j,k)+frst3d(j,k,l)
 330                frst3d(j,k,l)=0.0
 331             end do
 332             end do
 333             end do
 334          case ('_BYS') ;             ! South boundary
 335             do l=1,dims(3)
 336             do k=1,dims(2)
 337             do i=1,dims(1)
 338                full3d(i,l,k)=full3d(i,l,k)+frst3d(i,k,l)
 339                frst3d(i,k,l)=0.0
 340             end do
 341             end do
 342             end do
 343          case ('_BYE') ;             ! North boundary
 344             do l=1,dims(3)
 345             do k=1,dims(2)
 346             do i=1,dims(1)
 347                full3d(i,north_end-l,k)=full3d(i,north_end-l,k)+frst3d(i,k,l)
 348                frst3d(i,k,l)=0.0
 349             end do
 350             end do
 351             end do
 352          case default ;
 353             write(unit=stdout,fmt=*) 'It is impossible here.'
 354             write(unit=stdout,fmt=*) 'bdyname(', m, ')=', trim(bdyname(m))
 355             stop
 356          end select
 357
 358          deallocate(frst3d)
 359          deallocate(tend3d)
 360       end do
 361
 362       select case(trim(var3d(n)))
 363       case ('A_U') ;           ! U
 364          var_pref=trim(var3d(n))
 365          a_u(:,:,:)=a_u(:,:,:)+full3d(:,:,:)
 366          full3d(:,:,:)=0.0
 367       case ('A_V') ;           ! V
 368          var_pref=trim(var3d(n))
 369          a_v(:,:,:)=a_v(:,:,:)+full3d(:,:,:)
 370          full3d(:,:,:)=0.0
 371       case ('A_T') ;
 372          var_pref=trim(var3d(n))
 373
 374          call da_get_dims_cdf( init_sens, trim(var3d(n)), dims, ndims, debug)
 375          call da_get_var_3d_real_cdf( wrfinput, 'T', &
 376             full3d9, dims(1), dims(2), dims(3), 1, debug)
 377
 378          do k=1,dims(3)
 379             do j=1,dims(2)
 380                do i=1,dims(1)
 381                   a_mu(i,j)=a_mu(i,j)+full3d(i,j,k)*full3d9(i,j,k)
 382                   full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))
 383                end do
 384             end do
 385          end do
 386
 387          call da_get_var_3d_real_cdf( init_sens, 'A_T', &
 388             full3d9, dims(1), dims(2), dims(3), 1, debug)
 389
 390          full3d=full3d+full3d9
 391          full3d9=0.0
 392
 393          call da_put_var_3d_real_cdf( init_sens, trim(var3d(n)), &
 394             full3d, dims(1), dims(2), dims(3), 1, debug)
 395
 396          full3d=0.0
 397
 398       case ('A_PH') ;
 399          var_pref=trim(var3d(n))
 400
 401          call da_get_dims_cdf( init_sens, trim(var3d(n)), dims, ndims, debug)
 402          call da_get_var_3d_real_cdf( wrfinput, 'PH', &
 403             full3d9, dims(1), dims(2), dims(3), 1, debug)
 404
 405          do k=1,dims(3)
 406             do j=1,dims(2)
 407                do i=1,dims(1)
 408                   a_mu(i,j)=a_mu(i,j)+full3d(i,j,k)*full3d9(i,j,k)
 409                   full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))
 410                end do
 411             end do
 412          end do
 413
 414          call da_get_var_3d_real_cdf( init_sens, 'A_PH', &
 415             full3d9, dims(1), dims(2), dims(3), 1, debug)
 416
 417          full3d=full3d+full3d9
 418          full3d9=0.0
 419
 420          call da_put_var_3d_real_cdf( init_sens, trim(var3d(n)), &
 421             full3d, dims(1), dims(2), dims(3), 1, debug)
 422
 423          full3d=0.0
 424
 425       case ('A_QVAPOR') ;
 426          var_pref=var3d(n)
 427
 428          call da_get_dims_cdf( init_sens, trim(var3d(n)), dims, ndims, debug)
 429          call da_get_var_3d_real_cdf( wrfinput, 'QVAPOR', &
 430             full3d9, dims(1), dims(2), dims(3), 1, debug)
 431
 432          do k=1,dims(3)
 433             do j=1,dims(2)
 434                do i=1,dims(1)
 435                   a_mu(i,j)=a_mu(i,j)+full3d(i,j,k)*full3d9(i,j,k)
 436                   full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))
 437                end do
 438             end do
 439          end do
 440
 441          call da_get_var_3d_real_cdf( init_sens, 'A_QVAPOR', &
 442             full3d9, dims(1), dims(2), dims(3), 1, debug)
 443
 444          full3d=full3d+full3d9
 445          full3d9=0.0
 446
 447          call da_put_var_3d_real_cdf( init_sens, trim(var3d(n)), &
 448             full3d, dims(1), dims(2), dims(3), 1, debug)
 449
 450          full3d=0.0
 451
 452       case default ;
 453          write(unit=stdout,fmt=*) 'It is impossible here. var3d(', n, ')=', trim(var3d(n))
 454       end select
 455
 456       deallocate(full3d)
 457       deallocate(full3d9)
 458
 459    end do
 460
 461    !---------------------------------------------------------------------
 462    ! Couple u, v.
 463    call da_couple_uv_b ( u, a_u, v, a_v, mu, a_mu, mub, msfu, msfv, ids, ide, jds, jde, kds, kde)
 464
 465    call da_get_dims_cdf( init_sens, 'A_V', dims, ndims, debug)
 466    call da_get_var_3d_real_cdf( init_sens, 'A_V', v, &
 467                              dims(1), dims(2), dims(3), 1, debug)
 468    a_v=a_v+v
 469    v=0.0
 470    call da_put_var_3d_real_cdf( init_sens, 'A_V', a_v, &
 471                              dims(1), dims(2), dims(3), 1, debug)
 472    a_v=0.0
 473
 474    call da_get_dims_cdf( init_sens, 'A_U', dims, ndims, debug)
 475    call da_get_var_3d_real_cdf( init_sens, 'A_U', u, &
 476                              dims(1), dims(2), dims(3), 1, debug)
 477    a_u=a_u+u
 478    u=0.0
 479    call da_put_var_3d_real_cdf( init_sens, 'A_U', a_u, &
 480                              dims(1), dims(2), dims(3), 1, debug)
 481    a_u=0.0
 482
 483    call da_get_dims_cdf( wrfinput, 'MU', dims, &
 484       ndims, debug)
 485
 486    east_end=dims(1)+1
 487    north_end=dims(2)+1
 488
 489    do m=1,4
 490       var_name='A_MU' // trim(bdyname(m))
 491       vbt_name='A_MU' // trim(tenname(m))
 492
 493       call da_get_dims_cdf( bdy_sens, trim(var_name), dims, ndims, debug)
 494
 495       allocate(frst2d(dims(1), dims(2)))
 496       allocate(tend2d(dims(1), dims(2)))
 497
 498       ! get new variable at first time level
 499       call da_get_var_2d_real_cdf( bdy_sens, trim(var_name), frst2d, &
 500                                 dims(1), dims(2), 1, debug)
 501       ! output new tendancy
 502       call da_get_var_2d_real_cdf( bdy_sens, trim(vbt_name), tend2d, &
 503                                 dims(1), dims(2), 1, debug)
 504
 505       ! calculate new tendancy
 506       do l=1,dims(2)
 507          do i=1,dims(1)
 508             frst2d(i,l)=frst2d(i,l)-tend2d(i,l)/bdyfrq
 509             tend2d(i,l)=0.0
 510          end do
 511       end do
 512
 513       ! calculate variable at first time level
 514       select case(m)
 515       case (1) ;             ! West boundary
 516          do l=1,dims(2)
 517             do j=1,dims(1)
 518                a_mu(l,j)=a_mu(l,j)+frst2d(j,l)
 519                frst2d(j,l)=0.0
 520             end do
 521          end do
 522       case (2) ;             ! East boundary
 523          do l=1,dims(2)
 524             do j=1,dims(1)
 525                a_mu(east_end-l,j)=a_mu(east_end-l,j)+frst2d(j,l)
 526                frst2d(j,l)=0.0
 527             end do
 528          end do
 529       case (3) ;             ! South boundary
 530          do l=1,dims(2)
 531             do i=1,dims(1)
 532                a_mu(i,l)=a_mu(i,l)+frst2d(i,l)
 533                frst2d(i,l)=0.0
 534             end do
 535          end do
 536       case (4) ;             ! North boundary
 537          do l=1,dims(2)
 538             do i=1,dims(1)
 539                a_mu(i,north_end-l)=a_mu(i,north_end-l)+frst2d(i,l)
 540                frst2d(i,l)=0.0
 541             end do
 542          end do
 543       case default ;
 544          write(unit=stdout,fmt=*) 'It is impossible here. mu, m=', m
 545       end select
 546
 547       deallocate(frst2d)
 548       deallocate(tend2d)
 549    end do
 550
 551    call da_get_dims_cdf( init_sens, 'A_MU', dims, &
 552       ndims, debug)
 553    call da_get_var_2d_real_cdf( init_sens, &
 554       'A_MU', mu, dims(1), dims(2), 1, debug)
 555    a_mu=a_mu+mu
 556    mu=0.0
 557    call da_put_var_2d_real_cdf( init_sens, &
 558       'A_MU', a_mu, dims(1), dims(2), 1, debug)
 559    a_mu=0.0
 560
 561    !---------------------------------------------------------------------
 562
 563
 564    deallocate(mu)
 565    deallocate(u)
 566    deallocate(v)
 567    deallocate(a_mu)
 568    deallocate(a_u)
 569    deallocate(a_v)
 570    deallocate(mub)
 571    deallocate(msfu)
 572    deallocate(msfv)
 573
 574    deallocate(times)
 575    deallocate(thisbdytime)
 576    deallocate(nextbdytime)
 577
 578    write (unit=stdout,fmt=*) "*** Adjoint of update_bc completed successfully ***"
 579
 580 end program da_update_bc_ad
 581