subroutine da_balance_equation_adj(grid, xbx, phi_b, u, v) !--------------------------------------------------------------------------- ! Purpose: Adjoint of da_balance_equation !--------------------------------------------------------------------------- implicit none type(domain), intent(inout) :: grid type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars. real, intent(in) :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment. real, intent(inout) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp. (dot pts) real, intent(inout) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp. (dot pts) integer :: i, j, k ! Loop counters. integer :: is, ie ! 1st dim. end points. integer :: js, je ! 2nd dim. end points. real, dimension(ims:ime,jms:jme) :: coefx, & ! Multiplicative coefficient. coefy, & ! Multiplicative coefficient. term_x, & ! Balance eqn x term term_y ! Balance eqn y term real, dimension(ims:ime,jms:jme,kms:kme) :: del2phi_b ! Del**2 phi_b/M**2 real :: coeff1, coeff2 ! Multiplicative coefficient. if (trace_use) call da_trace_entry("da_balance_equation_adj") !--------------------------------------------------------------------------- ! [1.0] Initialise: !--------------------------------------------------------------------------- ! Computation to check for edge of domain: is = its-1; ie = ite+1; js = jts-1; je = jte+1 if (.not. global .and. its == ids) is = ids+1 if (.not. global .and. ite == ide) ie = ide-1 if (jts == jds ) js = jds+1 if (jte == jde ) je = jde-1 if (fg_format == fg_format_kma_global) then coefx(is:ie,js:je) = grid%xb%coefx(is:ie,js:je) coefy(is:ie,js:je) = grid%xb%coefy(is:ie,js:je) else if (fg_format == fg_format_wrf_arw_global) then write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format = ',fg_format call da_error(__FILE__,__LINE__,message(1:1)) else if (fg_format == fg_format_wrf_arw_regional) then coefx(is:ie,js:je) = grid%xb%coefz(is:ie,js:je) coefy(is:ie,js:je) = coefx(is:ie,js:je) else if (fg_format == fg_format_wrf_nmm_regional) then write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format = ',fg_format call da_error(__FILE__,__LINE__,message(1:1)) else write (unit=message(1),fmt='(A,I3)') ' Wrong FG_FORMAT = ',fg_format call da_error(__FILE__,__LINE__,message(1:1)) end if ! [1.1] Multiplicative coefficent for conversion RHS->Del**2 phi_b/M**2: del2phi_b(:,:,:) = 0.0 !--------------------------------------------------------------------------- ! [3.0] Solve Del**2 phi_b = RHS for phi_b: !--------------------------------------------------------------------------- call da_solve_poissoneqn_fst_adj(grid,xbx, phi_b, del2phi_b) !--------------------------------------------------------------------------- ! [2.0] Calculate RHS of balance equation in gridpt space: !--------------------------------------------------------------------------- do k = kts, kte ! [2.4] Del**2 Phi_b boundary conditions (null as zero boundary conditions): ! [2.3] Take divergence to get Del**2 phi_b/M**2: term_x(ims:ime,jms:jme) = 0.0 term_y(ims:ime,jms:jme) = 0.0 do j = je, js, -1 do i = ie, is, -1 coeff1 = coefx(i,j) * del2phi_b(i,j,k) coeff2 = coefy(i,j) * del2phi_b(i,j,k) term_x(i+1,j) = term_x(i+1,j) - coeff1 term_x(i-1,j) = term_x(i-1,j) + coeff1 term_y(i,j+1) = term_y(i,j+1) - coeff2 term_y(i,j-1) = term_y(i,j-1) + coeff2 end do end do ! [2.2] Include cyclostrophic terms in balance eqn if requested: if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then call da_balance_cycloterm_adj (grid%xb%rho(:,:,k),grid%xb%u(:,:,k),& grid%xb%v(:,:,k), u(:,:,k), v(:,:,k), grid%xb%coefx(:,:), grid%xb%coefy(:,:),& term_x(:,:), term_y(:,:)) end if ! [2.1] Calculate geostrophic terms in balance eqn: if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then call da_balance_geoterm_adj (grid%xb%cori, grid%xb%rho(:,:,k), term_x, term_y, & u(:,:,k), v(:,:,k)) end if end do if (trace_use) call da_trace_exit("da_balance_equation_adj") end subroutine da_balance_equation_adj