Changeset 1365 for palm/trunk/SOURCE

Timestamp:

Apr 22, 2014 3:03:56 PM (11 years ago)

Author:

boeske

Message:

large scale forcing enabled

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (5 diffs)
• buoyancy.f90 (modified) (4 diffs)
• calc_mean_profile.f90 (added)
• check_parameters.f90 (modified) (3 diffs)
• flow_statistics.f90 (modified) (18 diffs)
• header.f90 (modified) (10 diffs)
• ls_forcing.f90 (modified) (21 diffs)
• modules.f90 (modified) (7 diffs)
• nudging.f90 (modified) (22 diffs)
• parin.f90 (modified) (3 diffs)
• prognostic_equations.f90 (modified) (10 diffs)
• subsidence.f90 (modified) (21 diffs)
• time_integration.f90 (modified) (8 diffs)

palm/trunk/SOURCE/Makefile

@@ -20 +20 @@
 # Current revisions:
 # ------------------
-# 
+# Added new module calc_mean_profile, previously in module buoyancy,
+# removed buoyancy dependency from nudging
 # 
 # Former revisions:
@@ -162 +163 @@
         advec_u_pw.f90 advec_u_up.f90 advec_v_pw.f90 advec_v_up.f90 \
         advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
-        buoyancy.f90 calc_liquid_water_content.f90 calc_precipitation.f90 \
-        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
-        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
-        coriolis.f90 cpulog.f90 cuda_fft_interfaces.f90 data_log.f90 \
-        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
+        buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
+        calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
+        check_for_restart.f90 check_open.f90 check_parameters.f90 \
+        close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
+        cuda_fft_interfaces.f90 data_log.f90 data_output_dvrp.f90 \
+        data_output_mask.f90 data_output_profiles.f90 \
         data_output_ptseries.f90 data_output_spectra.f90 \
         data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \
@@ -256 +258 @@
 boundary_conds.o: modules.o mod_kinds.o
 buoyancy.o: modules.o mod_kinds.o
+calc_mean_profile.o: modules.o mod_kinds.o
 calc_liquid_water_content.o: modules.o mod_kinds.o
 calc_precipitation.o: modules.o mod_kinds.o
@@ -352 +355 @@
 mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
 netcdf.o: modules.o mod_kinds.o
-nudging.o: modules.o buoyancy.o cpulog.o mod_kinds.o
+nudging.o: modules.o cpulog.o mod_kinds.o
 package_parin.o: modules.o mod_kinds.o
 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
@@ -382 +385 @@
 swap_timelevel.o: modules.o cpulog.o mod_kinds.o
 temperton_fft.o: modules.o mod_kinds.o
-time_integration.o: modules.o advec_ws.o buoyancy.o cpulog.o interaction_droplets_ptq.o \
-        ls_forcing.o mod_kinds.o production_e.o prognostic_equations.o user_actions.o
+time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \
+        cpulog.o interaction_droplets_ptq.o ls_forcing.o mod_kinds.o \
+        production_e.o prognostic_equations.o user_actions.o
 time_to_string.o: mod_kinds.o
 timestep.o: modules.o cpulog.o mod_kinds.o

palm/trunk/SOURCE/buoyancy.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Calculation of reference state in subroutine calc_mean_profile moved to
+! subroutine time_integration,
+! subroutine calc_mean_profile moved to new file calc_mean_profile.f90
 ! 
 ! Former revisions:
@@ -78 +80 @@
 
     PRIVATE
-    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
+    PUBLIC buoyancy, buoyancy_acc
 
     INTERFACE buoyancy
@@ -88 +90 @@
        MODULE PROCEDURE buoyancy_acc
     END INTERFACE buoyancy_acc
-
-    INTERFACE calc_mean_profile
-       MODULE PROCEDURE calc_mean_profile
-    END INTERFACE calc_mean_profile
 
  CONTAINS
@@ -376 +374 @@
     END SUBROUTINE buoyancy_ij
 
-
-    SUBROUTINE calc_mean_profile( var, pr, loc )
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
-! of potential temperature, 44 in case of virtual potential temperature, and 64
-! in case of density (ocean runs)).
-!------------------------------------------------------------------------------!
-
-       USE arrays_3d,                                                          &
-           ONLY:  ref_state
-
-       USE control_parameters,                                                 &
-           ONLY:  intermediate_timestep_count, message_string
-
-       USE indices,                                                            &
-           ONLY:  ngp_2dh_s_inner, nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
-
-       USE kinds
-
-       USE pegrid
-
-       USE statistics,                                                         &
-           ONLY:  flow_statistics_called, hom, sums, sums_l
-
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN=*) ::  loc !:
-       
-       INTEGER(iwp) ::  i                  !:
-       INTEGER(iwp) ::  j                  !:
-       INTEGER(iwp) ::  k                  !:
-       INTEGER(iwp) ::  pr                 !: 
-       INTEGER(iwp) ::  omp_get_thread_num !:
-       INTEGER(iwp) ::  tn                 !:
-       
-#if defined( __nopointer )
-       REAL(wp), DIMENSION(:,:,:) ::  var  !:
-#else
-       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
-#endif
-
-!
-!--    Computation of the horizontally averaged profile of variable var, unless
-!--    already done by the relevant call from flow_statistics. The calculation
-!--    is done only for the first respective intermediate timestep in order to
-!--    spare communication time and to produce identical model results with jobs
-!--    which are calling flow_statistics at different time intervals.
-       IF ( .NOT. flow_statistics_called  .AND.                                &
-            intermediate_timestep_count == 1 )  THEN
-
-!
-!--       Horizontal average of variable var
-          tn           =   0  ! Default thread number in case of one thread
-          !$OMP PARALLEL PRIVATE( i, j, k, tn )
-!$        tn = omp_get_thread_num()
-          sums_l(:,pr,tn) = 0.0_wp
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-
-          DO  i = 1, threads_per_task-1
-             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
-          ENDDO
-
-#if defined( __parallel )
-
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb,       &
-                              MPI_REAL, MPI_SUM, comm2d, ierr )
-
-#else
-
-          sums(:,pr) = sums_l(:,pr,0)
-
-#endif
-
-          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
-
-       ENDIF
-
-       SELECT CASE ( loc )
-
-          CASE ( 'time_int' )
-
-             ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
-
-
-          CASE ( 'nudging' )
-             !nothing to be done
-
-
-          CASE DEFAULT
-             message_string = 'unknown location "' // loc // '"'
-             CALL message( 'calc_mean_profile', 'PA0379', 1, 2, 0, 6, 0 )
-
-       END SELECT
-
-
-
-    END SUBROUTINE calc_mean_profile
-
  END MODULE buoyancy_mod

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Usage of large scale forcing for pt and q enabled:
+! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q), 
+! large scale subsidence (td_sub_lpt, td_sub_q) 
+! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
+! check if use_subsidence_tendencies is used correctly
 ! 
 ! Former revisions:
@@ -1271 +1275 @@
        ENDIF
          
+    ENDIF
+
+!
+!-- Check if the control parameter use_subsidence_tendencies is used correctly
+    IF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_subsidence )  THEN
+       message_string = 'The usage of use_subsidence_tendencies ' // &
+                            'requires large_scale_subsidence = .T..'
+       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
+    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
+       message_string = 'The usage of use_subsidence_tendencies ' // &
+                            'requires large_scale_forcing = .T..'
+       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -2661 +2677 @@
              ENDIF
 
+          CASE ( 'td_lsa_lpt' )
+             IF ( .NOT. large_scale_forcing )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for large_scale_forcing = .FALSE.'
+                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 81
+                dopr_unit(i)  = 'K/s'
+                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_lsa_q' )
+             IF ( .NOT. large_scale_forcing )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for large_scale_forcing = .FALSE.'
+                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 82
+                dopr_unit(i)  = 'kg/kgs'
+                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_sub_lpt' )
+             IF ( .NOT. large_scale_forcing )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for large_scale_forcing = .FALSE.'
+                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 83
+                dopr_unit(i)  = 'K/s'
+                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_sub_q' )
+             IF ( .NOT. large_scale_forcing )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for large_scale_forcing = .FALSE.'
+                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 84
+                dopr_unit(i)  = 'kg/kgs'
+                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_nud_lpt' )
+             IF ( .NOT. nudging )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for nudging = .FALSE.'
+                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 85
+                dopr_unit(i)  = 'K/s'
+                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_nud_q' )
+             IF ( .NOT. nudging )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for nudging = .FALSE.'
+                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 86
+                dopr_unit(i)  = 'kg/kgs'
+                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_nud_u' )
+             IF ( .NOT. nudging )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for nudging = .FALSE.'
+                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 87
+                dopr_unit(i)  = 'm/s2'
+                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'td_nud_v' )
+             IF ( .NOT. nudging )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for nudging = .FALSE.'
+                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 88
+                dopr_unit(i)  = 'm/s2'
+                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+
           CASE DEFAULT
 

palm/trunk/SOURCE/flow_statistics.f90

@@ -21 +21 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of large scale advection, large scale subsidence and nudging tendencies
+! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
 ! 
 ! Former revisions:
@@ -95 +96 @@
 
     USE arrays_3d,                                                             &
-        ONLY :  ddzu, ddzw, e, hyp, km, kh,nr,  p, prho, pt, q, qc, ql, qr,    &
-                qs, qsws, qswst, rho, sa, saswsb, saswst, shf, ts, tswst, u,   &
-                ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, pt_lsa, pt_subs, q,&
+               qc, ql, qr, qs, qsws, qswst, q_lsa, q_subs, rho, sa, saswsb,    &
+               saswst, shf, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,&
+               v, vg, vpt, vswst, w, w_subs, zw
         
     USE cloud_parameters,                                                      &
@@ -104 +106 @@
     USE control_parameters,                                                    &
         ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
-                dt_3d, g, humidity, icloud_scheme, kappa, max_pr_user,         &
-                message_string, ocean, passive_scalar, precipitation,          &
-                use_surface_fluxes, use_top_fluxes, ws_scheme_mom, ws_scheme_sca
+                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
+                large_scale_subsidence, max_pr_user, message_string, ocean,    &
+                passive_scalar, precipitation, simulated_time,                 &
+                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
+                ws_scheme_mom, ws_scheme_sca
         
     USE cpulog,                                                                &
@@ -115 +119 @@
         
     USE indices,                                                               &
-        ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, nxl, &
-                nxr, nyn, nys, nzb, nzb_diff_s_inner, nzb_s_inner, nzt, nzt_diff
+        ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
+                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
+                nzb_s_inner, nzt, nzt_diff
         
     USE kinds
@@ -129 +134 @@
     INTEGER(iwp) ::  j                   !:
     INTEGER(iwp) ::  k                   !:
+    INTEGER(iwp) ::  nt                  !:
     INTEGER(iwp) ::  omp_get_thread_num  !:
     INTEGER(iwp) ::  sr                  !:
@@ -136 +142 @@
     
     REAL(wp) ::  dptdz_threshold  !: 
+    REAL(wp) ::  fac              !:
     REAL(wp) ::  height           !:
     REAL(wp) ::  pts              !:
@@ -982 +989 @@
 
 !
+!--    Collect current large scale advection and subsidence tendencies for
+!--    data output
+       IF ( large_scale_forcing )  THEN
+!
+!--       Interpolation in time of LSF_DATA 
+          nt = 1
+          DO WHILE ( simulated_time > time_vert(nt) )
+             nt = nt + 1
+          ENDDO
+          IF ( simulated_time /= time_vert(nt) )  THEN
+            nt = nt - 1
+          ENDIF
+
+          fac = ( simulated_time-time_vert(nt) )                               &
+                / ( time_vert(nt+1)-time_vert(nt) )
+
+
+          DO  k = nzb, nzt
+             sums_ls_l(k,0) = pt_lsa(k,nt)                                     &
+                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
+             sums_ls_l(k,1) = q_lsa(k,nt)                                      &
+                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
+          ENDDO
+
+          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
+
+             DO  k = nzb, nzt
+                sums_ls_l(k,2) = pt_subs(k,nt)                                 &
+                                 + fac * ( pt_subs(k,nt+1) - pt_subs(k,nt) )
+                sums_ls_l(k,3) = q_subs(k,nt)                                  &
+                                 + fac * ( q_subs(k,nt+1) - q_subs(k,nt) )
+             ENDDO
+
+          ENDIF
+
+       ENDIF
+
+!
 !--    Calculate the user-defined profiles
        CALL user_statistics( 'profiles', sr, tn )
@@ -1007 +1052 @@
 !--    Compute total sum from local sums
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
+       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
                            MPI_SUM, comm2d, ierr )
+       IF ( large_scale_forcing )  THEN
+          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
 #else
        sums = sums_l(:,:,0)
+       IF ( large_scale_forcing )  THEN
+          sums(:,81:88) = sums_ls_l
+       ENDIF
 #endif
 
@@ -1031 +1083 @@
           sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
           sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
-          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
+          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
+          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -1130 +1184 @@
        hom(:,1,80,sr) = w_subs         ! w_subs
 
+       IF ( large_scale_forcing )  THEN
+          hom(:,1,81,sr) = sums_ls_l(:,0)          ! pt_lsa
+          hom(:,1,82,sr) = sums_ls_l(:,1)          ! q_lsa
+          IF ( use_subsidence_tendencies )  THEN
+             hom(:,1,83,sr) = sums_ls_l(:,2)       ! pt_subs
+             hom(:,1,84,sr) = sums_ls_l(:,3)       ! q_subs
+          ELSE
+             hom(:,1,83,sr) = sums(:,83)           ! pt_subs
+             hom(:,1,84,sr) = sums(:,84)           ! q_subs
+          ENDIF
+          hom(:,1,85,sr) = sums(:,85)              ! pt_nudge
+          hom(:,1,86,sr) = sums(:,86)              ! q_nudge
+          hom(:,1,87,sr) = sums(:,87)              ! u_nudge
+          hom(:,1,88,sr) = sums(:,88)              ! v_nudge
+       ENDIF
+
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
                                        ! upstream-parts u_x, u_y, u_z, v_x,
@@ -1309 +1379 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh,nr,  p, prho, pt, q, qc, ql, qr,     &
-               qs, qsws, qswst, rho, sa, saswsb, saswst, shf, ts, tswst, u,    &
-               ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, pt_lsa, pt_subs, q,&
+               qc, ql, qr, qs, qsws, qswst, q_lsa, q_subs, rho, sa, saswsb,    &
+               saswst, shf, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,&
+               v, vg, vpt, vswst, w, w_subs, zw
+                 
         
     USE cloud_parameters,                                                      &
@@ -1317 +1389 @@
         
     USE control_parameters,                                                    &
-        ONLY:  average_count_pr, cloud_droplets, cloud_physics, do_sum,        &
-               dt_3d, g, humidity, icloud_scheme, kappa, max_pr_user,          &
-               message_string, ocean, passive_scalar, precipitation,           &
-               use_surface_fluxes, use_top_fluxes, ws_scheme_mom, ws_scheme_sca
+        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
+                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
+                large_scale_subsidence, max_pr_user, message_string, ocean,    &
+                passive_scalar, precipitation, simulated_time,                 &
+                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
+                ws_scheme_mom, ws_scheme_sca
         
     USE cpulog,                                                                &
@@ -1343 +1417 @@
     INTEGER(iwp) ::  j                   !:
     INTEGER(iwp) ::  k                   !:
+    INTEGER(iwp) ::  nt                  !:
     INTEGER(iwp) ::  omp_get_thread_num  !:
     INTEGER(iwp) ::  sr                  !:
@@ -1350 +1425 @@
     
     REAL(wp) ::  dptdz_threshold  !:
+    REAL(wp) ::  fac              !:
     REAL(wp) ::  height           !:
     REAL(wp) ::  pts              !:
@@ -2661 +2737 @@
 
 !
+!--    Collect current large scale advection and subsidence tendencies for
+!--    data output
+       IF ( large_scale_forcing )  THEN
+!
+!--       Interpolation in time of LSF_DATA 
+          nt = 1
+          DO WHILE ( simulated_time > time_vert(nt) )
+             nt = nt + 1
+          ENDDO
+          IF ( simulated_time /= time_vert(nt) )  THEN
+            nt = nt - 1
+          ENDIF
+
+          fac = ( simulated_time-time_vert(nt) )                               &
+                / ( time_vert(nt+1)-time_vert(nt) )
+
+
+          DO  k = nzb, nzt
+             sums_ls_l(k,0) = pt_lsa(k,nt)                                     &
+                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
+             sums_ls_l(k,1) = q_lsa(k,nt)                                      &
+                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
+          ENDDO
+
+          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
+
+             DO  k = nzb, nzt
+                sums_ls_l(k,2) = pt_subs(k,nt)                                 &
+                                 + fac * ( pt_subs(k,nt+1) - pt_subs(k,nt) )
+                sums_ls_l(k,3) = q_subs(k,nt)                                  &
+                                 + fac * ( q_subs(k,nt+1) - q_subs(k,nt) )
+             ENDDO
+
+          ENDIF
+
+       ENDIF
+
+!
 !--    Calculate the user-defined profiles
        CALL user_statistics( 'profiles', sr, tn )
@@ -2691 +2805 @@
        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
                            MPI_SUM, comm2d, ierr )
+       IF ( large_scale_forcing )  THEN
+          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
 #else
        sums = sums_l(:,:,0)
+       IF ( large_scale_forcing )  THEN
+          sums(:,81:88) = sums_ls_l
+       ENDIF
 #endif
 
@@ -2712 +2833 @@
           sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
           sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
-          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
-          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
+          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
+          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -2812 +2934 @@
        hom(:,1,80,sr) = w_subs         ! w_subs
 
+       IF ( large_scale_forcing )  THEN
+          hom(:,1,81,sr) = sums_ls_l(:,0)          ! pt_lsa
+          hom(:,1,82,sr) = sums_ls_l(:,1)          ! q_lsa
+          IF ( use_subsidence_tendencies )  THEN
+             hom(:,1,83,sr) = sums_ls_l(:,2)       ! pt_subs
+             hom(:,1,84,sr) = sums_ls_l(:,3)       ! q_subs
+          ELSE
+             hom(:,1,83,sr) = sums(:,83)           ! pt_subs
+             hom(:,1,84,sr) = sums(:,84)           ! q_subs
+          ENDIF
+          hom(:,1,85,sr) = sums(:,85)              ! pt_nudge
+          hom(:,1,86,sr) = sums(:,86)              ! q_nudge
+          hom(:,1,87,sr) = sums(:,87)              ! u_nudge
+          hom(:,1,88,sr) = sums(:,88)              ! v_nudge
+       END IF
+
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
                                        ! upstream-parts u_x, u_y, u_z, v_x,

palm/trunk/SOURCE/header.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! New section 'Large scale forcing and nudging':
+! output of large scale forcing and nudging information,
+! new section for initial profiles created
 ! 
 ! Former revisions:
@@ -98 +100 @@
 !
 ! 1015 2012-09-27 09:23:24Z raasch
-! output of Aajustment of mixing length to the Prandtl mixing length at first
+! output of Adjustment of mixing length to the Prandtl mixing length at first
 ! grid point above ground removed
 !
@@ -452 +454 @@
        ENDIF
     ENDIF
-    IF ( large_scale_subsidence )  THEN
-        WRITE ( io, 153 )
-        WRITE ( io, 154 )
-    ENDIF
-    IF ( nudging )  THEN
-        WRITE ( io, 155 )
-        WRITE ( io, 156 )
-    ENDIF
-    IF ( large_scale_forcing )  THEN
-        WRITE ( io, 157 )
-        WRITE ( io, 158 )
-    ENDIF
     WRITE ( io, 99 )
 
@@ -535 +525 @@
 
 !
+!-- Large scale forcing and nudging
+    WRITE ( io, 160 )
+    IF ( large_scale_forcing )  THEN
+       WRITE ( io, 162 )
+       WRITE ( io, 163 )
+
+       IF ( large_scale_subsidence )  THEN
+          IF ( .NOT. use_subsidence_tendencies )  THEN
+             WRITE ( io, 164 )
+          ELSE
+             WRITE ( io, 165 )
+          ENDIF 
+       ENDIF
+
+       IF ( bc_pt_b == 'dirichlet' )  THEN
+          WRITE ( io, 180 )
+       ELSEIF ( bc_pt_b == 'neumann' )  THEN
+          WRITE ( io, 181 )
+       ENDIF
+
+       IF ( bc_q_b == 'dirichlet' )  THEN
+          WRITE ( io, 182 )
+       ELSEIF ( bc_q_b == 'neumann' )  THEN
+          WRITE ( io, 183 )
+       ENDIF
+
+       WRITE ( io, 167 )
+       IF ( nudging )  THEN
+          WRITE ( io, 170 )
+       ENDIF
+    ELSE
+       WRITE ( io, 161 )
+       WRITE ( io, 171 )
+    ENDIF
+    IF ( large_scale_subsidence )  THEN
+       WRITE ( io, 168 )
+       WRITE ( io, 169 )
+    ENDIF
+
+!
+!-- Profile for the large scale vertial velocity
+!-- Building output strings, starting with surface value
+    IF ( large_scale_subsidence )  THEN
+       temperatures = '   0.0'
+       gradients = '------'
+       slices = '     0'
+       coordinates = '   0.0'
+       i = 1
+       DO  WHILE ( subs_vertical_gradient_level_i(i) /= -9999 )
+
+          WRITE (coor_chr,'(E10.2,7X)')  &
+                                w_subs(subs_vertical_gradient_level_i(i))
+          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(E10.2,7X)')  subs_vertical_gradient(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I10,7X)')  subs_vertical_gradient_level_i(i)
+          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  subs_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          IF ( i == 10 )  THEN
+             EXIT
+          ELSE
+             i = i + 1
+          ENDIF
+
+       ENDDO
+
+ 
+       IF ( .NOT. large_scale_forcing )  THEN
+          WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
+                             TRIM( gradients ), TRIM( slices )
+       ENDIF
+
+
+    ENDIF
+
+!-- Profile of the geostrophic wind (component ug)
+!-- Building output strings
+    WRITE ( ugcomponent, '(F6.2)' )  ug_surface
+    gradients = '------'
+    slices = '     0'
+    coordinates = '   0.0'
+    i = 1
+    DO  WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 )
+      
+       WRITE (coor_chr,'(F6.2,1X)')  ug(ug_vertical_gradient_level_ind(i))
+       ugcomponent = TRIM( ugcomponent ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F6.2,1X)')  ug_vertical_gradient(i)
+       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(I6,1X)')  ug_vertical_gradient_level_ind(i)
+       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F6.1,1X)')  ug_vertical_gradient_level(i)
+       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
+
+       IF ( i == 10 )  THEN
+          EXIT
+       ELSE
+          i = i + 1
+       ENDIF
+
+    ENDDO
+
+    IF ( .NOT. large_scale_forcing )  THEN
+       WRITE ( io, 423 )  TRIM( coordinates ), TRIM( ugcomponent ), &
+                          TRIM( gradients ), TRIM( slices )
+    ENDIF
+
+!-- Profile of the geostrophic wind (component vg)
+!-- Building output strings
+    WRITE ( vgcomponent, '(F6.2)' )  vg_surface
+    gradients = '------'
+    slices = '     0'
+    coordinates = '   0.0'
+    i = 1
+    DO  WHILE ( vg_vertical_gradient_level_ind(i) /= -9999 )
+
+       WRITE (coor_chr,'(F6.2,1X)')  vg(vg_vertical_gradient_level_ind(i))
+       vgcomponent = TRIM( vgcomponent ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F6.2,1X)')  vg_vertical_gradient(i)
+       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(I6,1X)')  vg_vertical_gradient_level_ind(i)
+       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F6.1,1X)')  vg_vertical_gradient_level(i)
+       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
+
+       IF ( i == 10 )  THEN
+          EXIT
+       ELSE
+          i = i + 1
+       ENDIF
+ 
+    ENDDO
+
+    IF ( .NOT. large_scale_forcing )  THEN
+       WRITE ( io, 424 )  TRIM( coordinates ), TRIM( vgcomponent ), &
+                          TRIM( gradients ), TRIM( slices )
+    ENDIF
+
+!
 !-- Topography
     WRITE ( io, 270 )  topography
@@ -805 +944 @@
        ENDIF
     ENDIF
+
+!
+!-- Initial Profiles
+    WRITE ( io, 321 )
+!
+!-- Initial wind profiles
+    IF ( u_profile(1) /= 9999999.9_wp )  WRITE ( io, 427 )
+
+!
+!-- Initial temperature profile
+!-- Building output strings, starting with surface temperature
+    WRITE ( temperatures, '(F6.2)' )  pt_surface
+    gradients = '------'
+    slices = '     0'
+    coordinates = '   0.0'
+    i = 1
+    DO  WHILE ( pt_vertical_gradient_level_ind(i) /= -9999 )
+
+       WRITE (coor_chr,'(F7.2)')  pt_init(pt_vertical_gradient_level_ind(i))
+       temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F7.2)')  pt_vertical_gradient(i)
+       gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(I7)')  pt_vertical_gradient_level_ind(i)
+       slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+       WRITE (coor_chr,'(F7.1)')  pt_vertical_gradient_level(i)
+       coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+       IF ( i == 10 )  THEN
+          EXIT
+       ELSE
+          i = i + 1
+       ENDIF
+
+    ENDDO
+
+    IF ( .NOT. nudging )  THEN
+       WRITE ( io, 420 )  TRIM( coordinates ), TRIM( temperatures ), &
+                          TRIM( gradients ), TRIM( slices )
+    ELSE
+       WRITE ( io, 428 ) 
+    ENDIF
+
+!
+!-- Initial humidity profile
+!-- Building output strings, starting with surface humidity
+    IF ( humidity  .OR.  passive_scalar )  THEN
+       WRITE ( temperatures, '(E8.1)' )  q_surface
+       gradients = '--------'
+       slices = '       0'
+       coordinates = '     0.0'
+       i = 1
+       DO  WHILE ( q_vertical_gradient_level_ind(i) /= -9999 )
+          
+          WRITE (coor_chr,'(E8.1,4X)')  q_init(q_vertical_gradient_level_ind(i))
+          temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(E8.1,4X)')  q_vertical_gradient(i)
+          gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
+          
+          WRITE (coor_chr,'(I8,4X)')  q_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
+          
+          WRITE (coor_chr,'(F8.1,4X)')  q_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
+
+          IF ( i == 10 )  THEN
+             EXIT
+          ELSE
+             i = i + 1
+          ENDIF
+
+       ENDDO
+
+       IF ( humidity )  THEN
+          IF ( .NOT. nudging )  THEN
+             WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
+                                TRIM( gradients ), TRIM( slices )
+          ENDIF
+       ELSE
+          WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
+                             TRIM( gradients ), TRIM( slices )
+       ENDIF
+    ENDIF
+
+!
+!-- Initial salinity profile
+!-- Building output strings, starting with surface salinity
+    IF ( ocean )  THEN
+       WRITE ( temperatures, '(F6.2)' )  sa_surface
+       gradients = '------'
+       slices = '     0'
+       coordinates = '   0.0'
+       i = 1
+       DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
+
+          WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
+          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          IF ( i == 10 )  THEN
+             EXIT
+          ELSE
+             i = i + 1
+          ENDIF
+
+       ENDDO
+
+       WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
+                          TRIM( gradients ), TRIM( slices )
+    ENDIF
+
 
 !
@@ -1268 +1529 @@
           IF ( precipitation )  WRITE ( io, 511 ) c_sedimentation
        ENDIF
-    ENDIF
-
-!-- Profile of the geostrophic wind (component ug)
-!-- Building output strings
-    WRITE ( ugcomponent, '(F6.2)' )  ug_surface
-    gradients = '------'
-    slices = '     0'
-    coordinates = '   0.0'
-    i = 1
-    DO  WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 )
-      
-       WRITE (coor_chr,'(F6.2,1X)')  ug(ug_vertical_gradient_level_ind(i))
-       ugcomponent = TRIM( ugcomponent ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.2,1X)')  ug_vertical_gradient(i)
-       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(I6,1X)')  ug_vertical_gradient_level_ind(i)
-       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.1,1X)')  ug_vertical_gradient_level(i)
-       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
-
-       IF ( i == 10 )  THEN
-          EXIT
-       ELSE
-          i = i + 1
-       ENDIF
-
-    ENDDO
-
-    IF ( .NOT. large_scale_forcing )  THEN
-       WRITE ( io, 423 )  TRIM( coordinates ), TRIM( ugcomponent ), &
-                          TRIM( gradients ), TRIM( slices )
-    ELSE
-       WRITE ( io, 429 ) 
-    ENDIF
-
-!-- Profile of the geostrophic wind (component vg)
-!-- Building output strings
-    WRITE ( vgcomponent, '(F6.2)' )  vg_surface
-    gradients = '------'
-    slices = '     0'
-    coordinates = '   0.0'
-    i = 1
-    DO  WHILE ( vg_vertical_gradient_level_ind(i) /= -9999 )
-
-       WRITE (coor_chr,'(F6.2,1X)')  vg(vg_vertical_gradient_level_ind(i))
-       vgcomponent = TRIM( vgcomponent ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.2,1X)')  vg_vertical_gradient(i)
-       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(I6,1X)')  vg_vertical_gradient_level_ind(i)
-       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F6.1,1X)')  vg_vertical_gradient_level(i)
-       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
-
-       IF ( i == 10 )  THEN
-          EXIT
-       ELSE
-          i = i + 1
-       ENDIF
- 
-    ENDDO
-
-    IF ( .NOT. large_scale_forcing )  THEN
-       WRITE ( io, 424 )  TRIM( coordinates ), TRIM( vgcomponent ), &
-                          TRIM( gradients ), TRIM( slices )
-    ENDIF
-
-!
-!-- Initial wind profiles
-    IF ( u_profile(1) /= 9999999.9_wp )  WRITE ( io, 427 )
-
-!
-!-- Initial temperature profile
-!-- Building output strings, starting with surface temperature
-    WRITE ( temperatures, '(F6.2)' )  pt_surface
-    gradients = '------'
-    slices = '     0'
-    coordinates = '   0.0'
-    i = 1
-    DO  WHILE ( pt_vertical_gradient_level_ind(i) /= -9999 )
-
-       WRITE (coor_chr,'(F7.2)')  pt_init(pt_vertical_gradient_level_ind(i))
-       temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F7.2)')  pt_vertical_gradient(i)
-       gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(I7)')  pt_vertical_gradient_level_ind(i)
-       slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
-
-       WRITE (coor_chr,'(F7.1)')  pt_vertical_gradient_level(i)
-       coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
-
-       IF ( i == 10 )  THEN
-          EXIT
-       ELSE
-          i = i + 1
-       ENDIF
-
-    ENDDO
-
-    IF ( .NOT. nudging )  THEN
-       WRITE ( io, 420 )  TRIM( coordinates ), TRIM( temperatures ), &
-                          TRIM( gradients ), TRIM( slices )
-    ELSE
-       WRITE ( io, 428 ) 
-    ENDIF
-
-!
-!-- Initial humidity profile
-!-- Building output strings, starting with surface humidity
-    IF ( humidity  .OR.  passive_scalar )  THEN
-       WRITE ( temperatures, '(E8.1)' )  q_surface
-       gradients = '--------'
-       slices = '       0'
-       coordinates = '     0.0'
-       i = 1
-       DO  WHILE ( q_vertical_gradient_level_ind(i) /= -9999 )
-          
-          WRITE (coor_chr,'(E8.1,4X)')  q_init(q_vertical_gradient_level_ind(i))
-          temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(E8.1,4X)')  q_vertical_gradient(i)
-          gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-          
-          WRITE (coor_chr,'(I8,4X)')  q_vertical_gradient_level_ind(i)
-          slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-          
-          WRITE (coor_chr,'(F8.1,4X)')  q_vertical_gradient_level(i)
-          coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
-
-          IF ( i == 10 )  THEN
-             EXIT
-          ELSE
-             i = i + 1
-          ENDIF
-
-       ENDDO
-
-       IF ( humidity )  THEN
-          IF ( .NOT. nudging )  THEN
-             WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
-                                TRIM( gradients ), TRIM( slices )
-          ENDIF
-       ELSE
-          WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
-                             TRIM( gradients ), TRIM( slices )
-       ENDIF
-    ENDIF
-
-!
-!-- Initial salinity profile
-!-- Building output strings, starting with surface salinity
-    IF ( ocean )  THEN
-       WRITE ( temperatures, '(F6.2)' )  sa_surface
-       gradients = '------'
-       slices = '     0'
-       coordinates = '   0.0'
-       i = 1
-       DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
-
-          WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
-          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
-          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
-          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
-          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
-
-          IF ( i == 10 )  THEN
-             EXIT
-          ELSE
-             i = i + 1
-          ENDIF
-
-       ENDDO
-
-       WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
-                          TRIM( gradients ), TRIM( slices )
-    ENDIF
-
-!
-!-- Profile for the large scale vertial velocity
-!-- Building output strings, starting with surface value
-    IF ( large_scale_subsidence )  THEN
-       temperatures = '   0.0'
-       gradients = '------'
-       slices = '     0'
-       coordinates = '   0.0'
-       i = 1
-       DO  WHILE ( subs_vertical_gradient_level_i(i) /= -9999 )
-
-          WRITE (coor_chr,'(E10.2,7X)')  &
-                                w_subs(subs_vertical_gradient_level_i(i))
-          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(E10.2,7X)')  subs_vertical_gradient(i)
-          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(I10,7X)')  subs_vertical_gradient_level_i(i)
-          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
-
-          WRITE (coor_chr,'(F10.2,7X)')  subs_vertical_gradient_level(i)
-          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
-
-          IF ( i == 10 )  THEN
-             EXIT
-          ELSE
-             i = i + 1
-          ENDIF
-
-       ENDDO
-
- 
-       IF ( .NOT. large_scale_forcing )  THEN
-          WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
-                             TRIM( gradients ), TRIM( slices )
-        ELSE
-           WRITE ( io, 460 ) 
-       ENDIF
-
-
     ENDIF
 
@@ -1727 +1757 @@
            /'     ',2(1X,E12.5),'Pa/m in x/y direction', &
            /'     starting from dp_level_b =', F8.3, 'm', A /)
-153 FORMAT (' --> Large-scale vertical motion is used in the ', &
+160 FORMAT (//' Large scale forcing and nudging:'/ &
+              ' -------------------------------'/)
+161 FORMAT (' --> No large scale forcing from external is used (default) ')
+162 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
+163 FORMAT ('     - large scale advection tendencies ')
+164 FORMAT ('     - large scale subsidence velocity w_subs ')
+165 FORMAT ('     - large scale subsidence tendencies ')
+167 FORMAT ('     - and geostrophic wind components ug and vg')
+168 FORMAT (' --> Large-scale vertical motion is used in the ', &
                   'prognostic equation(s) for')
-154 FORMAT ('     the scalar(s) only')
-155 FORMAT (' --> Nudging is used - initial profiles for u, v, pt and q ',& 
-                  'correspond to the')
-156 FORMAT ('     first profiles in NUDGING_DATA')
-157 FORMAT (' --> Large scale forcing from external file (LSF_DATA) is used: ')
-158 FORMAT ('     prescribed surfaces fluxes and geostrophic wind components')
+169 FORMAT ('     the scalar(s) only')
+170 FORMAT (' --> Nudging is used')
+171 FORMAT (' --> No nudging is used (default) ')
+180 FORMAT ('     - prescribed surface values for temperature')
+181 FORMAT ('     - prescribed surface fluxes for humidity')
+182 FORMAT ('     - prescribed surface values for temperature')
+183 FORMAT ('     - prescribed surface fluxes for humidity')
 200 FORMAT (//' Run time and time step information:'/ &
              ' ----------------------------------'/)
@@ -1832 +1871 @@
 320 FORMAT ('       Predefined constant momentumflux:  u: ',F9.6,' m**2/s**2'/ &
             '                                          v: ',F9.6,' m**2/s**2')
+321 FORMAT (//' Initial profiles:'/ &
+              ' ----------------')
 325 FORMAT (//' List output:'/ &
              ' -----------'//  &
@@ -1989 +2030 @@
 428 FORMAT (/'    Initial profiles (u, v, pt, q) are taken from file '/ &
              '    NUDGING_DATA')
-429 FORMAT (/'    Geostrophic wind profiles (ug, vg) are are taken from file '/ &
-             '    LSF_DATA')
 430 FORMAT (//' Cloud physics quantities / methods:'/ &
               ' ----------------------------------'/)
@@ -2017 +2056 @@
 454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
 455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
-460 FORMAT (/'    Profiles for large scale vertical velocity are '/ &
-             '    taken from file LSF_DATA')
 470 FORMAT (//' Actions during the simulation:'/ &
               ' -----------------------------'/)

palm/trunk/SOURCE/ls_forcing.f90

@@ -20 +20 @@
 ! 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 +65 @@
 
     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 +78 @@
 
        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 +86 @@
 
        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 +166 @@
 !
 !--    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 +186 @@
        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 +201 @@
 !--    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 +210 @@
 !--       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 +218 @@
           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 +228 @@
           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 +245 @@
                 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 +270 @@
              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 +288 @@
  
        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 +313 @@
        USE kinds
 
-
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  t                     !:
+       INTEGER(iwp) ::  nt                     !:
 
        REAL(wp)             :: fac            !:
@@ -276 +322 @@
 !
 !--    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 +336 @@
 !--       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 +342 @@
 !--       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 +351 @@
 !--       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 +374 @@
        USE kinds
 
-
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  t                     !:
+       INTEGER(iwp) ::  nt                     !:
 
        REAL(wp)             ::  fac           !:
@@ -338 +383 @@
 !
 !--    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 +404 @@
 
 
+!------------------------------------------------------------------------------!
+! 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

palm/trunk/SOURCE/modules.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Usage of large scale forcing enabled:
+! increase pr_palm from 90 to 100 to allow for more standard profiles
+! + ngp_sums_ls, pt_lsa, pt_subs, q_lsa, q_subs, tmp_tnudge, sums_ls_l, 
+! use_subsidence_tendencies
 ! 
 ! Former revisions:
@@ -273 +276 @@
           nc_1d, nr_1d, ptdf_x, ptdf_y, p_surf, pt_surf, pt_1d, pt_init,       &
           qsws_surf, q_1d, q_init, q_surf, qc_1d, qr_1d, rdf, rdf_sc,          & 
-          ref_state, sa_init, shf_surf, timenudge, time_surf, time_vert, ug,   &
-          u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
+          ref_state, sa_init, shf_surf, timenudge, time_surf, time_vert,       &
+          tmp_tnudge, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,                &
+          v_nzb_p1_for_vfc, w_subs, zu, zw
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
@@ -281 +285 @@
           flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
           flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
-          mean_inflow_profiles, nrs, nrsws, nrswst, ptnudge, pt_slope_ref,     &
-          qnudge, qs, qsws, qswst, qswst_remote, qrs, qrsws, qrswst, rif,      &
-          saswsb, saswst, shf, tnudge, total_2d_a, total_2d_o, ts, tswst,      &
-          ug_vert, unudge, us, usws, uswst, vnudge, vg_vert, vsws, vswst,      &
-          wnudge, wsubs_vert, z0, z0h
+          mean_inflow_profiles, nrs, nrsws, nrswst, ptnudge, pt_lsa,           &
+          pt_slope_ref, pt_subs, qnudge, qs, qsws, qswst, qswst_remote, qrs,   &
+          qrsws, qrswst, q_lsa, q_subs, rif, saswsb, saswst, shf, tnudge,      &
+          total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us, usws, uswst, &
+          vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, z0, z0h
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
@@ -638 +642 @@
                 use_prescribed_profile_data = .FALSE., &
                 use_single_reference_value = .FALSE., &
+                use_subsidence_tendencies = .FALSE., &
                 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
                 use_ug_for_galilei_tr = .TRUE., use_upstream_for_tke = .FALSE.,&
@@ -921 +926 @@
     USE kinds
 
-    INTEGER(iwp) ::  i_left, i_right, j_north, j_south, nbgp = 3, ngp_sums, nnx,    &
-                     nx = 0, nx_a, nx_o, nxl, nxlg, nxlu, nxr, nxrg, nx_on_file,    &
-                     nny, ny = 0, ny_a, ny_o, nyn, nyng, nys, nysg, nysv,           &
-                     ny_on_file, nnz, nz = 0, nzb, nzb_diff, nzb_max, nzt, nzt_diff
+    INTEGER(iwp) ::  i_left, i_right, j_north, j_south, nbgp = 3, ngp_sums,         &
+                     ngp_sums_ls, nnx, nx = 0, nx_a, nx_o, nxl, nxlg, nxlu, nxr,    &
+                     nxrg, nx_on_file, nny, ny = 0, ny_a, ny_o, nyn, nyng, nys,     &
+                     nysg, nysv, ny_on_file, nnz, nz = 0, nzb, nzb_diff, nzb_max,   &
+                     nzt, nzt_diff
 
 
@@ -1389 +1395 @@
 
     CHARACTER (LEN=40) ::  region(0:9)
-    INTEGER(iwp) ::  pr_palm = 90, statistic_regions = 0
+    INTEGER(iwp) ::  pr_palm = 100, statistic_regions = 0
     INTEGER(iwp) ::  u_max_ijk(3) = -1, v_max_ijk(3) = -1, w_max_ijk(3) = -1
     LOGICAL ::  flow_statistics_called = .FALSE.
@@ -1403 +1409 @@
                                                   sums_wssas_ws_l,                &
                                                   sums_wsqs_ws_l,                 &
-                                                  sums_wsqrs_ws_l
+                                                  sums_wsqrs_ws_l,                &
+                                                  sums_ls_l
                                               
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  hom_sum, rmask, spectrum_x, &

palm/trunk/SOURCE/nudging.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Variable t renamed nt, variable currtnudge renamed tmp_tnudge,
+! summation of nudging tendencies for data output added
+! +sums_ls_l, tmp_tend
+! Added new subroutine calc_tnudge, which calculates the current nudging time 
+! scale at each time step
 ! 
 ! Former revisions:
@@ -65 +69 @@
 
     PRIVATE
-    PUBLIC init_nudge, nudge 
+    PUBLIC init_nudge, calc_tnudge, nudge 
     SAVE
 
@@ -78 +82 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ptnudge, qnudge, timenudge, tnudge, unudge, vnudge, wnudge,  &
-                  zu
+           ONLY:  ptnudge, qnudge, timenudge, tmp_tnudge, tnudge, unudge,      &
+                  vnudge, wnudge, zu
 
        USE control_parameters,                                                 &
@@ -96 +100 @@
        INTEGER(iwp) ::  ierrn        !:
        INTEGER(iwp) ::  k            !:
-       INTEGER(iwp) ::  t            !:
+       INTEGER(iwp) ::  nt            !:
 
        CHARACTER(1) ::  hash     !:
@@ -122 +126 @@
                  vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
 
+       ALLOCATE( tmp_tnudge(nzb:nzt) )
+
        ALLOCATE( timenudge(0:ntnudge) )
-
 
        ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
        vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
-
-       t = 0
+!
+!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
+       tmp_tnudge = 21600.0_wp
+
+       nt = 0
        OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
               FORM='FORMATTED', IOSTAT=ierrn )
@@ -140 +148 @@
 
  rloop:DO
-          t = t + 1
+          nt = nt + 1
           hash = "#"
           ierrn = 1 ! not zero
@@ -148 +156 @@
           DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
           
-            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(t)
+            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
             IF ( ierrn < 0 )  EXIT rloop
 
@@ -202 +210 @@
              fac = ( highheight - zu(k) ) / ( highheight - lowheight )
 
-             tnudge(k,t)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
-             unudge(k,t)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
-             vnudge(k,t)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
-             wnudge(k,t)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
-             ptnudge(k,t) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
-             qnudge(k,t)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
+             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
+             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
+             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
+             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
+             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
+             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
           ENDDO
 
@@ -225 +233 @@
     END SUBROUTINE init_nudge
 
+
+    SUBROUTINE calc_tnudge ( time )
+
+       USE arrays_3d,                                                          &
+           ONLY:  timenudge, tmp_tnudge, tnudge
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_3d 
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+
+       REAL(wp) ::  dtm         !:
+       REAL(wp) ::  dtp         !:
+       REAL(wp) ::  time        !:
+
+       INTEGER(iwp) ::  k   !:
+       INTEGER(iwp) ::  nt  !:
+
+       nt = 1
+       DO WHILE ( time > timenudge(nt) )
+         nt = nt+1
+       ENDDO
+       IF ( time /= timenudge(1) ) THEN
+         nt = nt-1
+       ENDIF
+
+       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
+       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
+
+       DO  k = nzb, nzt
+          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
+       ENDDO
+
+    END SUBROUTINE calc_tnudge
+
 !------------------------------------------------------------------------------!
 ! Call for all grid points
@@ -231 +280 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tnudge, u, unudge,  &
-                  v, vnudge
-
-       USE buoyancy_mod,                                                       &
-           ONLY:  calc_mean_profile
+           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
+                  u, unudge, v, vnudge
 
        USE control_parameters,                                                 &
-           ONLY:  dt_3d, message_string
+           ONLY:  dt_3d, intermediate_timestep_count, message_string
 
        USE indices,                                                            &
            ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
 
-       USE kinds,                                                              &
-           ONLY:  iwp, wp
+       USE kinds
 
        USE statistics,                                                         &
-           ONLY:  hom
+           ONLY:  hom, sums_ls_l, weight_pres
 
        IMPLICIT NONE
@@ -253 +298 @@
        CHARACTER (LEN=*) ::  prog_var  !:
 
-       REAL(wp) ::  currtnudge  !:
+       REAL(wp) ::  tmp_tend    !:
        REAL(wp) ::  dtm         !:
        REAL(wp) ::  dtp         !:
@@ -261 +306 @@
        INTEGER(iwp) ::  j  !:
        INTEGER(iwp) ::  k  !:
-       INTEGER(iwp) ::  t  !:
-
-
-       t = 1
-       DO WHILE ( time > timenudge(t) )
-         t = t+1
+       INTEGER(iwp) ::  nt  !:
+
+
+       nt = 1
+       DO WHILE ( time > timenudge(nt) )
+         nt = nt+1
        ENDDO
        IF ( time /= timenudge(1) ) THEN
-         t = t-1
+         nt = nt-1
        ENDIF
 
-       dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
-       dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )
+       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
+       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
 
        SELECT CASE ( prog_var )
 
           CASE ( 'u' )
-
-             CALL calc_mean_profile( u, 1, 'nudging' )
 
              DO  i = nxl, nxr
@@ -285 +328 @@
                    DO  k = nzb_u_inner(j,i)+1, nzt
 
-                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
-                                               tnudge(k,t+1) * dtm )
-
-
-                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
-                                    ( unudge(k,t) * dtp + &
-                                       unudge(k,t+1) * dtm ) ) / currtnudge
+                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
+                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                      tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *             &
+                                       weight_pres(intermediate_timestep_count)
                    ENDDO
                 ENDDO
@@ -298 +341 @@
           CASE ( 'v' )
 
-             CALL calc_mean_profile( v, 2, 'nudging' )
-
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_u_inner(j,i)+1, nzt
 
-                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
-                                               tnudge(k,t+1) * dtm )
-
-
-                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
-                                    ( vnudge(k,t) * dtp + &
-                                       vnudge(k,t+1) * dtm ) ) / currtnudge
+                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
+                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                      tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *             &
+                                       weight_pres(intermediate_timestep_count)
                    ENDDO
                 ENDDO
@@ -317 +358 @@
           CASE ( 'pt' )
 
-             CALL calc_mean_profile( pt, 4, 'nudging' )
-
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_u_inner(j,i)+1, nzt
 
-                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
-                                               tnudge(k,t+1) * dtm )
-
-
-                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
-                                    ( ptnudge(k,t) * dtp + &
-                                       ptnudge(k,t+1) * dtm ) ) / currtnudge
+                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
+                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                      tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *             &
+                                       weight_pres(intermediate_timestep_count)
                    ENDDO
                 ENDDO
@@ -336 +375 @@
           CASE ( 'q' )
 
-             CALL calc_mean_profile( q, 41, 'nudging' )
-
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_u_inner(j,i)+1, nzt
 
-                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
-                                               tnudge(k,t+1) * dtm )
-
-
-                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
-                                    ( qnudge(k,t) * dtp + &
-                                       qnudge(k,t+1) * dtm ) ) / currtnudge
+                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
+                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                      tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *             &
+                                       weight_pres(intermediate_timestep_count)
                    ENDDO
                 ENDDO
@@ -369 +406 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tnudge, u, unudge,  &
-                  v, vnudge
-
-       USE buoyancy_mod,                                                       &
-           ONLY:  calc_mean_profile
+           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
+                  u, unudge, v, vnudge
 
        USE control_parameters,                                                 &
-           ONLY:  dt_3d, message_string
+           ONLY:  dt_3d, intermediate_timestep_count, message_string
 
        USE indices,                                                            &
            ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
 
-       USE kinds,                                                              &
-           ONLY:  iwp, wp
+       USE kinds
 
        USE statistics,                                                         &
-           ONLY:  hom
+           ONLY:  hom, sums_ls_l, weight_pres
 
        IMPLICIT NONE
@@ -392 +425 @@
        CHARACTER (LEN=*) ::  prog_var  !:
 
-       REAL(wp) ::  currtnudge  !:
+       REAL(wp) ::  tmp_tend    !:
        REAL(wp) ::  dtm         !:
        REAL(wp) ::  dtp         !:
@@ -400 +433 @@
        INTEGER(iwp) ::  j  !:
        INTEGER(iwp) ::  k  !:
-       INTEGER(iwp) ::  t  !:
-
-
-       t = 1
-       DO WHILE ( time > timenudge(t) )
-         t = t+1
+       INTEGER(iwp) ::  nt  !:
+
+
+       nt = 1
+       DO WHILE ( time > timenudge(nt) )
+         nt = nt+1
        ENDDO
        IF ( time /= timenudge(1) )  THEN
-         t = t-1
+         nt = nt-1
        ENDIF
 
-       dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
-       dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )
+       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
+       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
 
        SELECT CASE ( prog_var )
@@ -418 +451 @@
           CASE ( 'u' )
 
-             CALL calc_mean_profile( u, 1, 'nudging' )
-
              DO  k = nzb_u_inner(j,i)+1, nzt
 
-                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
-
-                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
-                               ( unudge(k,t) * dtp + &
-                                 unudge(k,t+1) * dtm ) ) / currtnudge
+                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
+                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
+                                 * weight_pres(intermediate_timestep_count)
              ENDDO
 
           CASE ( 'v' )
 
-             CALL calc_mean_profile( v, 2, 'nudging' )
-
              DO  k = nzb_u_inner(j,i)+1, nzt
 
-                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
-
-                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
-                               ( vnudge(k,t) * dtp + &
-                                 vnudge(k,t+1) * dtm ) ) / currtnudge
+                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
+                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
+                                 * weight_pres(intermediate_timestep_count)
              ENDDO
 
@@ -445 +478 @@
           CASE ( 'pt' )
 
-             CALL calc_mean_profile( pt, 4, 'nudging' )
-
              DO  k = nzb_u_inner(j,i)+1, nzt
 
-                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
-
-                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
-                               ( ptnudge(k,t) * dtp + &
-                                 ptnudge(k,t+1) * dtm ) ) / currtnudge
+                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
+                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
+                                 * weight_pres(intermediate_timestep_count)
              ENDDO
 
@@ -459 +492 @@
           CASE ( 'q' )
 
-             CALL calc_mean_profile( q, 41, 'nudging' )
-
              DO  k = nzb_u_inner(j,i)+1, nzt
 
-                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
-
-                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
-                               ( qnudge(k,t) * dtp + &
-                                 qnudge(k,t+1) * dtm ) ) / currtnudge
+                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
+                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
+
+                tend(k,j,i) = tend(k,j,i) + tmp_tend
+
+                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
+                                 * weight_pres(intermediate_timestep_count)
              ENDDO
 

palm/trunk/SOURCE/parin.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Usage of large scale forcing enabled:
+! +use_subsidence_tendencies
 ! 
 ! Former revisions:
@@ -211 +212 @@
                transpose_compute_overlap, turbulence, turbulent_inflow,        &
                ug_surface, ug_vertical_gradient, ug_vertical_gradient_level,   &
-               use_surface_fluxes, use_cmax, use_top_fluxes,                   &
-               use_ug_for_galilei_tr, use_upstream_for_tke, uv_heights,        &
-               u_bulk, u_profile, vg_surface, vg_vertical_gradient,            &
+               use_subsidence_tendencies, use_surface_fluxes, use_cmax,        &
+               use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke,    &
+               uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient,&
                vg_vertical_gradient_level, v_bulk, v_profile,                  &
                wall_adjustment, wall_heatflux, wall_humidityflux,              &
@@ -295 +296 @@
              top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
              transpose_compute_overlap, turbulence, turbulent_inflow, &
-             ug_surface, ug_vertical_gradient, &
+             use_subsidence_tendencies, ug_surface, ug_vertical_gradient, &
              ug_vertical_gradient_level, use_surface_fluxes, use_cmax, &
              use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &

palm/trunk/SOURCE/prognostic_equations.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Calls of ls_advec for large scale advection added, 
+! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
+! new argument ls_index added to the calls of subsidence
+! +ls_index
 ! 
 ! Former revisions:
@@ -143 +146 @@
                dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
                icloud_scheme, inflow_l, intermediate_timestep_count,           &
-               intermediate_timestep_count_max, large_scale_subsidence,        &
-               neutral, nudging, ocean, outflow_l, outflow_s, passive_scalar,  &
-               plant_canopy, precipitation, prho_reference, prho_reference,    &
-               prho_reference, pt_reference, pt_reference, pt_reference,       &
-               radiation, scalar_advec, scalar_advec, simulated_time,          &
-               sloping_surface, timestep_scheme, tsc, use_upstream_for_tke,    &
+               intermediate_timestep_count_max, large_scale_forcing,           &
+               large_scale_subsidence, neutral, nudging, ocean, outflow_l,     &
+               outflow_s, passive_scalar, plant_canopy, precipitation,         &
+               prho_reference, prho_reference, prho_reference, pt_reference,   &
+               pt_reference, pt_reference, radiation, scalar_advec,            &
+               scalar_advec, simulated_time, sloping_surface, timestep_scheme, &
+               tsc, use_subsidence_tendencies, use_upstream_for_tke,           &
                use_upstream_for_tke, use_upstream_for_tke, wall_heatflux,      &
                wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
@@ -222 +226 @@
 
     USE kinds
+
+    USE ls_forcing_mod,                                                        &
+        ONLY:  ls_advec
 
     USE microphysics_mod,                                                      &
@@ -307 +314 @@
           i_omp_start = i 
        ENDIF
- 
+
        DO  j = nys, nyn
 !
@@ -537 +544 @@
 
 !
+!--          Large scale advection
+             IF ( large_scale_forcing )  THEN
+                CALL ls_advec( i, j, simulated_time, 'pt' )
+             ENDIF     
+
+!
 !--          If required, compute effect of large-scale subsidence/ascent
-             IF ( large_scale_subsidence )  THEN
-                CALL subsidence( i, j, tend, pt, pt_init )
+             IF ( large_scale_subsidence  .AND.                                &
+                  .NOT. use_subsidence_tendencies )  THEN
+                CALL subsidence( i, j, tend, pt, pt_init, 2 )
              ENDIF
 
@@ -662 +676 @@
 
 !
+!--          Large scale advection
+             IF ( large_scale_forcing )  THEN
+                CALL ls_advec( i, j, simulated_time, 'q' )
+             ENDIF
+
+!
 !--          If required compute influence of large-scale subsidence/ascent
-             IF ( large_scale_subsidence )  THEN
-                CALL subsidence( i, j, tend, q, q_init )
+             IF ( large_scale_subsidence  .AND.                                &
+                  .NOT. use_subsidence_tendencies )  THEN
+                CALL subsidence( i, j, tend, q, q_init, 3 )
              ENDIF
 
@@ -1177 +1198 @@
 
 !
+!--    Large scale advection
+       IF ( large_scale_forcing )  THEN
+          CALL ls_advec( simulated_time, 'pt' )
+       ENDIF
+
+!
 !--    If required compute influence of large-scale subsidence/ascent
-       IF ( large_scale_subsidence )  THEN
-          CALL subsidence( tend, pt, pt_init )
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies )  THEN
+          CALL subsidence( tend, pt, pt_init, 2 )
        ENDIF
 
@@ -1365 +1393 @@
 !--    Sink or source of scalar concentration due to canopy elements
        IF ( plant_canopy ) CALL plant_canopy_model( 5 )
-       
+
+!
+!--    Large scale advection
+       IF ( large_scale_forcing )  THEN
+          CALL ls_advec( simulated_time, 'q' )
+       ENDIF
+
 !
 !--    If required compute influence of large-scale subsidence/ascent
-       IF ( large_scale_subsidence )  THEN
-         CALL subsidence( tend, q, q_init )
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies )  THEN
+         CALL subsidence( tend, q, q_init, 3 )
        ENDIF
 
@@ -1978 +2013 @@
 
 !
+!--    Large scale advection
+       IF ( large_scale_forcing )  THEN
+          CALL ls_advec( simulated_time, 'pt' )
+       ENDIF
+
+!
 !--    If required compute influence of large-scale subsidence/ascent
-       IF ( large_scale_subsidence )  THEN
-          CALL subsidence( tend, pt, pt_init )
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies )  THEN
+          CALL subsidence( tend, pt, pt_init, 2 )
        ENDIF
 
@@ -2142 +2184 @@
 
 !
+!--    Large scale advection
+       IF ( large_scale_forcing )  THEN
+          CALL ls_advec( simulated_time, 'q' )
+       ENDIF
+
+!
 !--    If required compute influence of large-scale subsidence/ascent
-       IF ( large_scale_subsidence )  THEN
-         CALL subsidence( tend, q, q_init )
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies )  THEN
+         CALL subsidence( tend, q, q_init, 3 )
        ENDIF
 

palm/trunk/SOURCE/subsidence.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Summation of subsidence tendencies for data output added 
+! +ls_index, sums_ls_l, tmp_tend
 ! 
 ! Former revisions:
@@ -90 +91 @@
        REAL(wp)     ::  ws_surface !: 
 
-       IF ( .NOT. ALLOCATED( w_subs )) THEN
+       IF ( .NOT. ALLOCATED( w_subs ))  THEN
           ALLOCATE( w_subs(nzb:nzt+1) )
           w_subs = 0.0_wp
        ENDIF 
 
-      IF ( ocean )  THEN
+       IF ( ocean )  THEN
           message_string = 'Applying large scale vertical motion is not ' // &
                            'allowed for ocean runs'
@@ -111 +112 @@
       subs_vertical_gradient_level_i(1) = 0
       DO  k = 1, nzt+1
-         IF ( i < 11 ) THEN
+         IF ( i < 11 )  THEN
             IF ( subs_vertical_gradient_level(i) < zu(k)  .AND. &
                  subs_vertical_gradient_level(i) >= 0.0_wp )  THEN
@@ -140 +141 @@
 
 
-    SUBROUTINE subsidence( tendency, var, var_init ) 
+    SUBROUTINE subsidence( tendency, var, var_init, ls_index ) 
 
        USE arrays_3d,                                                          &
@@ -146 +147 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dt_3d
+           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing
 
        USE indices,                                                            &
@@ -153 +154 @@
 
        USE kinds
+
+       USE statistics,                                                         &
+           ONLY:  sums_ls_l, weight_pres
 
        IMPLICIT NONE
@@ -159 +163 @@
        INTEGER(iwp) ::  j !: 
        INTEGER(iwp) ::  k !: 
-
+       INTEGER(iwp) ::  ls_index !:
+
+       REAL(wp)     ::  tmp_tend !:
        REAL(wp)     ::  tmp_grad !: 
     
@@ -174 +180 @@
           DO  j = nys, nyn
              DO  k = nzb_s_inner(j,i)+1, nzt 
-                IF ( w_subs(k) < 0.0_wp ) THEN    ! large-scale subsidence
-                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
-                                     ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
-                ELSE                           ! large-scale ascent
-                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
-                                     ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+                IF ( w_subs(k) < 0.0_wp )  THEN    ! large-scale subsidence
+                   tmp_tend = - w_subs(k) *                                    &
+                              ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
+                ELSE                               ! large-scale ascent
+                   tmp_tend = - w_subs(k) *                                    &
+                              ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+                ENDIF
+
+                tendency(k,j,i) = tendency(k,j,i) + tmp_tend
+
+                IF ( large_scale_forcing )  THEN
+                   sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend    &
+                                      * weight_pres(intermediate_timestep_count)
                 ENDIF
              ENDDO
@@ -190 +203 @@
 
        DO  k = nzb, nzt
-          IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
+          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
              var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
@@ -198 +211 @@
 !--   At the upper boundary, the initial profile is shifted with aid of 
 !--   the gradient tmp_grad. (This is ok if the gradients are linear.)
-      IF ( w_subs(nzt) < 0.0_wp ) THEN
+      IF ( w_subs(nzt) < 0.0_wp )  THEN
          tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
          var_mod(nzt+1) = var_init(nzt+1) -  &
@@ -206 +219 @@
 
       DO  k = nzt+1, nzb+1, -1
-         IF ( w_subs(k) >= 0.0_wp ) THEN  ! large-scale ascent
+         IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
             var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                ( var_init(k) - var_init(k-1) ) * ddzu(k) 
@@ -214 +227 @@
 !--   At the lower boundary shifting is not necessary because the 
 !--   subsidence velocity w_subs(nzb) vanishes.
-
-
-      IF ( w_subs(nzb+1) >= 0.0_wp ) THEN
+      IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
          var_mod(nzb) = var_init(nzb)
       ENDIF
@@ -225 +236 @@
  END SUBROUTINE subsidence
 
- SUBROUTINE subsidence_ij( i, j, tendency, var, var_init ) 
+ SUBROUTINE subsidence_ij( i, j, tendency, var, var_init, ls_index ) 
 
        USE arrays_3d,                                                          &
@@ -231 +242 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dt_3d
+           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing
 
        USE indices,                                                            &
@@ -237 +248 @@
 
        USE kinds
+
+       USE statistics,                                                         &
+           ONLY:  sums_ls_l, weight_pres
 
        IMPLICIT NONE
@@ -243 +257 @@
        INTEGER(iwp) ::  j !: 
        INTEGER(iwp) ::  k !: 
-
+       INTEGER(iwp) ::  ls_index !:
+
+       REAL(wp)     ::  tmp_tend !:
        REAL(wp)     ::  tmp_grad !: 
     
@@ -256 +272 @@
 !--    Influence of w_subsidence on the current tendency term
        DO  k = nzb_s_inner(j,i)+1, nzt 
-          IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
-             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
-                               ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
-          ELSE                             ! large-scale ascent
-             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
-                               ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
+             tmp_tend = - w_subs(k) * ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
+          ELSE                                 ! large-scale ascent
+             tmp_tend = - w_subs(k) * ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+          ENDIF
+
+          tendency(k,j,i) = tendency(k,j,i) + tmp_tend
+
+          IF ( large_scale_forcing )  THEN
+             sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend          &
+                                     * weight_pres(intermediate_timestep_count)
           ENDIF
        ENDDO
@@ -269 +290 @@
 !--    Shifting of the initial profile is especially necessary with Rayleigh 
 !--    damping switched on
-       IF ( i == nxl .AND. j == nys ) THEN ! shifting only once per PE 
+       IF ( i == nxl .AND. j == nys )  THEN ! shifting only once per PE 
 
           DO  k = nzb, nzt
-             IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
+             IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
                 var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                   ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
@@ -280 +301 @@
 !--       At the upper boundary, the initial profile is shifted with aid of 
 !--       the gradient tmp_grad. (This is ok if the gradients are linear.)
-          IF ( w_subs(nzt) < 0.0_wp ) THEN
+          IF ( w_subs(nzt) < 0.0_wp )  THEN
              tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
              var_mod(nzt+1) = var_init(nzt+1) -  &
@@ -288 +309 @@
 
           DO  k = nzt+1, nzb+1, -1
-             IF ( w_subs(k) >= 0.0_wp ) THEN  ! large-scale ascent
+             IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
                 var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                    ( var_init(k) - var_init(k-1) ) * ddzu(k)
@@ -296 +317 @@
 !--       At the lower boundary shifting is not necessary because the 
 !--       subsidence velocity w_subs(nzb) vanishes.
-
-
-          IF ( w_subs(nzb+1) >= 0.0_wp ) THEN
+          IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
              var_mod(nzb) = var_init(nzb)
           ENDIF

palm/trunk/SOURCE/time_integration.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+! Reset sums_ls_l to zero at each timestep
+! +sums_ls_l
+! Calculation of reference state (previously in subroutine calc_mean_profile)
+
 ! Former revisions:
 ! -----------------
@@ -126 +128 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  diss, e_p, nr_p, prho, pt, pt_p, ql, ql_c, ql_v, ql_vp, qr_p,   &
-               q_p, rho, sa_p, tend, u, u_p, v, vpt, v_p, w_p
-
-    USE buoyancy_mod,                                                          &
+        ONLY:  diss, e_p, nr_p, prho, pt, pt_p, q, ql, ql_c, ql_v, ql_vp, qr_p,&
+               q_p, ref_state, rho, sa_p, tend, u, u_p, v, vpt, v_p, w_p
+
+    USE calc_mean_profile_mod,                                                 &
         ONLY:  calc_mean_profile
 
@@ -149 +151 @@
                intermediate_timestep_count_max, large_scale_forcing,           &
                loop_optimization, lsf_surf, lsf_vert, masks, mid,              &
-               netcdf_data_format, neutral, nr_timesteps_this_run, ocean,      &
-               on_device, passive_scalar, prandtl_layer, precipitation,        &
+               netcdf_data_format, neutral, nr_timesteps_this_run, nudging,    &
+               ocean, on_device, passive_scalar, prandtl_layer, precipitation, &
                prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
                run_coupled, simulated_time, simulated_time_chr,                &
@@ -181 +183 @@
         ONLY:  ls_forcing_surf, ls_forcing_vert
 
+    USE nudge_mod,                                                             &
+        ONLY:  calc_tnudge
+
     USE particle_attributes,                                                   &
         ONLY:  particle_advection, particle_advection_start, wang_kernel
@@ -194 +199 @@
 
     USE statistics,                                                            &
-        ONLY:  flow_statistics_called, hom, pr_palm
+        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
 
     USE user_actions_mod,                                                      &
@@ -253 +258 @@
 !--    Determine ug, vg and w_subs in dependence on data from external file 
 !--    LSF_DATA
-       IF ( large_scale_forcing .AND. lsf_vert ) THEN
+       IF ( large_scale_forcing .AND. lsf_vert )  THEN
            CALL ls_forcing_vert ( simulated_time )
+           sums_ls_l = 0.0_wp
        ENDIF
 
@@ -283 +289 @@
 !--          buoyancy terms (WARNING: only the respective last call of
 !--          calc_mean_profile defines the reference state!)
-             IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4, 'time_int' )
-             IF ( ocean         )  CALL calc_mean_profile( rho, 64, 'time_int' )
-             IF ( humidity      )  CALL calc_mean_profile( vpt, 44, 'time_int' )
+             IF ( .NOT. neutral )  THEN
+                CALL calc_mean_profile( pt, 4 )
+                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
+             ENDIF
+             IF ( ocean )  THEN
+                CALL calc_mean_profile( rho, 64 )
+                ref_state(:)  = hom(:,1,64,0)
+             ENDIF
+             IF ( humidity )  THEN
+                CALL calc_mean_profile( vpt, 44 )
+                ref_state(:)  = hom(:,1,44,0)
+             ENDIF
+
           ENDIF
 
@@ -291 +307 @@
           IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
                intermediate_timestep_count == 1 )  CALL ws_statistics
+!
+!--       In case of nudging calculate current nudging time scale and horizontal
+!--       means of u,v,pt and q
+          IF ( nudging )  THEN
+             CALL calc_tnudge( simulated_time )
+             CALL calc_mean_profile( u, 1 )
+             CALL calc_mean_profile( v, 2 )
+             CALL calc_mean_profile( pt, 4 )
+             CALL calc_mean_profile( q, 41 )
+          ENDIF
 
 !