var/da/da_dynamics/da_divergence_constraint_adj.inc

   1 subroutine da_divergence_constraint_adj(grid, u, v, div)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Adjoint of the subroutine da_divergence_constraint
   5    !                        d   U      d   V
   6    !           Div = m^2 *[---(---) + ---(---) ]
   7    !                        dx  m      dy  M
   8    !---------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    type (domain), intent(inout) :: grid
  13
  14    real, intent(out)  :: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
  15    real, intent(out)  :: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.
  16    real, intent(inout):: div(ims:ime,jms:jme,kms:kme) ! Divergence.
  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    real               :: one_third                    ! 1/3.
  22
  23    real               :: coeff, inv_2ds
  24    real               :: um(ims:ime,jms:jme)          ! Temp. storage of mu*u/m.
  25    real               :: vm(ims:ime,jms:jme)          ! Temp. storage of mu*v/m
  26    real               :: mu(ims:ime,jms:jme)          ! Temp. storage of mu
  27
  28    if (trace_use) call da_trace_entry("da_divergence_constraint_adj")
  29
  30    !---------------------------------------------------------------------------
  31    ! [1.0] Initialise:
  32    !---------------------------------------------------------------------------
  33
  34    one_third = 1.0 / 3.0
  35
  36    is = its - 1; ie = ite + 1; js = jts - 1; je = jte + 1
  37
  38    if (.not. global .and. its == ids) is = ids+1
  39    if (.not. global .and. ite == ide) ie = ide-1
  40    if (jts == jds) js = jds+1; if (jte == jde) je = jde-1
  41
  42    mu = grid%xb%psfc - grid%xb%ptop
  43
  44    !---------------------------------------------------------------------------
  45    ! [2.0] Calculate divergence:
  46    !---------------------------------------------------------------------------
  47
  48    if (.not. global) then
  49       inv_2ds = 0.5 / grid%xb%ds
  50       do k =kds, kde
  51          um(ims:ime,jms:jme) = 0.0
  52          vm(ims:ime,jms:jme) = 0.0
  53          ! [2.2] Impose zero divergence gradient condition at boundaries:
  54
  55          ! [2.2.4] Right boundaries:
  56
  57          if (jte == jde) then
  58             j = jte
  59             do i = its, ite    ! This is different to original
  60                div(i,j-1,k)=div(i,j-1,k)+div(i,j,k)*one_third*4.0
  61                div(i,j-2,k)=div(i,j-2,k)-div(i,j,k)*one_third
  62                div(i,j,k)=0.0
  63             end do
  64          end if
  65
  66          ! [2.2.3] Left boundaries:
  67
  68          if (jts == jds) then
  69             j = jts
  70             do i = its, ite    ! This is different to original
  71                div(i,j+1,k)=div(i,j+1,k)+div(i,j,k)*one_third*4.0
  72                div(i,j+2,k)=div(i,j+2,k)-div(i,j,k)*one_third
  73                div(i,j,k)=0.0
  74             end do
  75          end if
  76
  77          ! [2.2.2] Top boundaries:
  78
  79          if (ite == ide) then
  80             i = ite
  81             do j = jts, jte
  82                div(i-1,j,k)=div(i-1,j,k)+div(i,j,k)*one_third*4.0
  83                div(i-2,j,k)=div(i-2,j,k)-div(i,j,k)*one_third
  84                div(i,j,k)=0.0
  85             end do
  86          end if
  87
  88          ! [2.2.1] Bottom boundaries:
  89
  90          if (its == ids) then
  91             i = its
  92             do j = jts, jte
  93                div(i+1,j,k)=div(i+1,j,k)+div(i,j,k)*one_third*4.0
  94                div(i+2,j,k)=div(i+2,j,k)-div(i,j,k)*one_third
  95                div(i,j,k)=0.0
  96             end do
  97          end if
  98
  99          ! [2.1] Compute fd divergence at interior points:
 100          ! Computation to check for edge of domain:
 101          ! This is only for adjoint, as we have to cross the processor boundary
 102          ! to get the contribution.
 103
 104          grid%xp%vxy(its:ite, jts:jte) = div(its:ite, jts:jte, k)
 105 #ifdef DM_PARALLEL
 106 #include "HALO_2D_WORK.inc"
 107 #endif
 108
 109          div(is, js:je, k) = grid%xp%vxy(is, js:je)
 110          div(ie, js:je, k) = grid%xp%vxy(ie, js:je)
 111          div(is:ie, js, k) = grid%xp%vxy(is:ie, js)
 112          div(is:ie, je, k) = grid%xp%vxy(is:ie, je)
 113
 114          do j = js, je
 115             do i = is, ie
 116                coeff = grid%xb%map_factor(i,j) * grid%xb%map_factor(i,j) * inv_2ds
 117                um(i+1,j)=um(i+1,j)+div(i,j,k)*coeff
 118                um(i-1,j)=um(i-1,j)-div(i,j,k)*coeff
 119                vm(i,j+1)=vm(i,j+1)+div(i,j,k)*coeff
 120                vm(i,j-1)=vm(i,j-1)-div(i,j,k)*coeff
 121             end do
 122          end do
 123
 124          u(is-1:ie+1,js-1:je+1,k) = u(is-1:ie+1,js-1:je+1,k) + &
 125                                     um(is-1:ie+1,js-1:je+1) * &
 126                                     mu(is-1:ie+1,js-1:je+1) / &
 127                                     grid%xb%map_factor(is-1:ie+1,js-1:je+1)
 128          v(is-1:ie+1,js-1:je+1,k) = v(is-1:ie+1,js-1:je+1,k) + &
 129                                     vm(is-1:ie+1,js-1:je+1) * &
 130                                     mu(is-1:ie+1,js-1:je+1) / &
 131                                     grid%xb%map_factor(is-1:ie+1,js-1:je+1)
 132       end do
 133
 134    else ! global
 135       call da_set_boundary_3d(div)
 136
 137       do k =kds, kde
 138          !-------------------------------------------------------------------------
 139          ! [2.1] Compute fd divergence at interior points:
 140          !-------------------------------------------------------------------------
 141
 142          do j = je, js, -1
 143             do i = ie, is, -1
 144                u(i+1,j,k) = u(i+1,j,k) + grid%xb%coefx(i,j) * div(i,j,k)
 145                u(i-1,j,k) = u(i-1,j,k) - grid%xb%coefx(i,j) * div(i,j,k)
 146                v(i,j+1,k) = v(i,j+1,k) + grid%xb%coefy(i,j) * div(i,j,k)
 147                v(i,j-1,k) = v(i,j-1,k) - grid%xb%coefy(i,j) * div(i,j,k)
 148             end do
 149          end do
 150       end do
 151    end if
 152
 153    div = 0.0
 154
 155    if (trace_use) call da_trace_exit("da_divergence_constraint_adj")
 156
 157 end subroutine da_divergence_constraint_adj
 158
 159