var/da/da_dynamics/da_balance_cycloterm.inc

   1 subroutine da_balance_cycloterm (xb, k, term_x, term_y)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculates cyclostrophic term in balance equation.
   5    !
   6    !  Method:  Term is rho (u.grad) u on a single level.
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type(xb_type), intent(in)    :: xb           ! First guess structure.
  12    integer,       intent(in)    :: k            ! Model level.
  13    real,          intent(inout) :: term_x(:,:)  ! x component of term.
  14    real,          intent(inout) :: term_y(:,:)  ! y component of term.
  15
  16    integer :: i, j                         ! Loop counters.
  17    integer :: is, ie                       ! 1st dim. end points.
  18    integer :: js, je                       ! 2nd dim. end points.
  19
  20    real    :: data(ims:ime,jms:jme)        ! Temporary storage.
  21
  22    real    :: varb
  23
  24    if (trace_use) call da_trace_entry("da_balance_cycloterm")
  25
  26    !---------------------------------------------------------------------------
  27    ! [1.0] Initialise:
  28    !---------------------------------------------------------------------------
  29
  30    ! Computation to check for edge of domain:
  31    is = its; ie = ite; js = jts; je = jte
  32    if (.not. global .and. its==ids) is = ids+1
  33    if (.not. global .and. ite==ide) ie = ide-1
  34    if (jts==jds) js = jds+1
  35    if (jte==jde) je = jde-1
  36
  37    !---------------------------------------------------------------------------
  38    ! [2.0] Calculate term_x = rho M (u.du/dx + v.du/dy):
  39    !---------------------------------------------------------------------------
  40
  41    ! [2.1] Interior points:
  42
  43    do j = js, je
  44       do i = is, ie
  45          data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%u(i+1,j,k) - &
  46             xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j)*(xb%u(i,j+1,k) - &
  47             xb%u(i,j-1,k))
  48       end do
  49    end do
  50
  51    if (.NOT. global) then ! For global only interior points
  52
  53       ! [2.2] Bottom boundaries:
  54
  55       if (its==ids) then
  56          i = its
  57
  58          do j = js, je
  59             varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i+1,j,k)-xb%u(i+2,j,k)
  60
  61             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
  62                xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
  63          end do
  64       end if
  65
  66       ! [2.3] Top boundaries:
  67
  68       if (ite==ide) then
  69          i = ite
  70
  71          do j = js, je
  72             varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i-1,j,k)+xb%u(i-2,j,k)
  73
  74             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
  75                xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
  76          end do
  77       end if
  78
  79       ! [2.4] Left boundaries:
  80
  81       if (jts==jds) then
  82          j = jts
  83
  84          do i = is, ie
  85             varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i,j+1,k)-xb%u(i,j+2,k)
  86
  87             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
  88                xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
  89          end do
  90       end if
  91
  92       ! [2.5] Right boundaries:
  93
  94       if (jte==jde) then
  95          j = jte
  96
  97          do i = is, ie
  98             varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i,j-1,k)+xb%u(i,j-2,k)
  99
 100             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
 101                xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
 102          end do
 103       end if
 104
 105       ! [2.6] Corner points:
 106
 107       if (its==ids .AND. jts==jds) then
 108          data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
 109       end if
 110
 111       if (ite==ide .AND. jts==jds) then
 112          data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
 113       end if
 114
 115       if (its==ids .AND. jte==jde) then
 116          data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
 117       end if
 118
 119       if (ite==ide .AND. jte==jde) then
 120          data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
 121       end if
 122    end if ! not global
 123
 124    !  [2.7] Multiply by rho  and add to term_x:
 125
 126    term_x(its:ite,jts:jte) = xb%rho(its:ite,jts:jte,k)*data(its:ite,jts:jte) + term_x(its:ite,jts:jte)
 127
 128    !---------------------------------------------------------------------------
 129    ! [3.0] Calculate term_y = rho M (u.dv/dx + v.dv/dy):
 130    !---------------------------------------------------------------------------
 131
 132    ! [3.1] Interior points:
 133
 134    do j = js, je
 135       do i = is, ie
 136          data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
 137                      xb%v(i,j,k) * xb%coefy(i,j)*(xb%v(i,j+1,k) - xb%v(i,j-1,k))
 138       end do
 139    end do
 140
 141
 142    if (.NOT. global) then ! For global only interior points
 143
 144       ! [3.2] Bottom boundaries:
 145
 146       if (its==ids) then
 147          i = its
 148
 149          do j = js, je
 150             varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i+1,j,k)-xb%v(i+2,j,k)
 151
 152             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
 153                         xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
 154          end do
 155       end if
 156
 157       !  [3.3] Top boundaries:
 158
 159       if (ite==ide) then
 160          i = ite
 161
 162          do j = js, je
 163             varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i-1,j,k)+xb%v(i-2,j,k)
 164
 165             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
 166                         xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
 167          end do
 168       end if
 169
 170       ! [3.4] Left boundaries:
 171
 172       if (jts==jds) then
 173          j = jts
 174
 175          do i = is, ie
 176             varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i,j+1,k)-xb%v(i,j+2,k)
 177
 178             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
 179                         xb%v(i,j,k) * xb%coefy(i,j)* varb
 180          end do
 181       end if
 182
 183       ! [3.5] Right boundaries:
 184
 185       if (jte==jde) then
 186          j = jte
 187
 188          do i = is, ie
 189             varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i,j+1,k)+xb%v(i,j+2,k)
 190
 191             data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
 192                         xb%v(i,j,k) * xb%coefy(i,j)* varb
 193          end do
 194       end if
 195
 196       ! [3.6] Corner points:
 197
 198       if (its==ids .AND. jts==jds) then
 199          data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
 200       end if
 201
 202       if (ite==ide .AND. jts==jds) then
 203          data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
 204       end if
 205
 206       if (its==ids .AND. jte==jde) then
 207          data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
 208       end if
 209
 210       if (ite==ide .AND. jte==jde) then
 211          data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
 212       end if
 213    end if ! not global
 214
 215    ! [3.7] Multiply by rho and add to term_y
 216
 217    term_y(its:ite,jts:jte)=xb%rho(its:ite,jts:jte,k)* data(its:ite,jts:jte) + &
 218       term_y(its:ite,jts:jte)
 219
 220    if (trace_use) call da_trace_exit("da_balance_cycloterm")
 221
 222 end subroutine da_balance_cycloterm
 223
 224