Index: palm/trunk/SOURCE/ls_forcing.f90
===================================================================
--- palm/trunk/SOURCE/ls_forcing.f90	(revision 1354)
+++ palm/trunk/SOURCE/ls_forcing.f90	(revision 1365)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! ------------------
-! 
+! Usage of large scale forcing for pt and q enabled:
+! Added new subroutine ls_advec for horizontal large scale advection and large 
+! scale subsidence, 
+! error message in init_ls_forcing specified,
+! variable t renamed nt
 ! 
 ! Former revisions:
@@ -61,7 +65,11 @@
 
     PRIVATE
-    PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert 
+    PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec
     SAVE
 
+    INTERFACE ls_advec
+       MODULE PROCEDURE ls_advec
+       MODULE PROCEDURE ls_advec_ij
+    END INTERFACE ls_advec
 
  CONTAINS
@@ -70,6 +78,7 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  p_surf, pt_surf, q_surf, qsws_surf, shf_surf, time_surf,     &
-                  time_vert, ug_vert, vg_vert, wsubs_vert, zu
+           ONLY:  p_surf, pt_lsa, pt_subs, pt_surf, q_lsa, q_subs, q_surf,     &
+                  qsws_surf, shf_surf, time_surf, time_vert, ug_vert, vg_vert, &
+                  wsubs_vert, zu
 
        USE control_parameters,                                                 &
@@ -77,40 +86,61 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzt
+           ONLY:  ngp_sums_ls, nzb, nz, nzt
 
        USE kinds
 
+       USE statistics,                                                         &
+           ONLY:  sums_ls_l
+
 
        IMPLICIT NONE
 
-       CHARACTER(100) ::  chmess              !:
-       CHARACTER(1)   ::  hash                !:
-
-       INTEGER(iwp) ::  ierrn                 !:
-       INTEGER(iwp) ::  finput = 90           !:
-       INTEGER(iwp) ::  k                     !:
-       INTEGER(iwp) ::  t                     !:
-
-       REAL(wp) ::  fac                       !:
-       REAL(wp) ::  highheight                !:
-       REAL(wp) ::  highug_vert               !:
-       REAL(wp) ::  highvg_vert               !:
-       REAL(wp) ::  highwsubs_vert            !:
-       REAL(wp) ::  lowheight                 !:
-       REAL(wp) ::  lowug_vert                !:
-       REAL(wp) ::  lowvg_vert                !:
-       REAL(wp) ::  lowwsubs_vert             !:
-       REAL(wp) ::  r_dummy                   !:
-
-       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),         &
-                 qsws_surf(0:nlsf), shf_surf(0:nlsf), time_vert(0:nlsf),  &
-                 time_surf(0:nlsf), ug_vert(nzb:nzt+1,0:nlsf),            &
-                 vg_vert(nzb:nzt+1,0:nlsf), wsubs_vert(nzb:nzt+1,0:nlsf) )
-
-       p_surf     = 0.0_wp;  pt_surf  = 0.0_wp;  q_surf    = 0.0_wp
-       qsws_surf  = 0.0_wp;  shf_surf = 0.0_wp;  time_vert = 0.0_wp
-       time_surf  = 0.0_wp;  ug_vert  = 0.0_wp;  vg_vert   = 0.0_wp
-       wsubs_vert = 0.0_wp
-
+       CHARACTER(100) ::  chmess      !:
+       CHARACTER(1)   ::  hash        !:
+
+       INTEGER(iwp) ::  ierrn         !:
+       INTEGER(iwp) ::  finput = 90   !:
+       INTEGER(iwp) ::  k             !:
+       INTEGER(iwp) ::  nt             !:
+
+       REAL(wp) ::  fac               !:
+       REAL(wp) ::  highheight        !:
+       REAL(wp) ::  highug_vert       !:
+       REAL(wp) ::  highvg_vert       !:
+       REAL(wp) ::  highwsubs_vert    !:
+       REAL(wp) ::  lowheight         !:
+       REAL(wp) ::  lowug_vert        !:
+       REAL(wp) ::  lowvg_vert        !:
+       REAL(wp) ::  lowwsubs_vert     !:
+       REAL(wp) ::  highpt_lsa        !:
+       REAL(wp) ::  lowpt_lsa         !:
+       REAL(wp) ::  highq_lsa         !:
+       REAL(wp) ::  lowq_lsa          !:
+       REAL(wp) ::  highpt_subs       !:
+       REAL(wp) ::  lowpt_subs        !:
+       REAL(wp) ::  highq_subs        !:
+       REAL(wp) ::  lowq_subs         !:
+       REAL(wp) ::  r_dummy           !:
+
+       ALLOCATE( p_surf(0:nlsf), pt_lsa(nzb:nzt+1,0:nlsf),                     &
+                 pt_subs(nzb:nzt+1,0:nlsf), pt_surf(0:nlsf),                   &
+                 q_lsa(nzb:nzt+1,0:nlsf), q_subs(nzb:nzt+1,0:nlsf),            &
+                 q_surf(0:nlsf), qsws_surf(0:nlsf), shf_surf(0:nlsf),          &
+                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
+                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
+                 wsubs_vert(nzb:nzt+1,0:nlsf) )
+
+       p_surf = 0.0_wp; pt_lsa = 0.0; pt_subs = 0.0; pt_surf = 0.0_wp
+       q_lsa = 0.0; q_subs = 0.0; q_surf = 0.0_wp; qsws_surf = 0.0_wp
+       shf_surf = 0.0_wp; time_vert = 0.0_wp; time_surf = 0.0_wp;
+       ug_vert = 0.0_wp; vg_vert = 0.0_wp; wsubs_vert = 0.0_wp
+
+!
+!--    Array for storing large scale forcing and nudging tendencies at each 
+!--    timestep for data output
+       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
+       sums_ls_l = 0.0_wp
+
+       ngp_sums_ls = (nz+2)*6
 
        OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
@@ -136,12 +166,12 @@
 !
 !--    Surface values are read in
-       t     = 0
+       nt     = 0
        ierrn = 0
 
-       DO WHILE ( time_surf(t) < end_time )
-          t = t + 1
-          READ ( finput, *, IOSTAT = ierrn ) time_surf(t), shf_surf(t), &
-                                             qsws_surf(t), pt_surf(t),  &
-                                             q_surf(t), p_surf(t)
+       DO WHILE ( time_surf(nt) < end_time )
+          nt = nt + 1
+          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
+                                             qsws_surf(nt), pt_surf(nt),       &
+                                             q_surf(nt), p_surf(nt)
 
           IF ( ierrn < 0 )  THEN
@@ -156,5 +186,5 @@
        IF ( time_surf(1) > end_time )  THEN
           WRITE ( message_string, * ) 'No time dependent surface variables in ',&
-                                     '&LSF_DATA for end of run found - ',  &
+                                     '&LSF_DATA for end of run found - ',      &
                                      'lsf_surf is set to FALSE'
           CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
@@ -171,7 +201,7 @@
 !--    Profiles of ug, vg and w_subs are read in (large scale forcing)
 
-       t = 0
-       DO WHILE ( time_vert(t) < end_time )
-          t = t + 1
+       nt = 0
+       DO WHILE ( time_vert(nt) < end_time )
+          nt = nt + 1
           hash = "#"
           ierrn = 1 ! not zero
@@ -180,5 +210,5 @@
 !--       from there onwards the profiles will be read
           DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
-             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(t)
+             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
              IF ( ierrn < 0 )  THEN 
                 WRITE( message_string, * ) 'No time dependent vertical profiles',&
@@ -188,8 +218,9 @@
           ENDDO
 
-          IF ( t == 1 .AND. time_vert(t) > end_time ) EXIT
-
-          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, &
-                                           lowwsubs_vert
+          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
+
+          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
+                                           lowwsubs_vert, lowpt_lsa,           &
+                                           lowq_lsa, lowpt_subs, lowq_subs
           IF ( ierrn /= 0 )  THEN
              message_string = 'errors in file LSF_DATA'
@@ -197,6 +228,8 @@
           ENDIF
 
-          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, highvg_vert, &
-                                           highwsubs_vert
+          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
+                                           highvg_vert, highwsubs_vert,        &
+                                           highpt_lsa, highq_lsa,              &
+                                           highpt_subs, highq_subs
       
           IF ( ierrn /= 0 )  THEN
@@ -212,12 +245,21 @@
                 lowvg_vert    = highvg_vert
                 lowwsubs_vert = highwsubs_vert
+                lowpt_lsa     = highpt_lsa
+                lowq_lsa      = highq_lsa
+                lowpt_subs    = highpt_subs
+                lowq_subs     = highq_subs
 
                 ierrn = 0
-                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
-                                                 highvg_vert, highwsubs_vert
+                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
+                                                 highvg_vert, highwsubs_vert,  &
+                                                 highpt_lsa, highq_lsa,        &
+                                                 highpt_subs, highq_subs
 
                 IF ( ierrn /= 0 )  THEN
-                   message_string = 'errors in file LSF_DATA'
-                   CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
+                   WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm ', &
+                        'is higher than the maximum height in LSF_DATA which ',&
+                        'is ', lowheight, 'm. Interpolation on PALM ',         &
+                        'grid is not possible.'
+                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
                 ENDIF
 
@@ -228,7 +270,12 @@
              fac = (highheight-zu(k))/(highheight - lowheight)
 
-             ug_vert(k,t)    = fac*lowug_vert    + (1-fac)*highug_vert
-             vg_vert(k,t)    = fac*lowvg_vert    + (1-fac)*highvg_vert
-             wsubs_vert(k,t) = fac*lowwsubs_vert + (1-fac)*highwsubs_vert
+             ug_vert(k,nt)    = fac*lowug_vert    + (1.0_wp-fac)*highug_vert
+             vg_vert(k,nt)    = fac*lowvg_vert    + (1.0_wp-fac)*highvg_vert
+             wsubs_vert(k,nt) = fac*lowwsubs_vert + (1.0_wp-fac)*highwsubs_vert
+
+             pt_lsa(k,nt)     = fac*lowpt_lsa     + (1.0_wp-fac)*highpt_lsa
+             q_lsa(k,nt)      = fac*lowq_lsa      + (1.0_wp-fac)*highq_lsa
+             pt_subs(k,nt)    = fac*lowpt_subs    + (1.0_wp-fac)*highpt_subs
+             q_subs(k,nt)     = fac*lowq_subs     + (1.0_wp-fac)*highq_subs
 
           ENDDO
@@ -241,6 +288,6 @@
  
        IF ( time_vert(1) > end_time )  THEN
-          WRITE ( message_string, * ) 'Time dependent large scale profile ',&
-                             'forcing from&LSF_DATA sets in after end of ' ,&
+          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
+                             'forcing from&LSF_DATA sets in after end of ' ,   &
                              'simulation - lsf_vert is set to FALSE'
           CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
@@ -266,8 +313,7 @@
        USE kinds
 
-
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  t                     !:
+       INTEGER(iwp) ::  nt                     !:
 
        REAL(wp)             :: fac            !:
@@ -276,13 +322,13 @@
 !
 !--    Interpolation in time of LSF_DATA at the surface
-       t = 1
-       DO WHILE ( time > time_surf(t) )
-          t = t + 1
+       nt = 1
+       DO WHILE ( time > time_surf(nt) )
+          nt = nt + 1
        ENDDO
-       IF ( time /= time_surf(t) )  THEN
-         t = t - 1
-       ENDIF
-
-       fac = ( time -time_surf(t) ) / ( time_surf(t+1) - time_surf(t) )
+       IF ( time /= time_surf(nt) )  THEN
+         nt = nt - 1
+       ENDIF
+
+       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
 
        IF ( ibc_pt_b == 0 )  THEN
@@ -290,5 +336,5 @@
 !--       In case of Dirichlet boundary condition shf must not 
 !--       be set - it is calculated via MOST in prandtl_fluxes
-          pt_surface = pt_surf(t) + fac * ( pt_surf(t+1) - pt_surf(t) )
+          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
 
        ELSEIF ( ibc_pt_b == 1 )  THEN
@@ -296,6 +342,6 @@
 !--       In case of Neumann boundary condition pt_surface is needed for 
 !--       calculation of reference density
-          shf        = shf_surf(t) + fac * ( shf_surf(t+1) - shf_surf(t) )
-          pt_surface = pt_surf(t) + fac * ( pt_surf(t+1) - pt_surf(t) )
+          shf        = shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) )
+          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
 
        ENDIF
@@ -305,13 +351,13 @@
 !--       In case of Dirichlet boundary condition qsws must not 
 !--       be set - it is calculated via MOST in prandtl_fluxes
-          q_surface = q_surf(t) + fac * ( q_surf(t+1) - q_surf(t) )
+          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
 
        ELSEIF ( ibc_q_b == 1 )  THEN
 
-          qsws = qsws_surf(t) + fac * ( qsws_surf(t+1) - qsws_surf(t) )
-
-       ENDIF
-
-       surface_pressure = p_surf(t) + fac * ( p_surf(t+1) - p_surf(t) )
+          qsws = qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) )
+
+       ENDIF
+
+       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
 
     END SUBROUTINE ls_forcing_surf 
@@ -328,8 +374,7 @@
        USE kinds
 
-
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  t                     !:
+       INTEGER(iwp) ::  nt                     !:
 
        REAL(wp)             ::  fac           !:
@@ -338,19 +383,20 @@
 !
 !--    Interpolation in time of LSF_DATA for ug, vg and w_subs
-       t = 1
-       DO WHILE ( time > time_vert(t) )
-          t = t + 1
+       nt = 1
+       DO WHILE ( time > time_vert(nt) )
+          nt = nt + 1
        ENDDO
-       IF ( time /= time_vert(t) )  THEN
-         t = t - 1
-       ENDIF
-
-       fac = ( time-time_vert(t) ) / ( time_vert(t+1)-time_vert(t) )
-
-       ug     = ug_vert(:,t) + fac * ( ug_vert(:,t+1) - ug_vert(:,t) )
-       vg     = vg_vert(:,t) + fac * ( vg_vert(:,t+1) - vg_vert(:,t) )
+       IF ( time /= time_vert(nt) )  THEN
+         nt = nt - 1
+       ENDIF
+
+       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
+
+       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
+       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
 
        IF ( large_scale_subsidence )  THEN
-          w_subs = wsubs_vert(:,t) + fac * ( wsubs_vert(:,t+1) - wsubs_vert(:,t) )
+          w_subs = wsubs_vert(:,nt)                                            &
+                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
        ENDIF
 
@@ -358,3 +404,191 @@
 
 
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE ls_advec ( time, prog_var )
+
+       USE arrays_3d,                                                          &
+           ONLY:  pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert        
+       
+       USE control_parameters,                                                 &
+           ONLY:  large_scale_subsidence, use_subsidence_tendencies
+       
+       USE indices
+       
+       USE kinds
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=*) ::  prog_var   !: 
+
+       REAL(wp), INTENT(in)  :: time    !: 
+       REAL(wp) :: fac                  !:  
+
+       INTEGER(iwp) ::  i               !: 
+       INTEGER(iwp) ::  j               !: 
+       INTEGER(iwp) ::  k               !: 
+       INTEGER(iwp) ::  nt               !: 
+
+!
+!--    Interpolation in time of LSF_DATA 
+       nt = 1
+       DO WHILE ( time > time_vert(nt) )
+          nt = nt + 1
+       ENDDO
+       IF ( time /= time_vert(nt) )  THEN
+         nt = nt - 1
+       ENDIF
+
+       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
+
+!
+!--    Add horizontal large scale advection tendencies of pt and q 
+       SELECT CASE ( prog_var )
+
+          CASE ( 'pt' )
+
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_u_inner(j,i)+1, nzt
+                      tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt)                 &
+                                    + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'q' )
+
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_u_inner(j,i)+1, nzt
+                      tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt)                  &
+                                    + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+       END SELECT
+
+!
+!--    Subsidence of pt and q with prescribed subsidence tendencies
+       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
+
+          SELECT CASE ( prog_var )
+
+             CASE ( 'pt' )
+
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_u_inner(j,i)+1, nzt
+                         tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac *     &
+                                       ( pt_subs(k,nt+1) - pt_subs(k,nt) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+  
+             CASE ( 'q' )
+
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_u_inner(j,i)+1, nzt
+                         tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac *      &
+                                       ( q_subs(k,nt+1) - q_subs(k,nt) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+
+          END SELECT
+
+       ENDIF
+
+    END SUBROUTINE ls_advec
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
+
+       USE arrays_3d,                                                          &
+           ONLY:  pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert        
+       
+       USE control_parameters,                                                 &
+           ONLY:  large_scale_subsidence, use_subsidence_tendencies
+       
+       USE indices
+       
+       USE kinds
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=*) ::  prog_var   !: 
+
+       REAL(wp), INTENT(in)  :: time    !: 
+       REAL(wp) :: fac                  !: 
+
+       INTEGER(iwp) ::  i               !: 
+       INTEGER(iwp) ::  j               !: 
+       INTEGER(iwp) ::  k               !: 
+       INTEGER(iwp) ::  nt               !: 
+
+!
+!--    Interpolation in time of LSF_DATA 
+       nt = 1
+       DO WHILE ( time > time_vert(nt) )
+          nt = nt + 1
+       ENDDO
+       IF ( time /= time_vert(nt) )  THEN
+         nt = nt - 1
+       ENDIF
+
+       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
+
+!
+!--    Add horizontal large scale advection tendencies of pt and q 
+       SELECT CASE ( prog_var )
+
+          CASE ( 'pt' )
+
+             DO  k = nzb_u_inner(j,i)+1, nzt
+                tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt)                       &
+                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
+             ENDDO
+
+          CASE ( 'q' )
+
+             DO  k = nzb_u_inner(j,i)+1, nzt
+                tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt)                        &
+                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
+             ENDDO
+
+       END SELECT
+
+!
+!--    Subsidence of pt and q with prescribed profiles
+       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
+
+          SELECT CASE ( prog_var )
+
+             CASE ( 'pt' )
+
+                DO  k = nzb_u_inner(j,i)+1, nzt
+                   tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac *           &
+                                 ( pt_subs(k,nt+1) - pt_subs(k,nt) )
+                ENDDO
+  
+             CASE ( 'q' )
+
+                DO  k = nzb_u_inner(j,i)+1, nzt
+                   tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac *            &
+                                 ( q_subs(k,nt+1) - q_subs(k,nt) )
+                ENDDO
+
+          END SELECT
+
+       ENDIF
+
+    END SUBROUTINE ls_advec_ij
+
+
  END MODULE ls_forcing_mod
