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