forces Subroutine

  subroutine forces

    !-----------------------------------------------------------------|
    !                                                                 |
    !      Hans Cuijpers   I.M.A.U.                                   |
    !      Pier Siebesma   K.N.M.I.     06/01/1995                    |
    !                                                                 |
    !     purpose.                                                    |
    !     --------                                                    |
    !                                                                 |
    !      Calculates all other terms in the N-S equation,            |
    !      except for the diffusion and advection terms.              |
    !                                                                 |
    !**   interface.                                                  |
    !     ----------                                                  |
    !                                                                 |
    !     *forces* is called from *program*.                          |
    !                                                                 |
    !-----------------------------------------------------------------|

    !  use modglobal, only : i1,j1,kmax,dzh,dzf,grav
    use modglobal, only : ib,ie,jb,je,kb,ke,kh,dzhi,dzf,grav,lbuoyancy
    use modfields, only : u0,v0,w0,up,vp,wp,thv0h,dpdxl,dpdyl,thlp,thlpcar,thvh
    use modibmdata, only : nxwallsnorm, xwallsnorm
    use modsurfdata,only : thvs
    use modmpi, only : myid
    implicit none

    real thvsi
    integer i, j, k, n, jm, jp, km, kp


    if (lbuoyancy ) then
    !ILS13 replace thvsi by thvh
    ! thvsi = 1./thvsi

 
       !write(*,*) 'thvh',thvh
       do k=kb+1,ke
          do j=jb,je
             do i=ib,ie
                up(i,j,k) = up(i,j,k) - dpdxl(k)
                vp(i,j,k) = vp(i,j,k) - dpdyl(k)
                wp(i,j,k) = wp(i,j,k) + grav * (thv0h(i,j,k)-thvh(k))/thvh(k)
             end do
          end do
       end do

    else

       do k=kb+1,ke
          do j=jb,je
             do i=ib,ie
                up(i,j,k) = up(i,j,k) - dpdxl(k)
                vp(i,j,k) = vp(i,j,k) - dpdyl(k)
                ! IS+HJ      wp(i,j,k) = wp(i,j,k)
             end do
          end do
       end do
    end if

    !     ----------------------------------------------
    !     add radiative heating to potential temperature
    !     ----------------------------------------------
    do k=kb,ke
       do j=jb,je
          do i=ib,ie
             thlp(i,j,k) = thlp(i,j,k)+thlpcar(k)
          end do
       end do
    end do

    !     --------------------------------------------
    !     special treatment for lowest full level: k=1
    !     --------------------------------------------

    do j=jb,je
       jp = j+1
       jm = j-1
       do i=ib,ie

          up(i,j,kb) = up(i,j,kb) - dpdxl(kb)

          vp(i,j,kb) = vp(i,j,kb) - dpdyl(kb)

          wp(i,j,kb) = 0.0

       end do
    end do
    !     ----------------------------------------------end i,j-loop


    return
  end subroutine forces