external/RSL_LITE/force_domain_em_part2.F

   1 #if ( EM_CORE == 1 && DA_CORE != 1 )
   2
   3 !------------------------------------------------------------------
   4
   5    SUBROUTINE force_domain_em_part2 ( grid, ngrid, pgrid, config_flags    &
   6 !
   7 #include "dummy_new_args.inc"
   8 !
   9                  )
  10       USE module_state_description
  11       USE module_domain, ONLY : domain, get_ijk_from_grid
  12       USE module_configure, ONLY : grid_config_rec_type
  13       USE module_dm, ONLY : ntasks, ntasks_x, ntasks_y, local_communicator, mytask, &
  14                             nest_pes_x, nest_pes_y ! ,                                 &
  15                             !push_communicators_for_domain,pop_communicators_for_domain
  16       USE module_comm_nesting_dm, ONLY : halo_force_down_sub
  17       USE module_model_constants
  18       IMPLICIT NONE
  19 !
  20       TYPE(domain), POINTER :: grid          ! name of the grid being dereferenced (must be "grid")
  21       TYPE(domain), POINTER :: ngrid
  22       TYPE(domain), POINTER :: pgrid         !KAL added for vertical nesting
  23 #include "dummy_new_decl.inc"
  24       INTEGER nlev, msize
  25       INTEGER i,j,pig,pjg,cm,cn,nig,njg,retval,k,kk
  26       TYPE (grid_config_rec_type)            :: config_flags
  27       REAL xv(2000)
  28       INTEGER       ::          cids, cide, cjds, cjde, ckds, ckde,    &
  29                                 cims, cime, cjms, cjme, ckms, ckme,    &
  30                                 cips, cipe, cjps, cjpe, ckps, ckpe
  31       INTEGER       ::          nids, nide, njds, njde, nkds, nkde,    &
  32                                 nims, nime, njms, njme, nkms, nkme,    &
  33                                 nips, nipe, njps, njpe, nkps, nkpe
  34       INTEGER       ::          ids, ide, jds, jde, kds, kde,    &
  35                                 ims, ime, jms, jme, kms, kme,    &
  36                                 ips, ipe, jps, jpe, kps, kpe
  37       INTEGER idim1,idim2,idim3,idim4,idim5,idim6,idim7,itrace
  38       REAL  dummy_xs, dummy_xe, dummy_ys, dummy_ye
  39
  40       !KAL variables for vertical nesting
  41       REAL :: p_top_m  , p_surf_m , mu_m , hsca_m , pre_c ,pre_n
  42       REAL, DIMENSION(pgrid%s_vert:pgrid%e_vert) :: alt_w_c
  43       REAL, DIMENSION(pgrid%s_vert:pgrid%e_vert+1) :: alt_u_c
  44       REAL, DIMENSION(ngrid%s_vert:ngrid%e_vert) :: alt_w_n
  45       REAL, DIMENSION(ngrid%s_vert:ngrid%e_vert+1) :: alt_u_n
  46
  47       REAL, DIMENSION(:,:,:), ALLOCATABLE :: p, al
  48       REAL :: pfu, pfd, phm, temp, qvf, qvf1, qvf2
  49
  50       !KAL change this for vertical nesting
  51       ! force_domain_em_part1 packs up the interpolation onto the coarse (vertical) grid
  52       ! therefore the message size is based on the coarse grid number of levels
  53       ! here it is unpacked onto the intermediate grid
  54       CALL get_ijk_from_grid (  pgrid ,                   &
  55                                 cids, cide, cjds, cjde, ckds, ckde,    &
  56                                 cims, cime, cjms, cjme, ckms, ckme,    &
  57                                 cips, cipe, cjps, cjpe, ckps, ckpe    )
  58
  59       !KAL this is the original WRF code
  60       !CALL get_ijk_from_grid (  grid ,                   &
  61       !                          cids, cide, cjds, cjde, ckds, ckde,    &
  62       !                          cims, cime, cjms, cjme, ckms, ckme,    &
  63       !                          cips, cipe, cjps, cjpe, ckps, ckpe    )
  64       CALL get_ijk_from_grid (  ngrid ,              &
  65                                 nids, nide, njds, njde, nkds, nkde,    &
  66                                 nims, nime, njms, njme, nkms, nkme,    &
  67                                 nips, nipe, njps, njpe, nkps, nkpe    )
  68
  69       nlev  = ckde - ckds + 1
  70
  71 #include "nest_interpdown_unpack.inc"
  72
  73 if (ngrid%vert_refine_method .NE. 0) then
  74
  75       !KAL calculating the vertical coordinate for parent and nest grid (code from ndown)
  76       ! assume that the parent and nest have the same p_top value (as in ndown)
  77
  78 !KAL ckde is equal to e_vert of the coarse grid. There are e_vert-1 u points.  The coarse 1D grid here is e_vert+1,
  79 !    so it is the e_vert-1 points from the coarse grid, plus a surface point plus a top point.  Extrapolation coefficients
  80 !    are used to get the surface and top points to fill out the pro_u_c 1D array of u values from the coarse grid.
  81
  82       hsca_m = 6.7 !KAL scale height of the atmosphere
  83       p_top_m = ngrid%p_top
  84       p_surf_m = 1.e5
  85       mu_m = p_surf_m - p_top_m
  86 !    parent
  87       do  k = 1,ckde
  88       pre_c = mu_m * pgrid%c3f(k) + p_top_m + pgrid%c4f(k)
  89       alt_w_c(k) =  -hsca_m * alog(pre_c/p_surf_m)
  90       enddo
  91       do  k = 1,ckde-1
  92       pre_c = mu_m * pgrid%c3h(k) + p_top_m + pgrid%c4h(k)
  93       alt_u_c(k+1) =  -hsca_m * alog(pre_c/p_surf_m)
  94       enddo
  95       alt_u_c(1) =  alt_w_c(1)
  96       alt_u_c(ckde+1) =  alt_w_c(ckde)
  97 !    nest
  98       do  k = 1,nkde
  99       pre_n = mu_m * ngrid%c3f(k) + p_top_m + ngrid%c4f(k)
 100       alt_w_n(k) =  -hsca_m * alog(pre_n/p_surf_m)
 101       enddo
 102       do  k = 1,nkde-1
 103       pre_n = mu_m * ngrid%c3h(k) + p_top_m + ngrid%c4h(k)
 104       alt_u_n(k+1) =  -hsca_m * alog(pre_n/p_surf_m)
 105       enddo
 106       alt_u_n(1) =  alt_w_n(1)
 107       alt_u_n(nkde+1) =  alt_w_n(nkde)
 108
 109 endif
 110
 111       !KAL added this call for vertical nesting (return coarse grid dimensions to intended values)
 112       CALL get_ijk_from_grid (  grid ,                   &
 113                                 cids, cide, cjds, cjde, ckds, ckde,    &
 114                                 cims, cime, cjms, cjme, ckms, ckme,    &
 115                                 cips, cipe, cjps, cjpe, ckps, ckpe    )
 116
 117       CALL get_ijk_from_grid (  grid ,              &
 118                                 ids, ide, jds, jde, kds, kde,    &
 119                                 ims, ime, jms, jme, kms, kme,    &
 120                                 ips, ipe, jps, jpe, kps, kpe    )
 121
 122       !  Vertical refinement is turned on.
 123
 124       IF (ngrid%vert_refine_method .NE. 0) THEN
 125
 126 #include "nest_forcedown_interp_vert.inc"
 127
 128          IF ( ngrid%this_is_an_ideal_run ) THEN
 129             IF ( SIZE( grid%t_init, 1 ) * SIZE( grid%t_init, 3 ) .GT. 1 ) THEN
 130                CALL vert_interp_vert_nesting( grid%t_init, & !CD field
 131                                               ids, ide, kds, kde, jds, jde, & !CD dims
 132                                               ims, ime, kms, kme, jms, jme, & !CD dims
 133                                               ips, ipe, kps, MIN( (kde-1), kpe ), jps, jpe, & !CD dims
 134                                               pgrid%s_vert, pgrid%e_vert, & !vertical dimension of the parent grid
 135                                               pgrid%cf1, pgrid%cf2, pgrid%cf3, pgrid%cfn, pgrid%cfn1, & !coarse grid extrapolation constants
 136                                               alt_u_c, alt_u_n ) !coordinates for parent and nest
 137             END IF  ! Check t_init is a fully allocated 3d array.
 138          END IF ! only for ideal runs
 139
 140
 141          !  Rebalance the grid on the intermediate grid.  The intermediate grid has the horizontal
 142          !  resolution of the parent grid, but at this point has been interpolated in the vertical
 143          !  to the resolution of the nest.  The base state (phb, pb, etc) from the parent grid is
 144          !  unpacked onto the intermediate grid every time this subroutine is called.  We need the
 145          !  base state of the nest, so it is recalculated here.
 146
 147          !  Additionally, we do not need to vertically interpolate the entire intermediate grid
 148          !  above, just the points that contribute to the boundary forcing.
 149
 150          !  Base state potential temperature and inverse density (alpha = 1/rho) from
 151          !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
 152          !  from equation of state.  The potential temperature is a perturbation from t0.
 153
 154          !  Uncouple the variables moist and t_2 that are used to calculate ph_2
 155
 156          DO j = MAX(jds,jps),MIN(jde-1,jpe)
 157             DO i = MAX(ids,ips),MIN(ide-1,ipe)
 158                DO k=kds,kde-1
 159                   grid%t_2(i,k,j) = grid%t_2(i,k,j)/((ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k)) + (ngrid%c1h(k)*grid%mu_2(i,j)))
 160                   moist(i,k,j,P_QV) = moist(i,k,j,P_QV)/((ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k)) + (ngrid%c1h(k)*grid%mu_2(i,j)))
 161                END DO
 162             END DO
 163          END DO
 164
 165          DO j = MAX(jds,jps),MIN(jde-1,jpe)
 166             DO i = MAX(ids,ips),MIN(ide-1,ipe)
 167
 168                DO k = 1, kpe-1
 169                   grid%pb(i,k,j) = ngrid%c3h(k) * grid%mub(i,j) + ngrid%c4h(k) + ngrid%p_top
 170
 171                   !  If this is a real run, recalc t_init.
 172
 173                   IF ( .NOT. ngrid%this_is_an_ideal_run ) THEN
 174                      temp = MAX ( ngrid%tiso, ngrid%t00 + ngrid%tlp*LOG(grid%pb(i,k,j)/ngrid%p00) )
 175                      IF ( grid%pb(i,k,j) .LT. ngrid%p_strat ) THEN
 176                         temp = ngrid%tiso + ngrid%tlp_strat * LOG ( grid%pb(i,k,j)/ngrid%p_strat )
 177                      END IF
 178                      grid%t_init(i,k,j) = temp*(ngrid%p00/grid%pb(i,k,j))**(r_d/cp) - t0
 179                   END IF
 180                   grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 181                END DO
 182
 183                !  Integrate base geopotential, starting at terrain elevation.  This assures that
 184                !  the base state is in exact hydrostatic balance with respect to the model equations.
 185                !  This field is on full levels.
 186
 187                grid%phb(i,1,j) = grid%ht(i,j) * g
 188                IF (grid%hypsometric_opt == 1) THEN
 189                   DO kk = 2,kpe
 190                      k = kk - 1
 191                      grid%phb(i,kk,j) = grid%phb(i,k,j) - ngrid%dnw(k)*(ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k))*grid%alb(i,k,j)
 192                   END DO
 193                ELSE IF (grid%hypsometric_opt == 2) THEN
 194                   DO k = 2,kpe
 195                      pfu = ngrid%c3f(k  )*grid%MUB(i,j) + ngrid%c4f(k  ) + ngrid%p_top
 196                      pfd = ngrid%c3f(k-1)*grid%MUB(i,j) + ngrid%c4f(k-1) + ngrid%p_top
 197                      phm = ngrid%c3h(k-1)*grid%MUB(i,j) + ngrid%c4h(k-1) + ngrid%p_top
 198                      grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
 199                   END DO
 200                ELSE
 201                   CALL wrf_error_fatal( 'module_dm: hypsometric_opt should be 1 or 2' )
 202                END IF  ! which hypsometric option
 203             END DO  ! i loop
 204          END DO  ! j loop
 205          !  Perturbation fields
 206          ALLOCATE( p (ips:ipe, kps:kpe, jps:jpe) )
 207          ALLOCATE( al(ips:ipe, kps:kpe, jps:jpe) )
 208          DO j = MAX(jds,jps),MIN(jde-1,jpe)
 209             DO i = MAX(ids,ips),MIN(ide-1,ipe)
 210                !  Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
 211                !  equation) down from the top to get the pressure perturbation.  First get the pressure
 212                !  perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
 213
 214                kk = kpe-1
 215                k = kk+1
 216
 217                qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk,j,P_QV))
 218                qvf2 = 1./(1.+qvf1)
 219                qvf1 = qvf1*qvf2
 220
 221                p(i,kk,j) = - 0.5*((ngrid%c1f(k)*grid%Mu_2(i,j))+qvf1*(ngrid%c1f(k)*grid%Mub(i,j)+ngrid%c2f(k)))/ngrid%rdnw(kk)/qvf2
 222                IF ( config_flags%use_theta_m == 0) THEN
 223                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
 224                ELSE
 225                   qvf = 1.
 226                ENDIF
 227                al(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
 228                            (((p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm) - grid%alb(i,kk,j)
 229
 230                !  Now, integrate down the column to compute the pressure perturbation, and diagnose the two
 231                !  inverse density fields (total and perturbation).
 232
 233                DO kk=kpe-2,1,-1
 234                   k = kk + 1
 235                   qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk+1,j,P_QV))
 236                   qvf2 = 1./(1.+qvf1)
 237                   qvf1 = qvf1*qvf2
 238                   p(i,kk,j) = p(i,kk+1,j) - ((ngrid%c1f(k)*grid%Mu_2(i,j)) + qvf1*(ngrid%c1f(k)*grid%Mub(i,j)+ngrid%c2f(k)))/qvf2/ngrid%rdn(kk+1)
 239                   IF ( config_flags%use_theta_m == 0) THEN
 240                      qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
 241                   ELSE
 242                      qvf = 1.
 243                   ENDIF
 244                   al(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
 245                               (((p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm) - grid%alb(i,kk,j)
 246                END DO
 247
 248                !  This is the hydrostatic equation used in the model after the small timesteps.  In
 249                !  the model, grid%al (inverse density) is computed from the geopotential.
 250
 251                IF (grid%hypsometric_opt == 1) THEN
 252                   DO kk = 2,kpe
 253                      k = kk - 1
 254                      grid%ph_2(i,kk,j) = grid%ph_2(i,kk-1,j) - &
 255                                         ngrid%dnw(kk-1) * ( ((ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k))+(ngrid%c1h(k)*grid%mu_2(i,j)))*al(i,kk-1,j) &
 256                                         + (ngrid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,kk-1,j) )
 257                   END DO
 258
 259                ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
 260                ! Note that al*p approximates Rd*T and dLOG(p) does z.
 261                ! Here T varies mostly linear with z, the first-order integration produces better result.
 262
 263                ELSE IF (grid%hypsometric_opt == 2) THEN
 264
 265                   grid%ph_2(i,1,j) = grid%phb(i,1,j)
 266                   DO k = 2,kpe
 267                      pfu = ngrid%c3f(k  )*( grid%MUB(i,j)+grid%MU_2(i,j) ) + ngrid%c4f(k  ) + ngrid%p_top
 268                      pfd = ngrid%c3f(k-1)*( grid%MUB(i,j)+grid%MU_2(i,j) ) + ngrid%c4f(k-1) + ngrid%p_top
 269                      phm = ngrid%c3h(k-1)*( grid%MUB(i,j)+grid%MU_2(i,j) ) + ngrid%c4h(k-1) + ngrid%p_top
 270                      grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + (grid%alb(i,k-1,j)+al(i,k-1,j))*phm*LOG(pfd/pfu)
 271                   END DO
 272
 273                   DO k = 1,kpe
 274                      grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
 275                   END DO
 276
 277                END IF
 278
 279             END DO ! i loop
 280          END DO ! j loop
 281
 282          DEALLOCATE(p)
 283          DEALLOCATE(al)
 284
 285          ! Couple the variables moist and t_2, and the newly calculated ph_2
 286          DO j = MAX(jds,jps),MIN(jde-1,jpe)
 287             DO i = MAX(ids,ips),MIN(ide-1,ipe)
 288                DO k=kps,kpe
 289                grid%ph_2(i,k,j) = grid%ph_2(i,k,j)*((ngrid%c1f(k)*grid%Mub(i,j)+ngrid%c2f(k)) + (ngrid%c1f(k)*grid%Mu_2(i,j)))
 290                END DO
 291             END DO
 292          END DO
 293          DO j = MAX(jds,jps),MIN(jde-1,jpe)
 294             DO i = MAX(ids,ips),MIN(ide-1,ipe)
 295                DO k=kps,kpe-1
 296                grid%t_2(i,k,j) = grid%t_2(i,k,j)*((ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k)) + (ngrid%c1h(k)*grid%mu_2(i,j)))
 297                moist(i,k,j,P_QV) = moist(i,k,j,P_QV)*((ngrid%c1h(k)*grid%mub(i,j)+ngrid%c2h(k)) + (ngrid%c1h(k)*grid%mu_2(i,j)))
 298                END DO
 299             END DO
 300          END DO
 301
 302
 303       END IF
 304
 305
 306 #include "HALO_FORCE_DOWN.inc"
 307
 308       ! code here to interpolate the data into the nested domain
 309 # include "nest_forcedown_interp.inc"
 310
 311       RETURN
 312    END SUBROUTINE force_domain_em_part2
 313
 314 #endif