var/da/da_dynamics/da_wpec_constraint_adj.inc

   1 subroutine da_wpec_constraint_adj(grid, xbx)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Calculates ADM of balance equation G(x)
   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
  13    real :: p(ims:ime,jms:jme,kms:kme) ! pressure increment.
  14    real :: geoh(ims:ime,jms:jme,kms:kme) ! geopotential height increment.
  15    real :: u(ims:ime,jms:jme,kms:kme) ! u wind comp. (dot pts)
  16    real :: v(ims:ime,jms:jme,kms:kme) ! v wind comp. (dot pts)
  17
  18    integer             :: i, j, k          ! Loop counters.
  19    integer             :: is, ie           ! 1st dim. end points.
  20    integer             :: js, je           ! 2nd dim. end points.
  21
  22    real, dimension(ims:ime,jms:jme) :: coefx, &   ! Multiplicative coefficient.
  23                                        coefy, &   ! Multiplicative coefficient.
  24                                        term_x, &  ! Balance eqn x term
  25                                        term_y     ! Balance eqn y term
  26    real    :: phi_b_x(ims:ime,jms:jme) ! Balance eqn x term
  27    real    :: phi_b_y(ims:ime,jms:jme) ! Balance eqn y term
  28
  29    real                 :: data(ims:ime,jms:jme)        ! Work array.
  30    real                 :: var, uar
  31
  32    if (trace_use) call da_trace_entry("da_wpec_constraint_adj")
  33
  34    !---------------------------------------------------------------------------
  35    ! [1.0] Initialise:
  36    !---------------------------------------------------------------------------
  37
  38    is = its; ie = ite; js = jts; je = jte
  39    if (.not. global .and. its == ids) is = ids+1
  40    if (.not. global .and. ite == ide) ie = ide-1
  41    if (jts == jds ) js = jds+1
  42    if (jte == jde ) je = jde-1
  43
  44    if (fg_format == fg_format_kma_global) then
  45       coefx = grid%xb%coefx
  46       coefy = grid%xb%coefy
  47    else if (fg_format == fg_format_wrf_arw_regional) then
  48       coefx = grid%xb%coefx
  49       coefy = grid%xb%coefy
  50    else if (fg_format == fg_format_wrf_arw_global) 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 if (fg_format == fg_format_wrf_nmm_regional) then
  54       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format  = ',fg_format
  55       call da_error(__FILE__,__LINE__,message(1:1))
  56    else
  57       write (unit=message(1),fmt='(A,I3)') ' Wrong FG_FORMAT = ',fg_format
  58       call da_error(__FILE__,__LINE__,message(1:1))
  59    end if
  60
  61    u       = 0.0
  62    v       = 0.0
  63    p       = 0.0
  64    geoh    = 0.0
  65
  66    do k = kts, kte
  67
  68       term_x = 0.0
  69       term_y = 0.0
  70       phi_b_x = 0.0
  71       phi_b_y = 0.0
  72
  73       !---------------------------------------------------------------------------
  74       ! [2.0] Solve Grad_p for balance eqn
  75       !---------------------------------------------------------------------------
  76
  77       phi_b_x = grid%xa%grad_p_x(:,:,k)
  78       phi_b_y = grid%xa%grad_p_y(:,:,k)
  79       term_x = grid%xa%grad_p_x(:,:,k)
  80       term_y = grid%xa%grad_p_y(:,:,k)
  81
  82       !---------------------------------------------------------------------------
  83       ! [3.0] Calculate RHS of balance equation in gridpt space
  84       !---------------------------------------------------------------------------
  85
  86       ! [3.1] Include phi_b terms in balance eqn
  87
  88       do j = je, js, -1
  89          do i = ie, is, -1
  90
  91             p(i+1,j,k) = p(i+1,j,k) + coefx(i,j) * phi_b_x(i,j)
  92             p(i-1,j,k) = p(i-1,j,k) - coefx(i,j) * phi_b_x(i,j)
  93             p(i,j+1,k) = p(i,j+1,k) + coefy(i,j) * phi_b_y(i,j)
  94             p(i,j-1,k) = p(i,j-1,k) - coefy(i,j) * phi_b_y(i,j)
  95
  96             geoh(i+1,j,k) = geoh(i+1,j,k) + coefx(i,j) * phi_b_x(i,j) * grid%xb % rho(i,j,k)
  97             geoh(i-1,j,k) = geoh(i-1,j,k) - coefx(i,j) * phi_b_x(i,j) * grid%xb % rho(i,j,k)
  98             geoh(i,j+1,k) = geoh(i,j+1,k) + coefy(i,j) * phi_b_y(i,j) * grid%xb % rho(i,j,k)
  99             geoh(i,j-1,k) = geoh(i,j-1,k) - coefy(i,j) * phi_b_y(i,j) * grid%xb % rho(i,j,k)
 100
 101          end do
 102       end do
 103
 104       ! [3.2] Include cyclostrophic terms in balance eqn if requested:
 105
 106       if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
 107
 108          ! [3.2.1] Calculate adjoint of term_y = rho M ( u'dv/dx + v'dv/dy + udv'/dx + vdv'/dy ):
 109          data = grid%xb%rho(:,:,k) * term_y
 110
 111          do j = je, js, -1
 112             do i = ie, is, -1
 113                uar = coefx(i,j) * grid%xb%u(i,j,k) * data(i,j)
 114                var = coefy(i,j) * grid%xb%v(i,j,k) * data(i,j)
 115
 116                 u(i,j,k) = u(i,j,k) + coefx(i,j)*data(i,j)*( grid%xb%v(i+1,j,k) - grid%xb%v(i-1,j,k) )
 117                 v(i,j,k) = v(i,j,k) + coefy(i,j)*data(i,j)*( grid%xb%v(i,j+1,k) - grid%xb%v(i,j-1,k) )
 118                 v(i+1,j,k) = v(i+1,j,k) + uar
 119                 v(i-1,j,k) = v(i-1,j,k) - uar
 120                 v(i,j+1,k) = v(i,j+1,k) + var
 121                 v(i,j-1,k) = v(i,j-1,k) - var
 122             end do
 123          end do
 124
 125          ! [3.2.2] Calculate adjoint of term_x = rho M ( u'du/dx + v'du/dy + udu'/dx + vdu'/dy ):
 126          data = grid%xb%rho(:,:,k) * term_x
 127
 128          do j = je, js, -1
 129             do i = ie, is, -1
 130                uar = coefx(i,j) * grid%xb%u(i,j,k) * data(i,j)
 131                var = coefy(i,j) * grid%xb%v(i,j,k) * data(i,j)
 132
 133                 u(i,j,k) = u(i,j,k) + coefx(i,j)*( grid%xb%u(i+1,j,k) - grid%xb%u(i-1,j,k) ) * data(i,j)
 134                 v(i,j,k) = v(i,j,k) + coefy(i,j)*( grid%xb%u(i,j+1,k) - grid%xb%u(i,j-1,k) ) * data(i,j)
 135                 u(i+1,j,k) = u(i+1,j,k) + uar
 136                 u(i-1,j,k) = u(i-1,j,k) - uar
 137                 u(i,j+1,k) = u(i,j+1,k) + var
 138                 u(i,j-1,k) = u(i,j-1,k) - var
 139             end do
 140          end do
 141       end if
 142
 143
 144       ! [3.3] Calculate geostrophic terms in balance eqn:
 145
 146       if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
 147          ! [3.3.1] Calculate term_y = f rho u~:
 148          u(:,:,k) = u(:,:,k) + grid%xb%rho(:,:,k) * grid%xb%cori * term_y
 149
 150          ! [3.3.2] Calculate term_x = -f rho v~:
 151          v(:,:,k) = v(:,:,k) - grid%xb%rho(:,:,k) * grid%xb%cori * term_x
 152
 153       end if
 154
 155    end do
 156
 157    !---------------------------------------------------------------------------
 158    ! [4.0] Store results in grid%xa structure
 159    !---------------------------------------------------------------------------
 160
 161
 162    grid%xa%u=u
 163    grid%xa%v=v
 164    grid%xa%p=p
 165    grid%xa%geoh=geoh
 166
 167
 168    if (trace_use) call da_trace_exit("da_wpec_constraint_adj")
 169
 170 end subroutine da_wpec_constraint_adj
 171
 172