var/da/da_physics/da_uvprho_to_w_adj.inc

   1 subroutine da_uvprho_to_w_adj(grid)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Calculates vertical velocity increments from Richardson's Eq.
   5    !
   6    ! Method: Richardson's Eq., which
   7    !         combines continuity Eq., thermodynamic Eq. and hrdrostatic Eq.
   8    !---------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    type (domain), intent(inout)     :: grid
  13
  14    integer                       :: is, ie       ! 1st dim. end points.
  15    integer                       :: js, je       ! 2nd dim. end points.
  16
  17    integer                       :: I,J,K
  18
  19    real, dimension(ims:ime,jms:jme,kms:kme) :: URHO, VRHO
  20    real, dimension(ims:ime,jms:jme,kms:kme) :: DIV, WZ
  21    real                                     :: TERM3
  22
  23    if (trace_use) call da_trace_entry("da_uvprho_to_w_adj")
  24
  25    ! initialize to zero for some variables because of the adjoint requirements.
  26
  27    WZ(:,:,:)   = 0.0
  28    URHO(:,:,:) = 0.0
  29    VRHO(:,:,:) = 0.0
  30    DIV(:,:,:)  = 0.0
  31    TERM3       = 0.0
  32
  33    ! integration to calculate the vertical velocity
  34
  35    do J=jts,jte
  36       do I=its,ite
  37          do K=kts,kte
  38             grid%xa%w(I,J,K+1)=grid%xa%w(I,J,K+1)+grid%xa%w(I,J,K)
  39             WZ(I,J,K)=grid%xa%w(I,J,K)*(grid%xb%hf(I,J,K)-grid%xb%hf(I,J,K+1))
  40             grid%xa%w(I,J,K)=0.0
  41          end do
  42          grid%xa%w(I,J,kte+1)=0.0
  43       end do
  44    end do
  45
  46    call da_w_adjustment_adj(grid%xb,WZ)
  47
  48    ! Divide by constant
  49
  50    WZ(its:ite,jts:jte,kts:kte)=WZ(its:ite,jts:jte,kts:kte)/(GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte))
  51
  52    ! Term 4: Derivative of basic vertical velocity with respect to z.
  53
  54    do J=jts,jte
  55       do I=its,ite
  56          do K=kte,kts,-1
  57             grid%xa%p(I,J,K)=grid%xa%p(I,J,K)-WZ(I,J,K)*GAMMA*  &
  58                      (grid%xb%w(I,J,K+1)-grid%xb%w(I,J,K))/  &
  59                      (grid%xb%hf(I,J,K+1)-grid%xb%hf(I,J,K))
  60          end do
  61       end do
  62    end do
  63
  64    ! Term 3.2: Vertical integration of the basic mass divergence
  65
  66    do J=jts,jte
  67       do I=its,ite
  68          do K=kts,kte-1
  69             TERM3=TERM3+WZ(I,J,K)
  70             DIV(I,J,K+1)=DIV(I,J,K+1)+  &
  71                       TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
  72             DIV(I,J,K)  =DIV(I,J,K)+  &
  73                       TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
  74          end do
  75         TERM3=0.0
  76       end do
  77    end do
  78
  79    call da_uv_to_divergence_adj(grid, URHO,VRHO, DIV)
  80
  81    ! Computation to check for edge of domain:
  82    if (test_transforms) then
  83       is = its-1; ie = ite+1; js = jts-1; je = jte+1
  84       if (its == ids) is = ids; if (ite == ide) ie = ide
  85       if (jts == jds) js = jds; if (jte == jde) je = jde
  86    else
  87       is = its
  88       ie = ite
  89       js = jts
  90       je = jte
  91    end if
  92
  93    grid%xa%rho(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)+VRHO(is:ie,js:je,kts:kte)*grid%xb%v(is:ie,js:je,kts:kte)
  94    grid%xa%rho(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)+URHO(is:ie,js:je,kts:kte)*grid%xb%u(is:ie,js:je,kts:kte)
  95
  96    URHO(:,:,:)=0.0
  97    VRHO(:,:,:)=0.0
  98
  99    ! Term 3.1: Vertical integration of the perturbed mass divergence
 100
 101    do J=jts,jte
 102       do I=its,ite
 103          do K=kts,kte-1
 104             TERM3=TERM3+WZ(I,J,K)
 105             DIV(I,J,K+1)=DIV(I,J,K+1)+  &
 106                       TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
 107             DIV(I,J,K)  =DIV(I,J,K)+  &
 108                       TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
 109          end do
 110          TERM3=0.0
 111       end do
 112    end do
 113
 114    call da_uv_to_divergence_adj(grid, URHO,VRHO, DIV)
 115
 116    grid%xa%v(is:ie,js:je,kts:kte)=grid%xa%v(is:ie,js:je,kts:kte)+VRHO(is:ie,js:je,kts:kte)*grid%xb%rho(is:ie,js:je,kts:kte)
 117    grid%xa%u(is:ie,js:je,kts:kte)=grid%xa%u(is:ie,js:je,kts:kte)+URHO(is:ie,js:je,kts:kte)*grid%xb%rho(is:ie,js:je,kts:kte)
 118
 119    URHO(:,:,:)=0.0
 120    VRHO(:,:,:)=0.0
 121
 122    ! Term 2.2: Divergence term from basic wind
 123
 124    call da_uv_to_divergence(grid%xb, grid%xb%u,grid%xb%v, DIV)
 125
 126    grid%xa%p(its:ite,jts:jte,kts:kte)=grid%xa%p(its:ite,jts:jte,kts:kte)-WZ(its:ite,jts:jte,kts:kte)*GAMMA*DIV(its:ite,jts:jte,kts:kte)
 127
 128    ! Term 2.1: Divergence term from perturbed wind
 129
 130    DIV(its:ite,jts:jte,kts:kte)=-WZ(its:ite,jts:jte,kts:kte)*GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte)  ! DIV redefined
 131
 132    call da_uv_to_divergence_adj(grid, grid%xa%u,grid%xa%v, DIV)
 133
 134    ! Computation to check for edge of domain:
 135    is = its
 136    ie = ite
 137    js = jts
 138    je = jte
 139    if (its == ids) is = ids+1
 140    if (ite == ide) ie = ide-1
 141    if (jts == jds) js = jds+1
 142    if (jte == jde) je = jde-1
 143
 144    ! Term 1.2: Basic pressure advection along the perturbed wind
 145
 146    if (jte == jde) then
 147       j = jte
 148       do K=kts,kte
 149          do I=its, ite
 150             WZ(I,J-1,K)=WZ(I,J-1,K)+WZ(I,J,K)
 151          end do
 152       end do
 153    end if
 154
 155    if (jts == jds) then
 156       j = jts
 157       do K=kts,kte
 158          do I=its, ite
 159             WZ(I,J+1,K)=WZ(I,J+1,K)+WZ(I,J,K)
 160          end do
 161       end do
 162    end if
 163
 164    if (ite == ide) then
 165       i = ite
 166       do K=kts,kte
 167          do J=js,je
 168             WZ(I-1,J,K)=WZ(I-1,J,K)+WZ(I,J,K)
 169          end do
 170       end do
 171    end if
 172
 173    if (its == ids) then
 174       i = its
 175       do K=kts,kte
 176          do J=js,je
 177             WZ(I+1,J,K)=WZ(I+1,J,K)+WZ(I,J,K)
 178          end do
 179       end do
 180    end if
 181
 182    do K=kts,kte
 183       do J=js,je
 184          do I=is,ie
 185             grid%xa%v(I,J,K)=grid%xa%v(I,J,K)-WZ(I,J,K)*  &
 186                (grid%xb%p(I,J+1,K)-grid%xb%p(I,J-1,K))*grid%xb%coefy(I,J)
 187             grid%xa%u(I,J,K)=grid%xa%u(I,J,K)-WZ(I,J,K)*  &
 188                (grid%xb%p(I+1,J,K)-grid%xb%p(I-1,J,K))*grid%xb%coefx(I,J)
 189          end do
 190       end do
 191    end do
 192
 193    !-------------------------------------------------------------------------
 194    ! Computation to check for edge of domain:
 195    ! This is only for adjoint, as we have to cross the processor boundary
 196    ! to get the contribution.
 197
 198    is = its - 1
 199    ie = ite + 1
 200    js = jts - 1
 201    je = jte + 1
 202
 203    grid%xp%v1z(its:ite, jts:jte, kts:kte) = wz(its:ite, jts:jte, kts:kte)
 204
 205 #ifdef DM_PARALLEL
 206 #include "HALO_BAL_EQN_ADJ.inc"
 207 #endif
 208
 209    if (its == ids) then
 210       is = ids+1
 211    else
 212       wz(is, js:je, kts:kte) = grid%xp%v1z(is, js:je, kts:kte)
 213    end if
 214
 215    if (ite == ide) then
 216       ie = ide-1
 217    else
 218       wz(ie, js:je, kts:kte) = grid%xp%v1z(ie, js:je, kts:kte)
 219    end if
 220
 221    if (jts == jds) then
 222       js = jds+1
 223    else
 224       wz(is:ie, js, kts:kte) = grid%xp%v1z(is:ie, js, kts:kte)
 225    end if
 226
 227    if (jte == jde) then
 228       je = jde-1
 229    else
 230       wz(is:ie, je, kts:kte) = grid%xp%v1z(is:ie, je, kts:kte)
 231    end if
 232
 233    ! Term 1.1: Perturbed pressure advection along the basic wind
 234
 235    do K=kts,kte
 236       do J=js,je
 237          do I=is,ie
 238             grid%xa%p(I,J+1,K)=grid%xa%p(I,J+1,K)-WZ(I,J,K)*grid%xb%v(I,J,K)*grid%xb%coefy(I,J)
 239             grid%xa%p(I,J-1,K)=grid%xa%p(I,J-1,K)+WZ(I,J,K)*grid%xb%v(I,J,K)*grid%xb%coefy(I,J)
 240             grid%xa%p(I+1,J,K)=grid%xa%p(I+1,J,K)-WZ(I,J,K)*grid%xb%u(I,J,K)*grid%xb%coefx(I,J)
 241             grid%xa%p(I-1,J,K)=grid%xa%p(I-1,J,K)+WZ(I,J,K)*grid%xb%u(I,J,K)*grid%xb%coefx(I,J)
 242          end do
 243       end do
 244    end do
 245
 246    WZ(:,:,:) = 0.0
 247
 248    if (trace_use) call da_trace_exit("da_uvprho_to_w_adj")
 249
 250 end subroutine da_uvprho_to_w_adj
 251
 252