var/da/da_dynamics/da_wpec_constraint_lin.inc

   1 subroutine da_wpec_constraint_lin(grid, xbx)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculates TLM of balance equation G(x)
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(domain),   intent(inout) :: grid
  10    type(xbx_type), intent(in)    :: xbx                            ! Header & non-gridded vars.
  11
  12    integer :: i, j, k                 ! Loop counters.
  13    integer :: is, ie                  ! 1st dim. end points.
  14    integer :: js, je                  ! 2nd dim. end points.
  15
  16    real    :: coefx(ims:ime,jms:jme)  ! Multiplicative coefficient.
  17    real    :: coefy(ims:ime,jms:jme)  ! Multiplicative coefficient.
  18    real    :: term_x(ims:ime,jms:jme) ! Balance eqn x term
  19    real    :: term_y(ims:ime,jms:jme) ! Balance eqn y term
  20    real    :: phi_b_x(ims:ime,jms:jme) ! Balance eqn x term
  21    real    :: phi_b_y(ims:ime,jms:jme) ! Balance eqn y term
  22
  23
  24    if (trace_use) call da_trace_entry("da_wpec_constraint_lin")
  25
  26    !---------------------------------------------------------------------------
  27    ! [1.0] Initialise i and set multiplicative constants
  28    !---------------------------------------------------------------------------
  29
  30    is = its; ie = ite; js = jts; je = jte
  31    if (.not.global .and. its == ids) is = ids+1
  32    if (.not.global .and. ite == ide) ie = ide-1
  33    if (jts == jds ) js = jds+1; if (jte == jde) je = jde-1
  34
  35    if (fg_format == fg_format_kma_global) then
  36       coefx = grid%xb%coefx
  37       coefy = grid%xb%coefy
  38    else if( fg_format == fg_format_wrf_arw_regional) then
  39       coefx = grid%xb%coefx
  40       coefy = grid%xb%coefy
  41    else if (fg_format == fg_format_wrf_arw_global) then
  42       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_arw_global  = ',fg_format
  43       call da_error(__FILE__,__LINE__,message(1:1))
  44    else if (fg_format == fg_format_wrf_nmm_regional) then
  45       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_nmm_regional = ',fg_format
  46       call da_error(__FILE__,__LINE__,message(1:1))
  47    else
  48       write(unit=message(1),fmt='(A,I3)') 'Wrong FG_FORMAT = ',fg_format
  49       call da_error(__FILE__,__LINE__,message(1:1))
  50    end if
  51
  52    ! [1.1]
  53
  54    phi_b_x = 0.0
  55    phi_b_y = 0.0
  56
  57
  58    do k = kts, kte
  59
  60       term_x = 0.0
  61       term_y = 0.0
  62
  63       !---------------------------------------------------------------------------
  64       ! [2.0] Calculate RHS of balance equation in gridpoint space:
  65       !---------------------------------------------------------------------------
  66
  67       ! [2.1] Include geostrophic terms in balance eqn if requested:
  68
  69       if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
  70
  71          ! [2.1.1] Calculate term_x = -f rho v~:
  72          term_x = term_x - grid%xb%rho(:,:,k) * grid%xb%cori * grid%xa%v(:,:,k)
  73
  74          ! [2.1.2] Calculate term_y = f rho u~:
  75          term_y = term_y + grid%xb%rho(:,:,k) * grid%xb%cori * grid%xa%u(:,:,k)
  76
  77       end if
  78
  79       ! [2.2] Include cyclostrophic terms in balance eqn if requested:
  80
  81       if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
  82
  83          do j = js, je
  84             do i = is, ie
  85                ! [2.2.1] Calculate term_x = rho M ( u'du/dx + v'du/dy + udu'/dx + vdu'/dy )
  86                term_x(i,j) = term_x(i,j) + grid%xb%rho(i,j,k) * &
  87                    ( coefx(i,j)*grid%xa%u(i,j,k) * ( grid%xb%u(i+1,j,k) - grid%xb%u(i-1,j,k)) + &
  88                      coefy(i,j)*grid%xa%v(i,j,k) * ( grid%xb%u(i,j+1,k) - grid%xb%u(i,j-1,k)) + &
  89                      coefx(i,j)*grid%xb%u(i,j,k) * ( grid%xa%u(i+1,j,k) - grid%xa%u(i-1,j,k)) + &
  90                      coefy(i,j)*grid%xb%v(i,j,k) * ( grid%xa%u(i,j+1,k) - grid%xa%u(i,j-1,k)))
  91
  92                ! [2.2.2] Calculate term_y = rho M ( u'dv/dx + v'dv/dy + udv'/dx + vdv'/dy )
  93                term_y(i,j) = term_y(i,j) + grid%xb%rho(i,j,k) * &
  94                    ( coefx(i,j)*grid%xa%u(i,j,k) * ( grid%xb%v(i+1,j,k) - grid%xb%v(i-1,j,k)) + &
  95                      coefy(i,j)*grid%xa%v(i,j,k) * ( grid%xb%v(i,j+1,k) - grid%xb%v(i,j-1,k)) + &
  96                      coefx(i,j)*grid%xb%u(i,j,k) * ( grid%xa%v(i+1,j,k) - grid%xa%v(i-1,j,k)) + &
  97                      coefy(i,j)*grid%xb%v(i,j,k) * ( grid%xa%v(i,j+1,k) - grid%xa%v(i,j-1,k)))
  98             end do
  99          end do
 100
 101       end if
 102
 103       ! [2.3] Include phi_b terms in balance eqn
 104       do j = js, je
 105          do i = is, ie
 106             phi_b_x(i,j) = coefx(i,j)*(grid%xa%p(i+1,j,k)-grid%xa%p(i-1,j,k) + grid%xb % rho(i,j,k)*(grid%xa%geoh(i+1,j,k)-grid%xa%geoh(i-1,j,k)) )
 107             phi_b_y(i,j) = coefy(i,j)*(grid%xa%p(i,j+1,k)-grid%xa%p(i,j-1,k) + grid%xb % rho(i,j,k)*(grid%xa%geoh(i,j+1,k)-grid%xa%geoh(i,j-1,k)) )
 108          end do
 109       end do
 110
 111    !------------------------------------------------------------------------------
 112    !  [3.0] Solve Grad_p for balance equation :
 113    !------------------------------------------------------------------------------
 114       do j = js, je
 115          do i = is, ie
 116             grid%xa%grad_p_x(i,j,k)=phi_b_x(i,j)+term_x(i,j)
 117             grid%xa%grad_p_y(i,j,k)=phi_b_y(i,j)+term_y(i,j)
 118          end do
 119       end do
 120
 121    end do
 122
 123    if (trace_use) call da_trace_exit("da_wpec_constraint_lin")
 124
 125 end subroutine da_wpec_constraint_lin
 126
 127