var/da/da_dynamics/da_balance_equation_adj.inc

   1 subroutine da_balance_equation_adj(grid, xbx, phi_b, u, v)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Adjoint of da_balance_equation
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(domain), intent(inout)               :: grid
  10
  11    type (xbx_type),intent(in) :: xbx              ! Header & non-gridded vars.
  12    real, intent(in)    :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment.
  13    real, intent(inout) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp. (dot pts)
  14    real, intent(inout) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp. (dot pts)
  15
  16    integer             :: i, j, k          ! Loop counters.
  17    integer             :: is, ie           ! 1st dim. end points.
  18    integer             :: js, je           ! 2nd dim. end points.
  19
  20    real, dimension(ims:ime,jms:jme) :: coefx, &   ! Multiplicative coefficient.
  21                                        coefy, &   ! Multiplicative coefficient.
  22                                        term_x, &  ! Balance eqn x term
  23                                        term_y     ! Balance eqn y term
  24
  25    real, dimension(ims:ime,jms:jme,kms:kme) :: del2phi_b  ! Del**2 phi_b/M**2
  26    real                       :: coeff1, coeff2   ! Multiplicative coefficient.
  27
  28    if (trace_use) call da_trace_entry("da_balance_equation_adj")
  29
  30    !---------------------------------------------------------------------------
  31    ! [1.0] Initialise:
  32    !---------------------------------------------------------------------------
  33
  34    ! Computation to check for edge of domain:
  35    is = its-1; ie = ite+1; js = jts-1; je = jte+1
  36    if (.not. global .and. its == ids) is = ids+1
  37    if (.not. global .and. ite == ide) ie = ide-1
  38    if (jts == jds ) js = jds+1
  39    if (jte == jde ) je = jde-1
  40
  41    if (fg_format == fg_format_kma_global) then
  42       coefx(is:ie,js:je) = grid%xb%coefx(is:ie,js:je)
  43       coefy(is:ie,js:je) = grid%xb%coefy(is:ie,js:je)
  44    else if (fg_format == fg_format_wrf_arw_global) then
  45       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format  = ',fg_format
  46       call da_error(__FILE__,__LINE__,message(1:1))
  47    else if (fg_format == fg_format_wrf_arw_regional) then
  48       coefx(is:ie,js:je) = grid%xb%coefz(is:ie,js:je)
  49       coefy(is:ie,js:je) = coefx(is:ie,js:je)
  50    else if (fg_format == fg_format_wrf_nmm_regional) then
  51       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format  = ',fg_format
  52       call da_error(__FILE__,__LINE__,message(1:1))
  53    else
  54       write (unit=message(1),fmt='(A,I3)') ' Wrong FG_FORMAT = ',fg_format
  55       call da_error(__FILE__,__LINE__,message(1:1))
  56    end if
  57
  58    ! [1.1] Multiplicative coefficent for conversion RHS->Del**2 phi_b/M**2:
  59
  60    del2phi_b(:,:,:) = 0.0
  61
  62    !---------------------------------------------------------------------------
  63    ! [3.0] Solve Del**2 phi_b = RHS for phi_b:
  64    !---------------------------------------------------------------------------
  65
  66    call da_solve_poissoneqn_fst_adj(grid,xbx, phi_b, del2phi_b)
  67
  68    !---------------------------------------------------------------------------
  69    ! [2.0] Calculate RHS of balance equation in gridpt space:
  70    !---------------------------------------------------------------------------
  71
  72    do k = kts, kte
  73       ! [2.4] Del**2 Phi_b boundary conditions (null as zero boundary conditions):
  74
  75       ! [2.3] Take divergence to get Del**2 phi_b/M**2:
  76
  77       term_x(ims:ime,jms:jme) = 0.0
  78       term_y(ims:ime,jms:jme) = 0.0
  79
  80       do j = je, js, -1
  81          do i = ie, is, -1
  82             coeff1 = coefx(i,j) * del2phi_b(i,j,k)
  83             coeff2 = coefy(i,j) * del2phi_b(i,j,k)
  84
  85             term_x(i+1,j) = term_x(i+1,j) - coeff1
  86             term_x(i-1,j) = term_x(i-1,j) + coeff1
  87             term_y(i,j+1) = term_y(i,j+1) - coeff2
  88             term_y(i,j-1) = term_y(i,j-1) + coeff2
  89          end do
  90       end do
  91
  92       ! [2.2] Include cyclostrophic terms in balance eqn if requested:
  93
  94       if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
  95          call da_balance_cycloterm_adj (grid%xb%rho(:,:,k),grid%xb%u(:,:,k),&
  96             grid%xb%v(:,:,k), u(:,:,k), v(:,:,k), grid%xb%coefx(:,:), grid%xb%coefy(:,:),&
  97             term_x(:,:), term_y(:,:))
  98       end if
  99
 100
 101       ! [2.1] Calculate geostrophic terms in balance eqn:
 102
 103       if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
 104          call da_balance_geoterm_adj (grid%xb%cori, grid%xb%rho(:,:,k), term_x, term_y, &
 105             u(:,:,k), v(:,:,k))
 106       end if
 107    end do
 108
 109    if (trace_use) call da_trace_exit("da_balance_equation_adj")
 110
 111 end subroutine da_balance_equation_adj
 112
 113