var/da/da_dynamics/da_uv_to_vorticity.inc

   1 subroutine da_uv_to_vorticity(xb , u, v, vor)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculate vorticity on a co-ordinate surface, given an input
   5    !           wind field.
   6    !
   7    !  Method:  vor = m**2 (d(v/m)/dx - d(u/m)/dy)
   8    !---------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    type (xb_type), intent(in)    :: xb         ! First guess structure.
  13    real,           intent(in)    :: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
  14    real,           intent(in)    :: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.
  15    real,           intent(inout) :: vor(ims:ime,jms:jme,kms:kme) ! Vorticity.
  16
  17    integer :: i, j, k                      ! Loop counters.
  18    integer :: is, ie                       ! 1st dim. end points.
  19    integer :: js, je                       ! 2nd dim. end points.
  20    real    :: one_third                    ! 1/3.
  21    real    :: inv_2ds                      ! 1/2ds.
  22    real    :: coeff(ims:ime,jms:jme)       ! Mult. coefficient.
  23    real    :: um(ims:ime,jms:jme)          ! Temp. storage of u/m.
  24    real    :: vm(ims:ime,jms:jme)          ! Temp. storage of v/m.
  25
  26    if (trace_use) call da_trace_entry("da_uv_to_vorticity")
  27
  28    !---------------------------------------------------------------------------
  29    ! [1.0] Initialise:
  30    !---------------------------------------------------------------------------
  31
  32    one_third = 1.0 / 3.0
  33
  34    vor = 0.0
  35
  36    !---------------------------------------------------------------------------
  37    ! Computation to check for edge of domain:
  38    !---------------------------------------------------------------------------
  39
  40    is = its
  41    ie = ite
  42    js = jts
  43    je = jte
  44    if (.not. global .and. its == ids) is = ids+1
  45    if (.not. global .and. ite == ide) ie = ide-1
  46    if (jts == jds) js = jds+1
  47    if (jte == jde) je = jde-1
  48
  49    if (.not.global) then
  50       inv_2ds = 0.5 / xb%ds
  51       coeff(its:ite,jts:jte) = xb%map_factor(its:ite,jts:jte) * &
  52          xb%map_factor(its:ite,jts:jte) * inv_2ds
  53    end if
  54
  55    !---------------------------------------------------------------------------
  56    ! [2.0] Calculate vorticity:
  57    !---------------------------------------------------------------------------
  58
  59    do k = kts, kte
  60       ! [2.1] Compute finite difference vorticity at interior points:
  61
  62       if (global) then
  63          do j = js, je
  64             do i = is, ie
  65                vor(i,j,k) = xb%coefy(i,j) * (-u(i,j+1,k) + u(i,j-1,k))  + &
  66                             xb%coefx(i,j) * ( v(i+1,j,k) - v(i-1,j,k))
  67             end do
  68          end do
  69
  70          ! if (global) cycle
  71       else
  72          um(is-1:ie+1,js-1:je+1) = u(is-1:ie+1,js-1:je+1,k) / xb%map_factor(is-1:ie+1,js-1:je+1)
  73          vm(is-1:ie+1,js-1:je+1) = v(is-1:ie+1,js-1:je+1,k) / xb%map_factor(is-1:ie+1,js-1:je+1)
  74
  75          ! [2.1] Compute finite difference vorticity at interior points:
  76
  77          do j = js, je
  78             do i = is, ie
  79                vor(i,j,k) = (-um(i,j+1) + um(i,j-1) + vm(i+1,j) - vm(i-1,j)) * coeff(i,j)
  80             end do
  81          end do
  82
  83          ! [2.2] Impose zero vorticity gradient condition at boundaries:
  84
  85          ! [2.2.1] Bottom boundaries:
  86
  87          if (its == ids) then
  88             i = its
  89             do j = jts, jte
  90                vor(i,j,k) = one_third * (4.0 * vor(i+1,j,k) - vor(i+2,j,k))
  91             end do
  92          end if
  93
  94          ! [2.2.2] Top boundaries:
  95
  96          if (ite == ide) then
  97             i = ite
  98             do j = jts, jte
  99                vor(i,j,k) = one_third * (4.0 * vor(i-1,j,k) - vor(i-2,j,k))
 100             end do
 101          end if
 102
 103          ! [2.2.3] Left boundaries:
 104
 105          if (jts == jds) then
 106             j = jts
 107             do i = its, ite
 108                vor(i,j,k) = one_third * (4.0 * vor(i,j+1,k) - vor(i,j+2,k))
 109             end do
 110          end if
 111
 112          ! [2.2.4] right boundaries:
 113
 114          if (jte == jde) then
 115             j = jte
 116             do i = its, ite
 117                vor(i,j,k) = one_third * (4.0 * vor(i,j-1,k) - vor(i,j-2,k))
 118             end do
 119          end if
 120       end if
 121    end do
 122
 123    if (global) then
 124       call da_set_boundary_3d(vor)
 125    end if
 126
 127    if (trace_use) call da_trace_exit("da_uv_to_vorticity")
 128
 129 end subroutine da_uv_to_vorticity
 130
 131