wrfv2_fire/dyn_em/module_force_scm.F

   1 MODULE module_force_scm
   2
   3 ! AUTHOR: Josh Hacker (NCAR/RAL)
   4 ! Forces a single-column (3x3) version of WRF
   5
   6 CONTAINS
   7
   8    SUBROUTINE force_scm(itimestep, dt, scm_force, dx, num_force_layers       &
   9                              , scm_th_adv, scm_qv_adv                        &
  10                              , scm_ql_adv                                    &
  11                              , scm_wind_adv, scm_vert_adv                    &
  12                              , scm_soilT_force, scm_soilQ_force              &
  13                              , scm_force_th_largescale                       &
  14                              , scm_force_qv_largescale                       &
  15                              , scm_force_ql_largescale                       &
  16                              , scm_force_wind_largescale                     &
  17                              , u_base, v_base, z_base                        &
  18                              , z_force, z_force_tend                         &
  19                              , u_g, v_g                                      &
  20                              , u_g_tend, v_g_tend                            &
  21                              , w_subs, w_subs_tend                           &
  22                              , th_upstream_x, th_upstream_x_tend             &
  23                              , th_upstream_y, th_upstream_y_tend             &
  24                              , qv_upstream_x, qv_upstream_x_tend             &
  25                              , qv_upstream_y, qv_upstream_y_tend             &
  26                              , ql_upstream_x, ql_upstream_x_tend             &
  27                              , ql_upstream_y, ql_upstream_y_tend             &
  28                              , u_upstream_x, u_upstream_x_tend               &
  29                              , u_upstream_y, u_upstream_y_tend               &
  30                              , v_upstream_x, v_upstream_x_tend               &
  31                              , v_upstream_y, v_upstream_y_tend               &
  32                              , tau_x, tau_x_tend                             &
  33                              , tau_y, tau_y_tend                             &
  34                              ,th_largescale                                  &
  35                              ,th_largescale_tend                             &
  36                              ,qv_largescale                                  &
  37                              ,qv_largescale_tend                             &
  38                              ,ql_largescale                                  &
  39                              ,ql_largescale_tend                             &
  40                              ,u_largescale                                   &
  41                              ,u_largescale_tend                              &
  42                              ,v_largescale                                   &
  43                              ,v_largescale_tend                              &
  44                              ,tau_largescale                                 &
  45                              ,tau_largescale_tend                            &
  46                              , num_force_soil_layers, num_soil_layers        &
  47                              , soil_depth_force, zs                          &
  48                              , tslb, smois                                   &
  49                              , t_soil_forcing_val, t_soil_forcing_tend       &
  50                              , q_soil_forcing_val, q_soil_forcing_tend       &
  51                              , tau_soil                                      &
  52                              , z, z_at_w, th, qv, ql, u, v                   &
  53                              , thten, qvten, qlten, uten, vten               &
  54                              , ids, ide, jds, jde, kds, kde                  &
  55                              , ims, ime, jms, jme, kms, kme                  &
  56                              , ips, ipe, jps, jpe, kps, kpe                  &
  57                              , kts, kte                                      &
  58                             )
  59
  60 ! adds forcing to bl tendencies and also to base state/geostrophic winds.
  61
  62    USE module_init_utilities, ONLY : interp_0
  63    IMPLICIT NONE
  64
  65
  66    INTEGER,    INTENT(IN   )                 :: itimestep
  67    INTEGER,    INTENT(IN   )                 :: num_force_layers, scm_force
  68    REAL,       INTENT(IN   )                 :: dt,dx
  69    LOGICAL,    INTENT(IN   )                 :: scm_th_adv, &
  70                                                 scm_qv_adv, &
  71                                                 scm_ql_adv, &
  72                                                 scm_wind_adv, &
  73                                                 scm_vert_adv, &
  74                                                 scm_soilT_force, &
  75                                                 scm_soilQ_force, &
  76                                                 scm_force_th_largescale, &
  77                                                 scm_force_qv_largescale, &
  78                                                 scm_force_ql_largescale, &
  79                                                 scm_force_wind_largescale
  80
  81    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN   ) :: z, th, qv, ql
  82    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN   ) :: u, v
  83    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN   ) :: z_at_w
  84    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: thten, qvten
  85    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: qlten
  86    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: uten, vten
  87    REAL, DIMENSION( kms:kme ), INTENT(INOUT)               :: u_base, v_base
  88    REAL, DIMENSION( kms:kme ), INTENT(INOUT)               :: z_base
  89    REAL, DIMENSION(num_force_layers), INTENT (INOUT)       :: z_force
  90    REAL, DIMENSION(num_force_layers), INTENT (INOUT)       :: u_g,v_g
  91
  92    REAL, DIMENSION(num_force_layers), INTENT (IN) :: z_force_tend
  93    REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_g_tend,v_g_tend
  94    REAL, DIMENSION(num_force_layers), INTENT (IN) :: w_subs_tend
  95    REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_x_tend
  96    REAL, DIMENSION(num_force_layers), INTENT (IN) :: th_upstream_y_tend
  97    REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_x_tend
  98    REAL, DIMENSION(num_force_layers), INTENT (IN) :: qv_upstream_y_tend
  99    REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_x_tend
 100    REAL, DIMENSION(num_force_layers), INTENT (IN) :: ql_upstream_y_tend
 101    REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_x_tend
 102    REAL, DIMENSION(num_force_layers), INTENT (IN) :: u_upstream_y_tend
 103    REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_x_tend
 104    REAL, DIMENSION(num_force_layers), INTENT (IN) :: v_upstream_y_tend
 105    REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_x_tend
 106    REAL, DIMENSION(num_force_layers), INTENT (IN) :: tau_y_tend
 107
 108    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_x
 109    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_upstream_y
 110    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_x
 111    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_upstream_y
 112    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_x
 113    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_upstream_y
 114    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_x
 115    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_upstream_y
 116    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_x
 117    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_upstream_y
 118    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: w_subs
 119    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_x
 120    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_y
 121
 122 ! WA 1/8/10 for large-scale forcing
 123    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale
 124    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: th_largescale_tend
 125    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale
 126    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: u_largescale_tend
 127    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale
 128    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: v_largescale_tend
 129    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale
 130    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: qv_largescale_tend
 131    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale
 132    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: ql_largescale_tend
 133    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale
 134    REAL, DIMENSION(num_force_layers), INTENT (INOUT) :: tau_largescale_tend
 135
 136 ! WA 1/3/10 For soil forcing
 137    INTEGER,    INTENT(IN   )         :: num_force_soil_layers, num_soil_layers
 138    REAL, DIMENSION(ims:ime,num_soil_layers,jms:jme),INTENT(INOUT) :: tslb, smois
 139    REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_val
 140    REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: t_soil_forcing_tend
 141    REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_val
 142    REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: q_soil_forcing_tend
 143    REAL, DIMENSION(num_force_soil_layers), INTENT (INOUT) :: tau_soil
 144    REAL, DIMENSION(num_force_soil_layers), INTENT (IN   ) :: soil_depth_force
 145    REAL, DIMENSION(num_soil_layers),       INTENT (IN   ) :: zs
 146
 147    INTEGER,    INTENT(IN   )    ::     ids,ide, jds,jde, kds,kde, &
 148                                        ims,ime, jms,jme, kms,kme, &
 149                                        ips,ipe, jps,jpe, kps,kpe, &
 150                                        kts,kte
 151
 152 ! Local
 153    INTEGER                      :: i,j,k
 154    LOGICAL                      :: debug = .false.
 155    REAL                         :: t_x, t_y, qv_x, qv_y, ql_x, ql_y
 156    REAL                         :: u_x, u_y, v_x, v_y
 157    REAL, DIMENSION(kms:kme)     :: th_adv_tend, qv_adv_tend, ql_adv_tend
 158    REAL, DIMENSION(kms:kme)     :: u_adv_tend, v_adv_tend
 159    REAL, DIMENSION(kms:kme)     :: dthdz, dudz, dvdz, dqvdz, dqldz
 160    REAL                         :: w
 161    REAL, DIMENSION(kms:kme)     :: w_dthdz, w_dudz, w_dvdz, w_dqvdz, w_dqldz
 162    REAL, DIMENSION(kms:kme)     :: adv_timescale_x, adv_timescale_y
 163    CHARACTER*256                :: message
 164 ! Large-scale forcing WA 1/8/10
 165    REAL                         :: t_ls, qv_ls, ql_ls
 166    REAL                         :: u_ls, v_ls
 167    REAL, DIMENSION(kms:kme)     :: th_ls_tend, qv_ls_tend, ql_ls_tend
 168    REAL, DIMENSION(kms:kme)     :: u_ls_tend, v_ls_tend
 169    REAL, DIMENSION(kms:kme)     :: ls_timescale
 170 ! Soil forcing WA 1/3/10
 171    INTEGER                      :: ks
 172    REAL                         :: t_soil, q_soil
 173    REAL, DIMENSION(num_soil_layers) :: t_soil_tend, q_soil_tend
 174    REAL, DIMENSION(num_soil_layers) :: timescale_soil
 175
 176    IF ( scm_force .EQ. 0 ) return
 177
 178 ! NOTES
 179 ! z is kts:kte
 180 ! z_at_w is kms:kme
 181
 182      ! this is a good place for checks on the configuration
 183      if ( z_force(1) > z(ids,1,jds) ) then
 184         CALL wrf_message("First forcing level must be lower than first WRF half-level")
 185         WRITE( message , * ) 'z forcing = ',z_force(1), 'z = ',z(ids,1,jds)
 186 !       print*,"z forcing = ",z_force(1), "z = ",z(ids,1,jds)
 187         CALL wrf_error_fatal( message )
 188      endif
 189
 190      z_force = z_force + dt*z_force_tend
 191      u_g = u_g + dt*u_g_tend
 192      v_g = v_g + dt*v_g_tend
 193      tau_x = tau_x + dt*tau_x_tend
 194      tau_y = tau_y + dt*tau_y_tend
 195      tau_largescale = tau_largescale + dt*tau_largescale_tend
 196
 197      if ( scm_th_adv ) then
 198        th_upstream_x = th_upstream_x + dt*th_upstream_x_tend
 199        th_upstream_y = th_upstream_y + dt*th_upstream_y_tend
 200      endif
 201      if ( scm_qv_adv) then
 202        qv_upstream_x = qv_upstream_x + dt*qv_upstream_x_tend
 203        qv_upstream_y = qv_upstream_y + dt*qv_upstream_y_tend
 204      endif
 205      if ( scm_ql_adv) then
 206        ql_upstream_x = ql_upstream_x + dt*ql_upstream_x_tend
 207        ql_upstream_y = ql_upstream_y + dt*ql_upstream_y_tend
 208      endif
 209      if ( scm_wind_adv ) then
 210        u_upstream_x = u_upstream_x + dt*u_upstream_x_tend
 211        u_upstream_y = u_upstream_y + dt*u_upstream_y_tend
 212        v_upstream_x = v_upstream_x + dt*v_upstream_x_tend
 213        v_upstream_y = v_upstream_y + dt*v_upstream_y_tend
 214      endif
 215      if ( scm_vert_adv ) then
 216        w_subs = w_subs + dt*w_subs_tend
 217      endif
 218
 219      if ( scm_force_th_largescale ) then
 220        th_largescale = th_largescale + dt*th_largescale_tend
 221      endif
 222      if ( scm_force_qv_largescale) then
 223        qv_largescale = qv_largescale + dt*qv_largescale_tend
 224      endif
 225      if ( scm_force_ql_largescale) then
 226        ql_largescale = ql_largescale + dt*ql_largescale_tend
 227      endif
 228      if ( scm_force_wind_largescale ) then
 229        u_largescale = u_largescale + dt*u_largescale_tend
 230        v_largescale = v_largescale + dt*v_largescale_tend
 231      endif
 232
 233      if ( scm_soilT_force ) then
 234        t_soil_forcing_val = t_soil_forcing_val + dt*t_soil_forcing_tend
 235      endif
 236      if ( scm_soilQ_force ) then
 237        q_soil_forcing_val = q_soil_forcing_val + dt*q_soil_forcing_tend
 238      endif
 239
 240 ! 0 everything in case we don't set it later
 241      th_adv_tend = 0.0
 242      qv_adv_tend = 0.0
 243      ql_adv_tend = 0.0
 244      u_adv_tend  = 0.0
 245      v_adv_tend  = 0.0
 246      th_ls_tend = 0.0
 247      qv_ls_tend = 0.0
 248      ql_ls_tend = 0.0
 249      u_ls_tend  = 0.0
 250      v_ls_tend  = 0.0
 251      w_dthdz     = 0.0
 252      w_dqvdz     = 0.0
 253      w_dudz      = 0.0
 254      w_dvdz      = 0.0
 255      adv_timescale_x = 0.0
 256      adv_timescale_y = 0.0
 257
 258 ! now interpolate forcing to model vertical grid
 259
 260 !    if ( debug ) print*,' z u_base v_base '
 261      CALL wrf_debug(100,'k z_base  u_base  v_base')
 262      do k = kms,kme-1
 263        z_base(k) = z(ids,k,jds)
 264        u_base(k) = interp_0(u_g,z_force,z_base(k),num_force_layers)
 265        v_base(k) = interp_0(v_g,z_force,z_base(k),num_force_layers)
 266 !      if ( debug ) print*,z_base(k),u_base(k),v_base(k)
 267        WRITE( message, '(i4,3f12.4)' ) k,z_base(k),u_base(k),v_base(k)
 268        CALL wrf_debug ( 100, message )
 269     enddo
 270
 271     if ( scm_th_adv .or. scm_qv_adv .or. scm_wind_adv ) then
 272        if ( scm_th_adv ) CALL wrf_debug ( 100, 'k  tau_x  tau_y t_ups_x t_ups_y t_m ' )
 273        do k = kms,kme-1
 274
 275           u_x = interp_0(u_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 276           u_y = interp_0(u_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 277
 278           v_x = interp_0(v_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 279           v_y = interp_0(v_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 280
 281           adv_timescale_x(k) = interp_0(tau_x,z_force,z(ids,k,jds),num_force_layers)
 282           adv_timescale_y(k) = interp_0(tau_y,z_force,z(ids,k,jds),num_force_layers)
 283
 284
 285 !          tau_u(k) = dx/abs(u(ids,k,jds))
 286 !          tau_v(k) = dx/abs(v(ids,k,jds))
 287 !          adv_timescale_x(k) = dx/abs(u(ids,k,jds))
 288 !          adv_timescale_y(k) = dx/abs(v(ids,k,jds))
 289
 290           if ( scm_wind_adv ) then
 291              u_adv_tend(k) = (u_x-u(ids,k,jds))/adv_timescale_x(k) + (u_y-u(ids,k,jds))/adv_timescale_y(k)
 292              v_adv_tend(k) = (v_x-v(ids,k,jds))/adv_timescale_x(k) + (v_y-v(ids,k,jds))/adv_timescale_y(k)
 293           endif
 294
 295        enddo
 296     endif
 297
 298
 299     if ( scm_th_adv ) then
 300        if ( th_upstream_x(kms) > 0.) then
 301           do k = kms,kme-1
 302              t_x = interp_0(th_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 303              t_y = interp_0(th_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 304
 305              th_adv_tend(k) = (t_x-th(ids,k,jds))/adv_timescale_x(k) + (t_y-th(ids,k,jds))/adv_timescale_y(k)
 306              WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
 307              CALL wrf_debug ( 100, message )
 308           enddo
 309        else ! WA if upstream is empty, use tendency only not value+tend
 310           do k = kms,kme-1
 311              t_x = interp_0(dt*th_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 312              t_y = interp_0(dt*th_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 313
 314              th_adv_tend(k) = t_x/adv_timescale_x(k) + t_y/adv_timescale_y(k)
 315              WRITE( message, '(i4,5f12.4)' ) k,adv_timescale_x(k), adv_timescale_y(k), t_x, t_y, th(ids,k,jds)
 316              CALL wrf_debug ( 100, message )
 317           enddo
 318        endif
 319     endif
 320      if (minval(tau_x) < 0) then
 321        print*,tau_x
 322        stop 'TAU_X'
 323      endif
 324      if (minval(tau_y) < 0) then
 325        print*,z_force
 326        print*,tau_y
 327        stop 'TAU_Y'
 328      endif
 329
 330     if ( scm_qv_adv ) then
 331        if ( qv_upstream_x(kms) > 0.) then
 332           do k = kms,kme-1
 333              qv_x = interp_0(qv_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 334              qv_y = interp_0(qv_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 335
 336              qv_adv_tend(k) = (qv_x-qv(ids,k,jds))/adv_timescale_x(k) + (qv_y-qv(ids,k,jds))/adv_timescale_y(k)
 337           enddo
 338        else ! WA if upstream is empty, use tendency only not value+tend
 339           do k = kms,kme-1
 340              qv_x = interp_0(dt*qv_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 341              qv_y = interp_0(dt*qv_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 342
 343              qv_adv_tend(k) = qv_x/adv_timescale_x(k) + qv_y/adv_timescale_y(k)
 344           enddo
 345        endif
 346     endif
 347
 348     if ( scm_ql_adv ) then
 349        if ( ql_upstream_x(kms) > 0.) then
 350           do k = kms,kme-1
 351              ql_x = interp_0(ql_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 352              ql_y = interp_0(ql_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 353
 354              ql_adv_tend(k) = (ql_x-ql(ids,k,jds))/adv_timescale_x(k) + (ql_y-ql(ids,k,jds))/adv_timescale_y(k)
 355           enddo
 356        else ! WA if upstream is empty, use tendency only not value+tend
 357           do k = kms,kme-1
 358              ql_x = interp_0(dt*ql_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 359              ql_y = interp_0(dt*ql_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 360
 361              ql_adv_tend(k) = ql_x/adv_timescale_x(k) + ql_y/adv_timescale_y(k)
 362           enddo
 363        endif
 364     endif
 365
 366     if ( scm_wind_adv ) then
 367        if ( u_upstream_x(kms) > -900.) then
 368           do k = kms,kme-1
 369              u_x = interp_0(u_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 370              u_y = interp_0(u_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 371
 372              v_x = interp_0(v_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 373              v_y = interp_0(v_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 374
 375              u_adv_tend(k) = (u_x-u(ids,k,jds))/adv_timescale_x(k) + (u_y-u(ids,k,jds))/adv_timescale_y(k)
 376              v_adv_tend(k) = (v_x-v(ids,k,jds))/adv_timescale_x(k) + (v_y-v(ids,k,jds))/adv_timescale_y(k)
 377           enddo
 378        else ! WA if upstream is empty, use tendency only not value+tend
 379           do k = kms,kme-1
 380              u_x = interp_0(dt*u_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 381              u_y = interp_0(dt*u_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 382
 383              v_x = interp_0(dt*v_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 384              v_y = interp_0(dt*v_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 385
 386              u_adv_tend(k) = u_x/adv_timescale_x(k) + u_y/adv_timescale_y(k)
 387              v_adv_tend(k) = v_x/adv_timescale_x(k) + v_y/adv_timescale_y(k)
 388           enddo
 389        endif
 390     endif
 391
 392 ! Large scale forcing starts here 1/8/10 WA
 393
 394     if ( scm_force_th_largescale .or. scm_force_qv_largescale .or. scm_force_ql_largescale .or. scm_force_wind_largescale ) then
 395        do k = kms,kme-1
 396           ls_timescale(k) = interp_0(tau_largescale,z_force,z(ids,k,jds),num_force_layers)
 397        enddo
 398     endif
 399
 400     if ( scm_force_th_largescale ) then
 401        if ( th_largescale(kms) > 0.) then
 402           do k = kms,kme-1
 403              t_ls = interp_0(th_largescale,z_force,z(ids,k,jds),num_force_layers)
 404              th_ls_tend(k) = (t_ls-th(ids,k,jds))/ls_timescale(k)
 405           enddo
 406        else ! WA if upstream is empty, use tendency only not value+tend
 407           do k = kms,kme-1
 408              t_ls = interp_0(dt*th_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 409              th_ls_tend(k) = t_ls/ls_timescale(k)
 410           enddo
 411        endif
 412     endif
 413
 414     if ( scm_force_qv_largescale ) then
 415        if ( qv_largescale(kms) > 0.) then
 416           do k = kms,kme-1
 417              qv_ls = interp_0(qv_largescale,z_force,z(ids,k,jds),num_force_layers)
 418              qv_ls_tend(k) = (qv_ls-qv(ids,k,jds))/ls_timescale(k)
 419           enddo
 420        else ! WA if upstream is empty, use tendency only not value+tend
 421           do k = kms,kme-1
 422              qv_ls = interp_0(dt*qv_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 423              qv_ls_tend(k) = qv_ls/ls_timescale(k)
 424           enddo
 425        endif
 426     endif
 427
 428     if ( scm_force_ql_largescale ) then
 429        if ( ql_largescale(kms) > 0.) then
 430           do k = kms,kme-1
 431              ql_ls = interp_0(ql_largescale,z_force,z(ids,k,jds),num_force_layers)
 432              ql_ls_tend(k) = (ql_ls-ql(ids,k,jds))/ls_timescale(k)
 433           enddo
 434        else ! WA if upstream is empty, use tendency only not value+tend
 435           do k = kms,kme-1
 436              ql_ls = interp_0(dt*ql_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 437              ql_ls_tend(k) = ql_ls/ls_timescale(k)
 438           enddo
 439        endif
 440     endif
 441
 442     if ( scm_force_wind_largescale ) then
 443        if ( u_largescale(kms) > -900.) then
 444           do k = kms,kme-1
 445              u_ls = interp_0(u_largescale,z_force,z(ids,k,jds),num_force_layers)
 446              v_ls = interp_0(v_largescale,z_force,z(ids,k,jds),num_force_layers)
 447              u_ls_tend(k) = (u_ls-u(ids,k,jds))/ls_timescale(k)
 448              v_ls_tend(k) = (v_ls-v(ids,k,jds))/ls_timescale(k)
 449           enddo
 450        else ! WA if upstream is empty, use tendency only not value+tend
 451           do k = kms,kme-1
 452              u_ls = interp_0(dt*u_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 453              v_ls = interp_0(dt*v_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 454              u_ls_tend(k) = u_ls/ls_timescale(k)
 455              v_ls_tend(k) = v_ls/ls_timescale(k)
 456           enddo
 457        endif
 458     endif
 459
 460 ! Now do vertical advection.  Note that no large-scale vertical advection
 461 ! is implemented at this time, may not make sense anyway (WA).
 462 ! loops are set so that the top and bottom (w=0) are handled correctly
 463 ! vertical derivatives
 464     do k = kms+1,kme-1
 465        dthdz(k) = (th(2,k,2)-th(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 466        dqvdz(k) = (qv(2,k,2)-qv(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 467        dqldz(k) = (ql(2,k,2)-ql(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 468        dudz(k)  = (u(2,k,2)-u(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 469        dvdz(k)  = (v(2,k,2)-v(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 470     enddo
 471
 472 ! w on full levels, then advect
 473     if ( scm_vert_adv ) then
 474        do k = kms+1,kme-1
 475           w = interp_0(w_subs,z_force,z_at_w(ids,k,jds),num_force_layers)
 476           w_dthdz(k) = w*dthdz(k)
 477           w_dqvdz(k) = w*dqvdz(k)
 478           w_dqldz(k) = w*dqldz(k)
 479           w_dudz(k)  = w*dudz(k)
 480           w_dvdz(k)  = w*dvdz(k)
 481        enddo
 482     endif
 483
 484 ! set tendencies for return
 485 ! vertical advection tendencies need to be interpolated back to half levels
 486     CALL wrf_debug ( 100, 'j, k, th_adv_ten, qv_adv_ten, u_adv_ten, v_adv_ten')
 487     do j = jms,jme
 488     do k = kms,kme-1
 489     if(j==1) WRITE( message, '(i4,4f12.4)' ) k,th_adv_tend(k),qv_adv_tend(k),u_adv_tend(k),v_adv_tend(k)
 490     if(j==1) CALL wrf_debug ( 100, message )
 491     do i = ims,ime
 492        thten(i,k,j) = thten(i,k,j) + th_adv_tend(k) +              &
 493                       0.5*(w_dthdz(k) + w_dthdz(k+1))              &
 494                       + th_ls_tend(k)
 495        qvten(i,k,j) = qvten(i,k,j) + qv_adv_tend(k) +              &
 496                       0.5*(w_dqvdz(k) + w_dqvdz(k+1))              &
 497                       + qv_ls_tend(k)
 498        qlten(i,k,j) = qlten(i,k,j) + ql_adv_tend(k) +              &
 499                       0.5*(w_dqldz(k) + w_dqldz(k+1))              &
 500                       + ql_ls_tend(k)
 501        uten(i,k,j)  = uten(i,k,j) + u_adv_tend(k) +                &
 502                       0.5*(w_dudz(k) + w_dudz(k+1))                &
 503                       + u_ls_tend(k)
 504        vten(i,k,j)  = vten(i,k,j) + v_adv_tend(k) +                &
 505                       0.5*(w_dvdz(k) + w_dvdz(k+1))                &
 506                       + v_ls_tend(k)
 507     enddo
 508     enddo
 509     enddo
 510
 511 ! soil forcing 1/3/10 WA
 512     if ( scm_soilT_force ) then
 513        do ks = 1,num_soil_layers
 514           t_soil = interp_0(t_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
 515           timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
 516           t_soil_tend(ks) = (t_soil-tslb(ids,ks,jds))/timescale_soil(ks)
 517        enddo
 518        do j = jms,jme
 519           do ks = 1,num_soil_layers
 520              do i = ims,ime
 521                 tslb(ids,ks,jds) = tslb(ids,ks,jds) + t_soil_tend(ks)
 522              enddo
 523           enddo
 524        enddo
 525     endif
 526     if ( scm_soilQ_force ) then
 527        do ks = 1,num_soil_layers
 528           q_soil = interp_0(q_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
 529           timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
 530           q_soil_tend(ks) = (q_soil-smois(ids,ks,jds))/timescale_soil(ks)
 531        enddo
 532        do j = jms,jme
 533           do ks = 1,num_soil_layers
 534              do i = ims,ime
 535                 smois(ids,ks,jds) = smois(ids,ks,jds) + q_soil_tend(ks)
 536              enddo
 537           enddo
 538        enddo
 539     endif
 540
 541     RETURN
 542
 543    END SUBROUTINE force_scm
 544
 545 END MODULE module_force_scm