var/da/da_dynamics/da_wpec_constraint.inc

   1 subroutine da_wpec_constraint(grid, xbx)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculates 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    real   :: u(ims:ime,jms:jme,kms:kme)     ! u wind comp.
  13    real   :: v(ims:ime,jms:jme,kms:kme)     ! v wind comp.
  14    real   :: phi_b(ims:ime,jms:jme,kms:kme)
  15
  16    integer :: i, j, k                 ! Loop counters.
  17    integer :: is, ie                  ! 1st dim. end points.
  18    integer :: js, je                  ! 2nd dim. end points.
  19
  20    real    :: coefx(ims:ime,jms:jme)  ! Multiplicative coefficient.
  21    real    :: coefy(ims:ime,jms:jme)  ! Multiplicative coefficient.
  22    real    :: term_x(ims:ime,jms:jme) ! Balance eqn x term
  23    real    :: term_y(ims:ime,jms:jme) ! Balance eqn y term
  24    real    :: phi_b_x(ims:ime,jms:jme) ! Balance eqn x term
  25    real    :: phi_b_y(ims:ime,jms:jme) ! Balance eqn y term
  26
  27
  28
  29    if (trace_use) call da_trace_entry("da_wpec_constraint")
  30
  31    !---------------------------------------------------------------------------
  32    ! [1.0] Initialise i and set multiplicative constants
  33    !---------------------------------------------------------------------------
  34
  35    ! Computation to check for edge of domain:
  36
  37    is = its; ie = ite; js = jts; je = jte
  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
  43    if (fg_format == fg_format_kma_global) then
  44       coefx = grid%xb%coefx
  45       coefy = grid%xb%coefy
  46    else if( fg_format == fg_format_wrf_arw_regional) then
  47       coefx = grid%xb%coefx
  48       coefy = grid%xb%coefy
  49    else if (fg_format == fg_format_wrf_arw_global) then
  50       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_arw_global  = ',fg_format
  51       call da_error(__FILE__,__LINE__,message(1:1))
  52    else if (fg_format == fg_format_wrf_nmm_regional) then
  53       write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_nmm_regional = ',fg_format
  54       call da_error(__FILE__,__LINE__,message(1:1))
  55    else
  56       write(unit=message(1),fmt='(A,I3)') 'Wrong FG_FORMAT = ',fg_format
  57       call da_error(__FILE__,__LINE__,message(1:1))
  58    end if
  59
  60    ! [1.1]
  61
  62    phi_b  =   grid%xb%p
  63
  64    do k = kts, kte
  65
  66       term_x(ims:ime,jms:jme) = 0.0
  67       term_y(ims:ime,jms:jme) = 0.0
  68       phi_b_x(ims:ime,jms:jme) = 0.0
  69       phi_b_y(ims:ime,jms:jme) = 0.0
  70
  71       !---------------------------------------------------------------------------
  72       ! [2.0] Calculate RHS of balance equation in gridpoint space:
  73       !---------------------------------------------------------------------------
  74
  75       ! [2.1] Include geostrophic terms in balance eqn if requested:
  76
  77       if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
  78          call da_wpec_constraint_geoterm (grid%xb % cori, grid%xb % rho(:,:,k), grid%xb %u(:,:,k), grid%xb %v(:,:,k), &
  79             term_x, term_y)
  80       end if
  81
  82       ! [2.2] Include cyclostrophic terms in balance eqn if requested:
  83
  84       if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
  85          call da_wpec_constraint_cycloterm (grid%xb % rho(:,:,k), grid%xb % u(:,:,k),   &
  86             grid%xb % v(:,:,k), grid%xb % coefx(:,:), grid%xb % coefy(:,:), &
  87             term_x, term_y)
  88       end if
  89
  90       ! [2.3] Include phi_b terms in balance eqn
  91       do j = js, je
  92          do i = is, ie
  93             phi_b_x(i,j) = coefx(i,j)*(phi_b(i+1,j,k)-phi_b(i-1,j,k))
  94             phi_b_y(i,j) = coefy(i,j)*(phi_b(i,j+1,k)-phi_b(i,j-1,k))
  95          end do
  96       end do
  97
  98    !------------------------------------------------------------------------------
  99    !  [3.0] Solve Grad_p for balance eqn :
 100    !------------------------------------------------------------------------------
 101       do j = js, je
 102          do i = is, ie
 103             grid%xb%xb_p_x(i,j,k)=phi_b_x(i,j)+term_x(i,j)
 104             grid%xb%xb_p_y(i,j,k)=phi_b_y(i,j)+term_y(i,j)
 105          end do
 106       end do
 107
 108     end do
 109
 110
 111    if (trace_use) call da_trace_exit("da_wpec_constraint")
 112
 113 end subroutine da_wpec_constraint
 114
 115