var/da/da_dynamics/da_balance_equation_lin.inc

   1 subroutine da_balance_equation_lin(grid, xbx, u, v, phi_b)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculates balanced mass increment phi_b from wind increments.
   5    !
   6    !  Method:  Solve grad**2 Phi_b = - div[ k x rho f u + rho (u.grad ) u ] by
   7    !           1) Calculate RHS of balance equation in gridpt space.
   8    !           2) Solve Del**2 phi_b = RHS for phi_b using double (FCT).
   9    !
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type(domain),   intent(inout) :: grid
  15    type(xbx_type), intent(in)    :: xbx                            ! Header & non-gridded vars.
  16    real,           intent(in)    :: u(ims:ime,jms:jme,kms:kme)     ! u wind comp.
  17    real,           intent(in)    :: v(ims:ime,jms:jme,kms:kme)     ! v wind comp.
  18    real,           intent(out)   :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment.
  19
  20    integer :: i, j, k                 ! Loop counters.
  21    integer :: is, ie                  ! 1st dim. end points.
  22    integer :: js, je                  ! 2nd dim. end points.
  23
  24    real    :: coefx(ims:ime,jms:jme)  ! Multiplicative coefficient.
  25    real    :: coefy(ims:ime,jms:jme)  ! Multiplicative coefficient.
  26    real    :: term_x(ims:ime,jms:jme) ! Balance eqn x term
  27    real    :: term_y(ims:ime,jms:jme) ! Balance eqn y term
  28
  29    real    :: del2phi_b(ims:ime,jms:jme,kms:kme)  ! Del**2 phi_b/M**2
  30
  31    if (trace_use) call da_trace_entry("da_balance_equation_lin")
  32
  33    !---------------------------------------------------------------------------
  34    ! [1.0] Initialise iand set Multipilactive constants
  35    !---------------------------------------------------------------------------
  36
  37    ! Computation to check for edge of domain:
  38
  39    is = its; ie = ite; js = jts; je = jte
  40    if (.not.global .and. its == ids) is = ids+1
  41    if (.not.global .and. ite == ide) ie = ide-1
  42    if (jts == jds ) js = jds+1; 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%coefz
  49       coefy = coefx
  50    else if (fg_format == fg_format_wrf_arw_global) then
  51       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_arw_global  = ',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_wrf_nmm_regional = ',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    ! [1.1] Multiplicative coefficent for conversion RHS->Del**2 phi_b/M**2:
  62
  63    do k = kts, kte
  64       term_x(ims:ime,jms:jme) = 0.0
  65       term_y(ims:ime,jms:jme) = 0.0
  66
  67       !---------------------------------------------------------------------------
  68       ! [2.0] Calculate RHS of balance equation in gridpt space:
  69       !---------------------------------------------------------------------------
  70
  71       ! [2.1] Include geostrophic terms in balance eqn if requested:
  72
  73       if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
  74          call da_balance_geoterm_lin (grid%xb % cori, grid%xb % rho(:,:,k), u(:,:,k), v(:,:,k), &
  75             term_x, term_y)
  76       end if
  77
  78       ! [2.2] Include cyclostrophic terms in balance eqn if requested:
  79
  80       if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
  81          call da_balance_cycloterm_lin (grid%xb % rho(:,:,k), grid%xb % u(:,:,k),   &
  82             grid%xb % v(:,:,k), u(:,:,k), v(:,:,k), grid%xb % coefx(:,:), grid%xb % coefy(:,:), &
  83             term_x(:,:), term_y(:,:))
  84       end if
  85
  86       ! [2.3] Take divergence to get Del**2 phi_b/M**2:
  87
  88       do j = js, je
  89          do i = is, ie
  90           del2phi_b(i,j,k) = -coefx(i,j)*( term_x(i+1,j) - term_x(i-1,j)) &
  91                              -coefy(i,j)*( term_y(i,j+1) - term_y(i,j-1))
  92          end do
  93       end do
  94
  95       ! [2.4] Del**2 Phi_b  boundary conditions:
  96
  97       if (.not. global .and. its == ids ) del2phi_b(its,js:je,k) = 0.0
  98       if (.not. global .and. ite == ide ) del2phi_b(ite,js:je,k) = 0.0
  99       if (jts == jds ) del2phi_b(is:ie,jts,k) = 0.0
 100       if (jte == jde ) del2phi_b(is:ie,jte,k) = 0.0
 101
 102       if (.not. global .and. (its == ids .and. jts == jds) ) del2phi_b(its,jts,k) = 0.0
 103       if (.not. global .and. (its == ids .and. jte == jde) ) del2phi_b(its,jte,k) = 0.0
 104       if (.not. global .and. (ite == ide .and. jts == jds) ) del2phi_b(ite,jts,k) = 0.0
 105       if (.not. global .and. (ite == ide .and. jte == jde) ) del2phi_b(ite,jte,k) = 0.0
 106    end do
 107
 108    !------------------------------------------------------------------------------
 109    !  [3.0] Solve Del**2 phi_b = RHS for phi_b:
 110    !------------------------------------------------------------------------------
 111
 112    call da_solve_poissoneqn_fst(grid, xbx, del2phi_b, phi_b)
 113
 114    if (trace_use) call da_trace_exit("da_balance_equation_lin")
 115
 116 end subroutine da_balance_equation_lin
 117
 118