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), &
 328                                   qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
 329              CALL wrf_debug ( 100, message )
 330           enddo
 331        else ! WA if upstream is empty, use tendency only not value+tend
 332           do k = kms,kme-1
 333              qv_x = interp_0(dt*qv_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 334              qv_y = interp_0(dt*qv_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 335
 336              qv_adv_tend(k) = qv_x/adv_timescale_x(k) + qv_y/adv_timescale_y(k)
 337              WRITE( message, * ) 'qv_adv_tend branch 2',k,adv_timescale_x(k), adv_timescale_y(k), qv_upstream_x(k), &
 338                                   qv_x, qv_y, qv(ids,k,jds), qv_adv_tend(k)
 339              CALL wrf_debug ( 100, message )
 340           enddo
 341        endif
 342     endif
 343
 344     if ( scm_ql_adv ) then
 345        if ( ql_upstream_x(1) > 0.) then
 346           do k = kms,kme-1
 347              ql_x = interp_0(ql_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 348              ql_y = interp_0(ql_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 349
 350              ql_adv_tend(k) = (ql_x-ql(ids,k,jds))/adv_timescale_x(k) + (ql_y-ql(ids,k,jds))/adv_timescale_y(k)
 351           enddo
 352        else ! WA if upstream is empty, use tendency only not value+tend
 353           do k = kms,kme-1
 354              ql_x = interp_0(dt*ql_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 355              ql_y = interp_0(dt*ql_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 356
 357              ql_adv_tend(k) = ql_x/adv_timescale_x(k) + ql_y/adv_timescale_y(k)
 358           enddo
 359        endif
 360     endif
 361
 362     if ( scm_wind_adv ) then
 363        if ( u_upstream_x(1) > -900.) then
 364           do k = kms,kme-1
 365              u_x = interp_0(u_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 366              u_y = interp_0(u_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 367
 368              v_x = interp_0(v_upstream_x,z_force,z(ids,k,jds),num_force_layers)
 369              v_y = interp_0(v_upstream_y,z_force,z(ids,k,jds),num_force_layers)
 370
 371              u_adv_tend(k) = (u_x-u(ids,k,jds))/adv_timescale_x(k) + (u_y-u(ids,k,jds))/adv_timescale_y(k)
 372              v_adv_tend(k) = (v_x-v(ids,k,jds))/adv_timescale_x(k) + (v_y-v(ids,k,jds))/adv_timescale_y(k)
 373           enddo
 374        else ! WA if upstream is empty, use tendency only not value+tend
 375           do k = kms,kme-1
 376              u_x = interp_0(dt*u_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 377              u_y = interp_0(dt*u_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 378
 379              v_x = interp_0(dt*v_upstream_x_tend,z_force,z(ids,k,jds),num_force_layers)
 380              v_y = interp_0(dt*v_upstream_y_tend,z_force,z(ids,k,jds),num_force_layers)
 381
 382              u_adv_tend(k) = u_x/adv_timescale_x(k) + u_y/adv_timescale_y(k)
 383              v_adv_tend(k) = v_x/adv_timescale_x(k) + v_y/adv_timescale_y(k)
 384           enddo
 385        endif
 386     endif
 387
 388
 389
 390     if ( scm_th_t_tend ) then
 391        do k = kms,kme-1
 392           th_t_tend_interp(k) = interp_0(th_t_tend,z_force,z(ids,k,jds),num_force_layers)
 393        enddo
 394     endif
 395
 396     if ( scm_qv_t_tend ) then
 397        do k = kms,kme-1
 398           qv_t_tend_interp(k) = interp_0(qv_t_tend,z_force,z(ids,k,jds),num_force_layers)
 399           write(*,'(i3, f20.15)') k, qv_t_tend_interp(k)
 400        enddo
 401     endif
 402
 403
 404 ! Large scale forcing starts here 1/8/10 WA
 405     if ( scm_force_th_largescale .or. scm_force_qv_largescale .or. scm_force_ql_largescale .or. scm_force_wind_largescale ) then
 406        do k = kms,kme-1
 407           ls_timescale(k) = interp_0(tau_largescale,z_force,z(ids,k,jds),num_force_layers)
 408        enddo
 409     endif
 410
 411     if ( scm_force_th_largescale ) then
 412        if ( th_largescale(1) > 0.) then
 413           do k = kms,kme-1
 414              t_ls = interp_0(th_largescale,z_force,z(ids,k,jds),num_force_layers)
 415              th_ls_tend(k) = (t_ls-th(ids,k,jds))/ls_timescale(k)
 416           enddo
 417        else ! WA if upstream is empty, use tendency only not value+tend
 418           do k = kms,kme-1
 419              t_ls = interp_0(dt*th_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 420              th_ls_tend(k) = t_ls/ls_timescale(k)
 421           enddo
 422        endif
 423     endif
 424
 425     if ( scm_force_qv_largescale ) then
 426        if ( qv_largescale(1) > 0.) then
 427           do k = kms,kme-1
 428              qv_ls = interp_0(qv_largescale,z_force,z(ids,k,jds),num_force_layers)
 429              qv_ls_tend(k) = (qv_ls-qv(ids,k,jds))/ls_timescale(k)
 430           enddo
 431        else ! WA if upstream is empty, use tendency only not value+tend
 432           do k = kms,kme-1
 433              qv_ls = interp_0(dt*qv_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 434              qv_ls_tend(k) = qv_ls/ls_timescale(k)
 435           enddo
 436        endif
 437     endif
 438
 439     if ( scm_force_ql_largescale ) then
 440        if ( ql_largescale(1) > 0.) then
 441           do k = kms,kme-1
 442              ql_ls = interp_0(ql_largescale,z_force,z(ids,k,jds),num_force_layers)
 443              ql_ls_tend(k) = (ql_ls-ql(ids,k,jds))/ls_timescale(k)
 444           enddo
 445        else ! WA if upstream is empty, use tendency only not value+tend
 446           do k = kms,kme-1
 447              ql_ls = interp_0(dt*ql_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 448              ql_ls_tend(k) = ql_ls/ls_timescale(k)
 449           enddo
 450        endif
 451     endif
 452
 453     if ( scm_force_wind_largescale ) then
 454        if ( u_largescale(1) > -900.) then
 455           do k = kms,kme-1
 456              u_ls = interp_0(u_largescale,z_force,z(ids,k,jds),num_force_layers)
 457              v_ls = interp_0(v_largescale,z_force,z(ids,k,jds),num_force_layers)
 458              u_ls_tend(k) = (u_ls-u(ids,k,jds))/ls_timescale(k)
 459              v_ls_tend(k) = (v_ls-v(ids,k,jds))/ls_timescale(k)
 460           enddo
 461        else ! WA if upstream is empty, use tendency only not value+tend
 462           do k = kms,kme-1
 463              u_ls = interp_0(dt*u_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 464              v_ls = interp_0(dt*v_largescale_tend,z_force,z(ids,k,jds),num_force_layers)
 465              u_ls_tend(k) = u_ls/ls_timescale(k)
 466              v_ls_tend(k) = v_ls/ls_timescale(k)
 467           enddo
 468        endif
 469     endif
 470
 471 ! Now do vertical advection.  Note that no large-scale vertical advection
 472 ! is implemented at this time, may not make sense anyway (WA).
 473 ! loops are set so that the top and bottom (w=0) are handled correctly
 474 ! vertical derivatives
 475     do k = kms+1,kme-1
 476        dthdz(k) = (th(2,k,2)-th(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 477        dqvdz(k) = (qv(2,k,2)-qv(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 478        dqldz(k) = (ql(2,k,2)-ql(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 479        dudz(k)  = (u(2,k,2)-u(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 480        dvdz(k)  = (v(2,k,2)-v(2,k-1,2))/(z(2,k,2)-z(2,k-1,2))
 481     enddo
 482
 483 ! w on full levels, then advect
 484     if ( scm_vert_adv ) then
 485        do k = kms+1,kme-1
 486           w = interp_0(w_subs,z_force,z_at_w(ids,k,jds),num_force_layers)
 487           w_dthdz(k) = -w*dthdz(k)
 488           w_dqvdz(k) = -w*dqvdz(k)
 489           w_dqldz(k) = -w*dqldz(k)
 490           w_dudz(k)  = -w*dudz(k)
 491           w_dvdz(k)  = -w*dvdz(k)
 492        enddo
 493     endif
 494
 495 ! set tendencies for return
 496 ! vertical advection tendencies need to be interpolated back to half levels
 497     CALL wrf_debug ( 100, 'j, k, th_adv_ten, qv_adv_ten, ql_adv_ten, u_adv_ten, v_adv_ten')
 498     do j = jms,jme
 499     do k = kms,kme-1
 500     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)
 501     if(j==1) CALL wrf_debug ( 100, message )
 502     do i = ims,ime
 503        thten(i,k,j) = thten(i,k,j) + th_adv_tend(k) +              &
 504                       0.5*(w_dthdz(k) + w_dthdz(k+1)) + th_t_tend_interp(k)&
 505                       + th_ls_tend(k)
 506        qvten(i,k,j) = qvten(i,k,j) + qv_adv_tend(k) +              &
 507                       0.5*(w_dqvdz(k) + w_dqvdz(k+1)) + qv_t_tend_interp(k)&
 508                       + qv_ls_tend(k)
 509        qlten(i,k,j) = qlten(i,k,j) + ql_adv_tend(k) +              &
 510                       0.5*(w_dqldz(k) + w_dqldz(k+1))              &
 511                       + ql_ls_tend(k)
 512        uten(i,k,j)  = uten(i,k,j) + u_adv_tend(k) +                &
 513                       0.5*(w_dudz(k) + w_dudz(k+1))                &
 514                       + u_ls_tend(k)
 515        vten(i,k,j)  = vten(i,k,j) + v_adv_tend(k) +                &
 516                       0.5*(w_dvdz(k) + w_dvdz(k+1))                &
 517                       + v_ls_tend(k)
 518     enddo
 519     enddo
 520     enddo
 521
 522 ! soil forcing 1/3/10 WA
 523     if ( scm_soilT_force ) then
 524        do ks = 1,num_soil_layers
 525           t_soil = interp_0(t_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
 526           timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
 527           t_soil_tend(ks) = (t_soil-tslb(ids,ks,jds))/timescale_soil(ks)
 528        enddo
 529        do j = jms,jme
 530           do ks = 1,num_soil_layers
 531              do i = ims,ime
 532                 tslb(ids,ks,jds) = tslb(ids,ks,jds) + t_soil_tend(ks)
 533              enddo
 534           enddo
 535        enddo
 536     endif
 537     if ( scm_soilQ_force ) then
 538        do ks = 1,num_soil_layers
 539           q_soil = interp_0(q_soil_forcing_val,soil_depth_force,zs(ks),num_force_soil_layers)
 540           timescale_soil(ks) = interp_0(tau_soil,soil_depth_force,zs(ks),num_force_soil_layers)
 541           q_soil_tend(ks) = (q_soil-smois(ids,ks,jds))/timescale_soil(ks)
 542        enddo
 543        do j = jms,jme
 544           do ks = 1,num_soil_layers
 545              do i = ims,ime
 546                 smois(ids,ks,jds) = smois(ids,ks,jds) + q_soil_tend(ks)
 547              enddo
 548           enddo
 549        enddo
 550     endif
 551
 552     RETURN
 553
 554    END SUBROUTINE force_scm
 555
 556 END MODULE module_force_scm