Changeset 1365


Ignore:
Timestamp:
Apr 22, 2014 3:03:56 PM (8 years ago)
Author:
boeske
Message:

large scale forcing enabled

Location:
palm/trunk/SOURCE
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1364 r1365  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Added new module calc_mean_profile, previously in module buoyancy,
     23# removed buoyancy dependency from nudging
    2324#
    2425# Former revisions:
     
    162163        advec_u_pw.f90 advec_u_up.f90 advec_v_pw.f90 advec_v_up.f90 \
    163164        advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
    164         buoyancy.f90 calc_liquid_water_content.f90 calc_precipitation.f90 \
    165         calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
    166         check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
    167         coriolis.f90 cpulog.f90 cuda_fft_interfaces.f90 data_log.f90 \
    168         data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
     165        buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
     166        calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
     167        check_for_restart.f90 check_open.f90 check_parameters.f90 \
     168        close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
     169        cuda_fft_interfaces.f90 data_log.f90 data_output_dvrp.f90 \
     170        data_output_mask.f90 data_output_profiles.f90 \
    169171        data_output_ptseries.f90 data_output_spectra.f90 \
    170172        data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \
     
    256258boundary_conds.o: modules.o mod_kinds.o
    257259buoyancy.o: modules.o mod_kinds.o
     260calc_mean_profile.o: modules.o mod_kinds.o
    258261calc_liquid_water_content.o: modules.o mod_kinds.o
    259262calc_precipitation.o: modules.o mod_kinds.o
     
    352355mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    353356netcdf.o: modules.o mod_kinds.o
    354 nudging.o: modules.o buoyancy.o cpulog.o mod_kinds.o
     357nudging.o: modules.o cpulog.o mod_kinds.o
    355358package_parin.o: modules.o mod_kinds.o
    356359palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
     
    382385swap_timelevel.o: modules.o cpulog.o mod_kinds.o
    383386temperton_fft.o: modules.o mod_kinds.o
    384 time_integration.o: modules.o advec_ws.o buoyancy.o cpulog.o interaction_droplets_ptq.o \
    385         ls_forcing.o mod_kinds.o production_e.o prognostic_equations.o user_actions.o
     387time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \
     388        cpulog.o interaction_droplets_ptq.o ls_forcing.o mod_kinds.o \
     389        production_e.o prognostic_equations.o user_actions.o
    386390time_to_string.o: mod_kinds.o
    387391timestep.o: modules.o cpulog.o mod_kinds.o
  • palm/trunk/SOURCE/buoyancy.f90

    r1354 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Calculation of reference state in subroutine calc_mean_profile moved to
     23! subroutine time_integration,
     24! subroutine calc_mean_profile moved to new file calc_mean_profile.f90
    2325!
    2426! Former revisions:
     
    7880
    7981    PRIVATE
    80     PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
     82    PUBLIC buoyancy, buoyancy_acc
    8183
    8284    INTERFACE buoyancy
     
    8890       MODULE PROCEDURE buoyancy_acc
    8991    END INTERFACE buoyancy_acc
    90 
    91     INTERFACE calc_mean_profile
    92        MODULE PROCEDURE calc_mean_profile
    93     END INTERFACE calc_mean_profile
    9492
    9593 CONTAINS
     
    376374    END SUBROUTINE buoyancy_ij
    377375
    378 
    379     SUBROUTINE calc_mean_profile( var, pr, loc )
    380 
    381 !------------------------------------------------------------------------------!
    382 ! Description:
    383 ! ------------
    384 ! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
    385 ! of potential temperature, 44 in case of virtual potential temperature, and 64
    386 ! in case of density (ocean runs)).
    387 !------------------------------------------------------------------------------!
    388 
    389        USE arrays_3d,                                                          &
    390            ONLY:  ref_state
    391 
    392        USE control_parameters,                                                 &
    393            ONLY:  intermediate_timestep_count, message_string
    394 
    395        USE indices,                                                            &
    396            ONLY:  ngp_2dh_s_inner, nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
    397 
    398        USE kinds
    399 
    400        USE pegrid
    401 
    402        USE statistics,                                                         &
    403            ONLY:  flow_statistics_called, hom, sums, sums_l
    404 
    405 
    406        IMPLICIT NONE
    407 
    408        CHARACTER (LEN=*) ::  loc !:
    409        
    410        INTEGER(iwp) ::  i                  !:
    411        INTEGER(iwp) ::  j                  !:
    412        INTEGER(iwp) ::  k                  !:
    413        INTEGER(iwp) ::  pr                 !:
    414        INTEGER(iwp) ::  omp_get_thread_num !:
    415        INTEGER(iwp) ::  tn                 !:
    416        
    417 #if defined( __nopointer )
    418        REAL(wp), DIMENSION(:,:,:) ::  var  !:
    419 #else
    420        REAL(wp), DIMENSION(:,:,:), POINTER ::  var
    421 #endif
    422 
    423 !
    424 !--    Computation of the horizontally averaged profile of variable var, unless
    425 !--    already done by the relevant call from flow_statistics. The calculation
    426 !--    is done only for the first respective intermediate timestep in order to
    427 !--    spare communication time and to produce identical model results with jobs
    428 !--    which are calling flow_statistics at different time intervals.
    429        IF ( .NOT. flow_statistics_called  .AND.                                &
    430             intermediate_timestep_count == 1 )  THEN
    431 
    432 !
    433 !--       Horizontal average of variable var
    434           tn           =   0  ! Default thread number in case of one thread
    435           !$OMP PARALLEL PRIVATE( i, j, k, tn )
    436 !$        tn = omp_get_thread_num()
    437           sums_l(:,pr,tn) = 0.0_wp
    438           !$OMP DO
    439           DO  i = nxl, nxr
    440              DO  j =  nys, nyn
    441                 DO  k = nzb_s_inner(j,i), nzt+1
    442                    sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
    443                 ENDDO
    444              ENDDO
    445           ENDDO
    446           !$OMP END PARALLEL
    447 
    448           DO  i = 1, threads_per_task-1
    449              sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
    450           ENDDO
    451 
    452 #if defined( __parallel )
    453 
    454           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    455           CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb,       &
    456                               MPI_REAL, MPI_SUM, comm2d, ierr )
    457 
    458 #else
    459 
    460           sums(:,pr) = sums_l(:,pr,0)
    461 
    462 #endif
    463 
    464           hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
    465 
    466        ENDIF
    467 
    468        SELECT CASE ( loc )
    469 
    470           CASE ( 'time_int' )
    471 
    472              ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
    473 
    474 
    475           CASE ( 'nudging' )
    476              !nothing to be done
    477 
    478 
    479           CASE DEFAULT
    480              message_string = 'unknown location "' // loc // '"'
    481              CALL message( 'calc_mean_profile', 'PA0379', 1, 2, 0, 6, 0 )
    482 
    483        END SELECT
    484 
    485 
    486 
    487     END SUBROUTINE calc_mean_profile
    488 
    489376 END MODULE buoyancy_mod
  • palm/trunk/SOURCE/check_parameters.f90

    r1362 r1365  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Usage of large scale forcing for pt and q enabled:
     23! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
     24! large scale subsidence (td_sub_lpt, td_sub_q)
     25! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
     26! check if use_subsidence_tendencies is used correctly
    2327!
    2428! Former revisions:
     
    12711275       ENDIF
    12721276         
     1277    ENDIF
     1278
     1279!
     1280!-- Check if the control parameter use_subsidence_tendencies is used correctly
     1281    IF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_subsidence )  THEN
     1282       message_string = 'The usage of use_subsidence_tendencies ' // &
     1283                            'requires large_scale_subsidence = .T..'
     1284       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
     1285    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
     1286       message_string = 'The usage of use_subsidence_tendencies ' // &
     1287                            'requires large_scale_forcing = .T..'
     1288       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
    12731289    ENDIF
    12741290
     
    26612677             ENDIF
    26622678
     2679          CASE ( 'td_lsa_lpt' )
     2680             IF ( .NOT. large_scale_forcing )  THEN
     2681                message_string = 'data_output_pr = ' // &
     2682                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2683                                 'lemented for large_scale_forcing = .FALSE.'
     2684                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
     2685             ELSE
     2686                dopr_index(i) = 81
     2687                dopr_unit(i)  = 'K/s'
     2688                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
     2689             ENDIF
     2690
     2691          CASE ( 'td_lsa_q' )
     2692             IF ( .NOT. large_scale_forcing )  THEN
     2693                message_string = 'data_output_pr = ' // &
     2694                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2695                                 'lemented for large_scale_forcing = .FALSE.'
     2696                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
     2697             ELSE
     2698                dopr_index(i) = 82
     2699                dopr_unit(i)  = 'kg/kgs'
     2700                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
     2701             ENDIF
     2702
     2703          CASE ( 'td_sub_lpt' )
     2704             IF ( .NOT. large_scale_forcing )  THEN
     2705                message_string = 'data_output_pr = ' // &
     2706                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2707                                 'lemented for large_scale_forcing = .FALSE.'
     2708                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
     2709             ELSE
     2710                dopr_index(i) = 83
     2711                dopr_unit(i)  = 'K/s'
     2712                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
     2713             ENDIF
     2714
     2715          CASE ( 'td_sub_q' )
     2716             IF ( .NOT. large_scale_forcing )  THEN
     2717                message_string = 'data_output_pr = ' // &
     2718                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2719                                 'lemented for large_scale_forcing = .FALSE.'
     2720                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
     2721             ELSE
     2722                dopr_index(i) = 84
     2723                dopr_unit(i)  = 'kg/kgs'
     2724                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
     2725             ENDIF
     2726
     2727          CASE ( 'td_nud_lpt' )
     2728             IF ( .NOT. nudging )  THEN
     2729                message_string = 'data_output_pr = ' // &
     2730                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2731                                 'lemented for nudging = .FALSE.'
     2732                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
     2733             ELSE
     2734                dopr_index(i) = 85
     2735                dopr_unit(i)  = 'K/s'
     2736                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
     2737             ENDIF
     2738
     2739          CASE ( 'td_nud_q' )
     2740             IF ( .NOT. nudging )  THEN
     2741                message_string = 'data_output_pr = ' // &
     2742                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2743                                 'lemented for nudging = .FALSE.'
     2744                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
     2745             ELSE
     2746                dopr_index(i) = 86
     2747                dopr_unit(i)  = 'kg/kgs'
     2748                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
     2749             ENDIF
     2750
     2751          CASE ( 'td_nud_u' )
     2752             IF ( .NOT. nudging )  THEN
     2753                message_string = 'data_output_pr = ' // &
     2754                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2755                                 'lemented for nudging = .FALSE.'
     2756                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
     2757             ELSE
     2758                dopr_index(i) = 87
     2759                dopr_unit(i)  = 'm/s2'
     2760                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
     2761             ENDIF
     2762
     2763          CASE ( 'td_nud_v' )
     2764             IF ( .NOT. nudging )  THEN
     2765                message_string = 'data_output_pr = ' // &
     2766                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2767                                 'lemented for nudging = .FALSE.'
     2768                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
     2769             ELSE
     2770                dopr_index(i) = 88
     2771                dopr_unit(i)  = 'm/s2'
     2772                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
     2773             ENDIF
     2774
     2775
    26632776          CASE DEFAULT
    26642777
  • palm/trunk/SOURCE/flow_statistics.f90

    r1354 r1365  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Output of large scale advection, large scale subsidence and nudging tendencies
     24! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
    2425!
    2526! Former revisions:
     
    9596
    9697    USE arrays_3d,                                                             &
    97         ONLY :  ddzu, ddzw, e, hyp, km, kh,nr,  p, prho, pt, q, qc, ql, qr,    &
    98                 qs, qsws, qswst, rho, sa, saswsb, saswst, shf, ts, tswst, u,   &
    99                 ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     98        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, pt_lsa, pt_subs, q,&
     99               qc, ql, qr, qs, qsws, qswst, q_lsa, q_subs, rho, sa, saswsb,    &
     100               saswst, shf, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,&
     101               v, vg, vpt, vswst, w, w_subs, zw
    100102       
    101103    USE cloud_parameters,                                                      &
     
    104106    USE control_parameters,                                                    &
    105107        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    106                 dt_3d, g, humidity, icloud_scheme, kappa, max_pr_user,         &
    107                 message_string, ocean, passive_scalar, precipitation,          &
    108                 use_surface_fluxes, use_top_fluxes, ws_scheme_mom, ws_scheme_sca
     108                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
     109                large_scale_subsidence, max_pr_user, message_string, ocean,    &
     110                passive_scalar, precipitation, simulated_time,                 &
     111                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
     112                ws_scheme_mom, ws_scheme_sca
    109113       
    110114    USE cpulog,                                                                &
     
    115119       
    116120    USE indices,                                                               &
    117         ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, nxl, &
    118                 nxr, nyn, nys, nzb, nzb_diff_s_inner, nzb_s_inner, nzt, nzt_diff
     121        ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
     122                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
     123                nzb_s_inner, nzt, nzt_diff
    119124       
    120125    USE kinds
     
    129134    INTEGER(iwp) ::  j                   !:
    130135    INTEGER(iwp) ::  k                   !:
     136    INTEGER(iwp) ::  nt                  !:
    131137    INTEGER(iwp) ::  omp_get_thread_num  !:
    132138    INTEGER(iwp) ::  sr                  !:
     
    136142   
    137143    REAL(wp) ::  dptdz_threshold  !:
     144    REAL(wp) ::  fac              !:
    138145    REAL(wp) ::  height           !:
    139146    REAL(wp) ::  pts              !:
     
    982989
    983990!
     991!--    Collect current large scale advection and subsidence tendencies for
     992!--    data output
     993       IF ( large_scale_forcing )  THEN
     994!
     995!--       Interpolation in time of LSF_DATA
     996          nt = 1
     997          DO WHILE ( simulated_time > time_vert(nt) )
     998             nt = nt + 1
     999          ENDDO
     1000          IF ( simulated_time /= time_vert(nt) )  THEN
     1001            nt = nt - 1
     1002          ENDIF
     1003
     1004          fac = ( simulated_time-time_vert(nt) )                               &
     1005                / ( time_vert(nt+1)-time_vert(nt) )
     1006
     1007
     1008          DO  k = nzb, nzt
     1009             sums_ls_l(k,0) = pt_lsa(k,nt)                                     &
     1010                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
     1011             sums_ls_l(k,1) = q_lsa(k,nt)                                      &
     1012                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
     1013          ENDDO
     1014
     1015          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
     1016
     1017             DO  k = nzb, nzt
     1018                sums_ls_l(k,2) = pt_subs(k,nt)                                 &
     1019                                 + fac * ( pt_subs(k,nt+1) - pt_subs(k,nt) )
     1020                sums_ls_l(k,3) = q_subs(k,nt)                                  &
     1021                                 + fac * ( q_subs(k,nt+1) - q_subs(k,nt) )
     1022             ENDDO
     1023
     1024          ENDIF
     1025
     1026       ENDIF
     1027
     1028!
    9841029!--    Calculate the user-defined profiles
    9851030       CALL user_statistics( 'profiles', sr, tn )
     
    10071052!--    Compute total sum from local sums
    10081053       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1009        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
     1054       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
    10101055                           MPI_SUM, comm2d, ierr )
     1056       IF ( large_scale_forcing )  THEN
     1057          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
     1058                              MPI_REAL, MPI_SUM, comm2d, ierr )
     1059       ENDIF
    10111060#else
    10121061       sums = sums_l(:,:,0)
     1062       IF ( large_scale_forcing )  THEN
     1063          sums(:,81:88) = sums_ls_l
     1064       ENDIF
    10131065#endif
    10141066
     
    10311083          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
    10321084          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
    1033           sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     1085          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
     1086          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
     1087          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    10341088       ENDDO
    10351089
     
    11301184       hom(:,1,80,sr) = w_subs         ! w_subs
    11311185
     1186       IF ( large_scale_forcing )  THEN
     1187          hom(:,1,81,sr) = sums_ls_l(:,0)          ! pt_lsa
     1188          hom(:,1,82,sr) = sums_ls_l(:,1)          ! q_lsa
     1189          IF ( use_subsidence_tendencies )  THEN
     1190             hom(:,1,83,sr) = sums_ls_l(:,2)       ! pt_subs
     1191             hom(:,1,84,sr) = sums_ls_l(:,3)       ! q_subs
     1192          ELSE
     1193             hom(:,1,83,sr) = sums(:,83)           ! pt_subs
     1194             hom(:,1,84,sr) = sums(:,84)           ! q_subs
     1195          ENDIF
     1196          hom(:,1,85,sr) = sums(:,85)              ! pt_nudge
     1197          hom(:,1,86,sr) = sums(:,86)              ! q_nudge
     1198          hom(:,1,87,sr) = sums(:,87)              ! u_nudge
     1199          hom(:,1,88,sr) = sums(:,88)              ! v_nudge
     1200       ENDIF
     1201
    11321202       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
    11331203                                       ! upstream-parts u_x, u_y, u_z, v_x,
     
    13091379
    13101380    USE arrays_3d,                                                             &
    1311         ONLY:  ddzu, ddzw, e, hyp, km, kh,nr,  p, prho, pt, q, qc, ql, qr,     &
    1312                qs, qsws, qswst, rho, sa, saswsb, saswst, shf, ts, tswst, u,    &
    1313                ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     1381        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, pt_lsa, pt_subs, q,&
     1382               qc, ql, qr, qs, qsws, qswst, q_lsa, q_subs, rho, sa, saswsb,    &
     1383               saswst, shf, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,&
     1384               v, vg, vpt, vswst, w, w_subs, zw
     1385                 
    13141386       
    13151387    USE cloud_parameters,                                                      &
     
    13171389       
    13181390    USE control_parameters,                                                    &
    1319         ONLY:  average_count_pr, cloud_droplets, cloud_physics, do_sum,        &
    1320                dt_3d, g, humidity, icloud_scheme, kappa, max_pr_user,          &
    1321                message_string, ocean, passive_scalar, precipitation,           &
    1322                use_surface_fluxes, use_top_fluxes, ws_scheme_mom, ws_scheme_sca
     1391        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
     1392                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
     1393                large_scale_subsidence, max_pr_user, message_string, ocean,    &
     1394                passive_scalar, precipitation, simulated_time,                 &
     1395                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
     1396                ws_scheme_mom, ws_scheme_sca
    13231397       
    13241398    USE cpulog,                                                                &
     
    13431417    INTEGER(iwp) ::  j                   !:
    13441418    INTEGER(iwp) ::  k                   !:
     1419    INTEGER(iwp) ::  nt                  !:
    13451420    INTEGER(iwp) ::  omp_get_thread_num  !:
    13461421    INTEGER(iwp) ::  sr                  !:
     
    13501425   
    13511426    REAL(wp) ::  dptdz_threshold  !:
     1427    REAL(wp) ::  fac              !:
    13521428    REAL(wp) ::  height           !:
    13531429    REAL(wp) ::  pts              !:
     
    26612737
    26622738!
     2739!--    Collect current large scale advection and subsidence tendencies for
     2740!--    data output
     2741       IF ( large_scale_forcing )  THEN
     2742!
     2743!--       Interpolation in time of LSF_DATA
     2744          nt = 1
     2745          DO WHILE ( simulated_time > time_vert(nt) )
     2746             nt = nt + 1
     2747          ENDDO
     2748          IF ( simulated_time /= time_vert(nt) )  THEN
     2749            nt = nt - 1
     2750          ENDIF
     2751
     2752          fac = ( simulated_time-time_vert(nt) )                               &
     2753                / ( time_vert(nt+1)-time_vert(nt) )
     2754
     2755
     2756          DO  k = nzb, nzt
     2757             sums_ls_l(k,0) = pt_lsa(k,nt)                                     &
     2758                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
     2759             sums_ls_l(k,1) = q_lsa(k,nt)                                      &
     2760                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
     2761          ENDDO
     2762
     2763          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
     2764
     2765             DO  k = nzb, nzt
     2766                sums_ls_l(k,2) = pt_subs(k,nt)                                 &
     2767                                 + fac * ( pt_subs(k,nt+1) - pt_subs(k,nt) )
     2768                sums_ls_l(k,3) = q_subs(k,nt)                                  &
     2769                                 + fac * ( q_subs(k,nt+1) - q_subs(k,nt) )
     2770             ENDDO
     2771
     2772          ENDIF
     2773
     2774       ENDIF
     2775
     2776!
    26632777!--    Calculate the user-defined profiles
    26642778       CALL user_statistics( 'profiles', sr, tn )
     
    26912805       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
    26922806                           MPI_SUM, comm2d, ierr )
     2807       IF ( large_scale_forcing )  THEN
     2808          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
     2809                              MPI_REAL, MPI_SUM, comm2d, ierr )
     2810       ENDIF
    26932811#else
    26942812       sums = sums_l(:,:,0)
     2813       IF ( large_scale_forcing )  THEN
     2814          sums(:,81:88) = sums_ls_l
     2815       ENDIF
    26952816#endif
    26962817
     
    27122833          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
    27132834          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
    2714           sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
    2715           sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     2835          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
     2836          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
     2837          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    27162838       ENDDO
    27172839
     
    28122934       hom(:,1,80,sr) = w_subs         ! w_subs
    28132935
     2936       IF ( large_scale_forcing )  THEN
     2937          hom(:,1,81,sr) = sums_ls_l(:,0)          ! pt_lsa
     2938          hom(:,1,82,sr) = sums_ls_l(:,1)          ! q_lsa
     2939          IF ( use_subsidence_tendencies )  THEN
     2940             hom(:,1,83,sr) = sums_ls_l(:,2)       ! pt_subs
     2941             hom(:,1,84,sr) = sums_ls_l(:,3)       ! q_subs
     2942          ELSE
     2943             hom(:,1,83,sr) = sums(:,83)           ! pt_subs
     2944             hom(:,1,84,sr) = sums(:,84)           ! q_subs
     2945          ENDIF
     2946          hom(:,1,85,sr) = sums(:,85)              ! pt_nudge
     2947          hom(:,1,86,sr) = sums(:,86)              ! q_nudge
     2948          hom(:,1,87,sr) = sums(:,87)              ! u_nudge
     2949          hom(:,1,88,sr) = sums(:,88)              ! v_nudge
     2950       END IF
     2951
    28142952       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
    28152953                                       ! upstream-parts u_x, u_y, u_z, v_x,
  • palm/trunk/SOURCE/header.f90

    r1360 r1365  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New section 'Large scale forcing and nudging':
     23! output of large scale forcing and nudging information,
     24! new section for initial profiles created
    2325!
    2426! Former revisions:
     
    98100!
    99101! 1015 2012-09-27 09:23:24Z raasch
    100 ! output of Aajustment of mixing length to the Prandtl mixing length at first
     102! output of Adjustment of mixing length to the Prandtl mixing length at first
    101103! grid point above ground removed
    102104!
     
    452454       ENDIF
    453455    ENDIF
    454     IF ( large_scale_subsidence )  THEN
    455         WRITE ( io, 153 )
    456         WRITE ( io, 154 )
    457     ENDIF
    458     IF ( nudging )  THEN
    459         WRITE ( io, 155 )
    460         WRITE ( io, 156 )
    461     ENDIF
    462     IF ( large_scale_forcing )  THEN
    463         WRITE ( io, 157 )
    464         WRITE ( io, 158 )
    465     ENDIF
    466456    WRITE ( io, 99 )
    467457
     
    535525
    536526!
     527!-- Large scale forcing and nudging
     528    WRITE ( io, 160 )
     529    IF ( large_scale_forcing )  THEN
     530       WRITE ( io, 162 )
     531       WRITE ( io, 163 )
     532
     533       IF ( large_scale_subsidence )  THEN
     534          IF ( .NOT. use_subsidence_tendencies )  THEN
     535             WRITE ( io, 164 )
     536          ELSE
     537             WRITE ( io, 165 )
     538          ENDIF
     539       ENDIF
     540
     541       IF ( bc_pt_b == 'dirichlet' )  THEN
     542          WRITE ( io, 180 )
     543       ELSEIF ( bc_pt_b == 'neumann' )  THEN
     544          WRITE ( io, 181 )
     545       ENDIF
     546
     547       IF ( bc_q_b == 'dirichlet' )  THEN
     548          WRITE ( io, 182 )
     549       ELSEIF ( bc_q_b == 'neumann' )  THEN
     550          WRITE ( io, 183 )
     551       ENDIF
     552
     553       WRITE ( io, 167 )
     554       IF ( nudging )  THEN
     555          WRITE ( io, 170 )
     556       ENDIF
     557    ELSE
     558       WRITE ( io, 161 )
     559       WRITE ( io, 171 )
     560    ENDIF
     561    IF ( large_scale_subsidence )  THEN
     562       WRITE ( io, 168 )
     563       WRITE ( io, 169 )
     564    ENDIF
     565
     566!
     567!-- Profile for the large scale vertial velocity
     568!-- Building output strings, starting with surface value
     569    IF ( large_scale_subsidence )  THEN
     570       temperatures = '   0.0'
     571       gradients = '------'
     572       slices = '     0'
     573       coordinates = '   0.0'
     574       i = 1
     575       DO  WHILE ( subs_vertical_gradient_level_i(i) /= -9999 )
     576
     577          WRITE (coor_chr,'(E10.2,7X)')  &
     578                                w_subs(subs_vertical_gradient_level_i(i))
     579          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
     580
     581          WRITE (coor_chr,'(E10.2,7X)')  subs_vertical_gradient(i)
     582          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     583
     584          WRITE (coor_chr,'(I10,7X)')  subs_vertical_gradient_level_i(i)
     585          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     586
     587          WRITE (coor_chr,'(F10.2,7X)')  subs_vertical_gradient_level(i)
     588          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     589
     590          IF ( i == 10 )  THEN
     591             EXIT
     592          ELSE
     593             i = i + 1
     594          ENDIF
     595
     596       ENDDO
     597
     598 
     599       IF ( .NOT. large_scale_forcing )  THEN
     600          WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
     601                             TRIM( gradients ), TRIM( slices )
     602       ENDIF
     603
     604
     605    ENDIF
     606
     607!-- Profile of the geostrophic wind (component ug)
     608!-- Building output strings
     609    WRITE ( ugcomponent, '(F6.2)' )  ug_surface
     610    gradients = '------'
     611    slices = '     0'
     612    coordinates = '   0.0'
     613    i = 1
     614    DO  WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 )
     615     
     616       WRITE (coor_chr,'(F6.2,1X)')  ug(ug_vertical_gradient_level_ind(i))
     617       ugcomponent = TRIM( ugcomponent ) // '  ' // TRIM( coor_chr )
     618
     619       WRITE (coor_chr,'(F6.2,1X)')  ug_vertical_gradient(i)
     620       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
     621
     622       WRITE (coor_chr,'(I6,1X)')  ug_vertical_gradient_level_ind(i)
     623       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
     624
     625       WRITE (coor_chr,'(F6.1,1X)')  ug_vertical_gradient_level(i)
     626       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
     627
     628       IF ( i == 10 )  THEN
     629          EXIT
     630       ELSE
     631          i = i + 1
     632       ENDIF
     633
     634    ENDDO
     635
     636    IF ( .NOT. large_scale_forcing )  THEN
     637       WRITE ( io, 423 )  TRIM( coordinates ), TRIM( ugcomponent ), &
     638                          TRIM( gradients ), TRIM( slices )
     639    ENDIF
     640
     641!-- Profile of the geostrophic wind (component vg)
     642!-- Building output strings
     643    WRITE ( vgcomponent, '(F6.2)' )  vg_surface
     644    gradients = '------'
     645    slices = '     0'
     646    coordinates = '   0.0'
     647    i = 1
     648    DO  WHILE ( vg_vertical_gradient_level_ind(i) /= -9999 )
     649
     650       WRITE (coor_chr,'(F6.2,1X)')  vg(vg_vertical_gradient_level_ind(i))
     651       vgcomponent = TRIM( vgcomponent ) // '  ' // TRIM( coor_chr )
     652
     653       WRITE (coor_chr,'(F6.2,1X)')  vg_vertical_gradient(i)
     654       gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
     655
     656       WRITE (coor_chr,'(I6,1X)')  vg_vertical_gradient_level_ind(i)
     657       slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
     658
     659       WRITE (coor_chr,'(F6.1,1X)')  vg_vertical_gradient_level(i)
     660       coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
     661
     662       IF ( i == 10 )  THEN
     663          EXIT
     664       ELSE
     665          i = i + 1
     666       ENDIF
     667 
     668    ENDDO
     669
     670    IF ( .NOT. large_scale_forcing )  THEN
     671       WRITE ( io, 424 )  TRIM( coordinates ), TRIM( vgcomponent ), &
     672                          TRIM( gradients ), TRIM( slices )
     673    ENDIF
     674
     675!
    537676!-- Topography
    538677    WRITE ( io, 270 )  topography
     
    805944       ENDIF
    806945    ENDIF
     946
     947!
     948!-- Initial Profiles
     949    WRITE ( io, 321 )
     950!
     951!-- Initial wind profiles
     952    IF ( u_profile(1) /= 9999999.9_wp )  WRITE ( io, 427 )
     953
     954!
     955!-- Initial temperature profile
     956!-- Building output strings, starting with surface temperature
     957    WRITE ( temperatures, '(F6.2)' )  pt_surface
     958    gradients = '------'
     959    slices = '     0'
     960    coordinates = '   0.0'
     961    i = 1
     962    DO  WHILE ( pt_vertical_gradient_level_ind(i) /= -9999 )
     963
     964       WRITE (coor_chr,'(F7.2)')  pt_init(pt_vertical_gradient_level_ind(i))
     965       temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
     966
     967       WRITE (coor_chr,'(F7.2)')  pt_vertical_gradient(i)
     968       gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     969
     970       WRITE (coor_chr,'(I7)')  pt_vertical_gradient_level_ind(i)
     971       slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     972
     973       WRITE (coor_chr,'(F7.1)')  pt_vertical_gradient_level(i)
     974       coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     975
     976       IF ( i == 10 )  THEN
     977          EXIT
     978       ELSE
     979          i = i + 1
     980       ENDIF
     981
     982    ENDDO
     983
     984    IF ( .NOT. nudging )  THEN
     985       WRITE ( io, 420 )  TRIM( coordinates ), TRIM( temperatures ), &
     986                          TRIM( gradients ), TRIM( slices )
     987    ELSE
     988       WRITE ( io, 428 )
     989    ENDIF
     990
     991!
     992!-- Initial humidity profile
     993!-- Building output strings, starting with surface humidity
     994    IF ( humidity  .OR.  passive_scalar )  THEN
     995       WRITE ( temperatures, '(E8.1)' )  q_surface
     996       gradients = '--------'
     997       slices = '       0'
     998       coordinates = '     0.0'
     999       i = 1
     1000       DO  WHILE ( q_vertical_gradient_level_ind(i) /= -9999 )
     1001         
     1002          WRITE (coor_chr,'(E8.1,4X)')  q_init(q_vertical_gradient_level_ind(i))
     1003          temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
     1004
     1005          WRITE (coor_chr,'(E8.1,4X)')  q_vertical_gradient(i)
     1006          gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
     1007         
     1008          WRITE (coor_chr,'(I8,4X)')  q_vertical_gradient_level_ind(i)
     1009          slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
     1010         
     1011          WRITE (coor_chr,'(F8.1,4X)')  q_vertical_gradient_level(i)
     1012          coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
     1013
     1014          IF ( i == 10 )  THEN
     1015             EXIT
     1016          ELSE
     1017             i = i + 1
     1018          ENDIF
     1019
     1020       ENDDO
     1021
     1022       IF ( humidity )  THEN
     1023          IF ( .NOT. nudging )  THEN
     1024             WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
     1025                                TRIM( gradients ), TRIM( slices )
     1026          ENDIF
     1027       ELSE
     1028          WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
     1029                             TRIM( gradients ), TRIM( slices )
     1030       ENDIF
     1031    ENDIF
     1032
     1033!
     1034!-- Initial salinity profile
     1035!-- Building output strings, starting with surface salinity
     1036    IF ( ocean )  THEN
     1037       WRITE ( temperatures, '(F6.2)' )  sa_surface
     1038       gradients = '------'
     1039       slices = '     0'
     1040       coordinates = '   0.0'
     1041       i = 1
     1042       DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
     1043
     1044          WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
     1045          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
     1046
     1047          WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
     1048          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     1049
     1050          WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
     1051          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     1052
     1053          WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
     1054          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     1055
     1056          IF ( i == 10 )  THEN
     1057             EXIT
     1058          ELSE
     1059             i = i + 1
     1060          ENDIF
     1061
     1062       ENDDO
     1063
     1064       WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
     1065                          TRIM( gradients ), TRIM( slices )
     1066    ENDIF
     1067
    8071068
    8081069!
     
    12681529          IF ( precipitation )  WRITE ( io, 511 ) c_sedimentation
    12691530       ENDIF
    1270     ENDIF
    1271 
    1272 !-- Profile of the geostrophic wind (component ug)
    1273 !-- Building output strings
    1274     WRITE ( ugcomponent, '(F6.2)' )  ug_surface
    1275     gradients = '------'
    1276     slices = '     0'
    1277     coordinates = '   0.0'
    1278     i = 1
    1279     DO  WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 )
    1280      
    1281        WRITE (coor_chr,'(F6.2,1X)')  ug(ug_vertical_gradient_level_ind(i))
    1282        ugcomponent = TRIM( ugcomponent ) // '  ' // TRIM( coor_chr )
    1283 
    1284        WRITE (coor_chr,'(F6.2,1X)')  ug_vertical_gradient(i)
    1285        gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
    1286 
    1287        WRITE (coor_chr,'(I6,1X)')  ug_vertical_gradient_level_ind(i)
    1288        slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
    1289 
    1290        WRITE (coor_chr,'(F6.1,1X)')  ug_vertical_gradient_level(i)
    1291        coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
    1292 
    1293        IF ( i == 10 )  THEN
    1294           EXIT
    1295        ELSE
    1296           i = i + 1
    1297        ENDIF
    1298 
    1299     ENDDO
    1300 
    1301     IF ( .NOT. large_scale_forcing )  THEN
    1302        WRITE ( io, 423 )  TRIM( coordinates ), TRIM( ugcomponent ), &
    1303                           TRIM( gradients ), TRIM( slices )
    1304     ELSE
    1305        WRITE ( io, 429 )
    1306     ENDIF
    1307 
    1308 !-- Profile of the geostrophic wind (component vg)
    1309 !-- Building output strings
    1310     WRITE ( vgcomponent, '(F6.2)' )  vg_surface
    1311     gradients = '------'
    1312     slices = '     0'
    1313     coordinates = '   0.0'
    1314     i = 1
    1315     DO  WHILE ( vg_vertical_gradient_level_ind(i) /= -9999 )
    1316 
    1317        WRITE (coor_chr,'(F6.2,1X)')  vg(vg_vertical_gradient_level_ind(i))
    1318        vgcomponent = TRIM( vgcomponent ) // '  ' // TRIM( coor_chr )
    1319 
    1320        WRITE (coor_chr,'(F6.2,1X)')  vg_vertical_gradient(i)
    1321        gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
    1322 
    1323        WRITE (coor_chr,'(I6,1X)')  vg_vertical_gradient_level_ind(i)
    1324        slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
    1325 
    1326        WRITE (coor_chr,'(F6.1,1X)')  vg_vertical_gradient_level(i)
    1327        coordinates = TRIM( coordinates ) // '  ' // TRIM( coor_chr )
    1328 
    1329        IF ( i == 10 )  THEN
    1330           EXIT
    1331        ELSE
    1332           i = i + 1
    1333        ENDIF
    1334  
    1335     ENDDO
    1336 
    1337     IF ( .NOT. large_scale_forcing )  THEN
    1338        WRITE ( io, 424 )  TRIM( coordinates ), TRIM( vgcomponent ), &
    1339                           TRIM( gradients ), TRIM( slices )
    1340     ENDIF
    1341 
    1342 !
    1343 !-- Initial wind profiles
    1344     IF ( u_profile(1) /= 9999999.9_wp )  WRITE ( io, 427 )
    1345 
    1346 !
    1347 !-- Initial temperature profile
    1348 !-- Building output strings, starting with surface temperature
    1349     WRITE ( temperatures, '(F6.2)' )  pt_surface
    1350     gradients = '------'
    1351     slices = '     0'
    1352     coordinates = '   0.0'
    1353     i = 1
    1354     DO  WHILE ( pt_vertical_gradient_level_ind(i) /= -9999 )
    1355 
    1356        WRITE (coor_chr,'(F7.2)')  pt_init(pt_vertical_gradient_level_ind(i))
    1357        temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
    1358 
    1359        WRITE (coor_chr,'(F7.2)')  pt_vertical_gradient(i)
    1360        gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
    1361 
    1362        WRITE (coor_chr,'(I7)')  pt_vertical_gradient_level_ind(i)
    1363        slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
    1364 
    1365        WRITE (coor_chr,'(F7.1)')  pt_vertical_gradient_level(i)
    1366        coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
    1367 
    1368        IF ( i == 10 )  THEN
    1369           EXIT
    1370        ELSE
    1371           i = i + 1
    1372        ENDIF
    1373 
    1374     ENDDO
    1375 
    1376     IF ( .NOT. nudging )  THEN
    1377        WRITE ( io, 420 )  TRIM( coordinates ), TRIM( temperatures ), &
    1378                           TRIM( gradients ), TRIM( slices )
    1379     ELSE
    1380        WRITE ( io, 428 )
    1381     ENDIF
    1382 
    1383 !
    1384 !-- Initial humidity profile
    1385 !-- Building output strings, starting with surface humidity
    1386     IF ( humidity  .OR.  passive_scalar )  THEN
    1387        WRITE ( temperatures, '(E8.1)' )  q_surface
    1388        gradients = '--------'
    1389        slices = '       0'
    1390        coordinates = '     0.0'
    1391        i = 1
    1392        DO  WHILE ( q_vertical_gradient_level_ind(i) /= -9999 )
    1393          
    1394           WRITE (coor_chr,'(E8.1,4X)')  q_init(q_vertical_gradient_level_ind(i))
    1395           temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
    1396 
    1397           WRITE (coor_chr,'(E8.1,4X)')  q_vertical_gradient(i)
    1398           gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
    1399          
    1400           WRITE (coor_chr,'(I8,4X)')  q_vertical_gradient_level_ind(i)
    1401           slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
    1402          
    1403           WRITE (coor_chr,'(F8.1,4X)')  q_vertical_gradient_level(i)
    1404           coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
    1405 
    1406           IF ( i == 10 )  THEN
    1407              EXIT
    1408           ELSE
    1409              i = i + 1
    1410           ENDIF
    1411 
    1412        ENDDO
    1413 
    1414        IF ( humidity )  THEN
    1415           IF ( .NOT. nudging )  THEN
    1416              WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
    1417                                 TRIM( gradients ), TRIM( slices )
    1418           ENDIF
    1419        ELSE
    1420           WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
    1421                              TRIM( gradients ), TRIM( slices )
    1422        ENDIF
    1423     ENDIF
    1424 
    1425 !
    1426 !-- Initial salinity profile
    1427 !-- Building output strings, starting with surface salinity
    1428     IF ( ocean )  THEN
    1429        WRITE ( temperatures, '(F6.2)' )  sa_surface
    1430        gradients = '------'
    1431        slices = '     0'
    1432        coordinates = '   0.0'
    1433        i = 1
    1434        DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
    1435 
    1436           WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
    1437           temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
    1438 
    1439           WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
    1440           gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
    1441 
    1442           WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
    1443           slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
    1444 
    1445           WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
    1446           coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
    1447 
    1448           IF ( i == 10 )  THEN
    1449              EXIT
    1450           ELSE
    1451              i = i + 1
    1452           ENDIF
    1453 
    1454        ENDDO
    1455 
    1456        WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
    1457                           TRIM( gradients ), TRIM( slices )
    1458     ENDIF
    1459 
    1460 !
    1461 !-- Profile for the large scale vertial velocity
    1462 !-- Building output strings, starting with surface value
    1463     IF ( large_scale_subsidence )  THEN
    1464        temperatures = '   0.0'
    1465        gradients = '------'
    1466        slices = '     0'
    1467        coordinates = '   0.0'
    1468        i = 1
    1469        DO  WHILE ( subs_vertical_gradient_level_i(i) /= -9999 )
    1470 
    1471           WRITE (coor_chr,'(E10.2,7X)')  &
    1472                                 w_subs(subs_vertical_gradient_level_i(i))
    1473           temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
    1474 
    1475           WRITE (coor_chr,'(E10.2,7X)')  subs_vertical_gradient(i)
    1476           gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
    1477 
    1478           WRITE (coor_chr,'(I10,7X)')  subs_vertical_gradient_level_i(i)
    1479           slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
    1480 
    1481           WRITE (coor_chr,'(F10.2,7X)')  subs_vertical_gradient_level(i)
    1482           coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
    1483 
    1484           IF ( i == 10 )  THEN
    1485              EXIT
    1486           ELSE
    1487              i = i + 1
    1488           ENDIF
    1489 
    1490        ENDDO
    1491 
    1492  
    1493        IF ( .NOT. large_scale_forcing )  THEN
    1494           WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
    1495                              TRIM( gradients ), TRIM( slices )
    1496         ELSE
    1497            WRITE ( io, 460 )
    1498        ENDIF
    1499 
    1500 
    15011531    ENDIF
    15021532
     
    17271757           /'     ',2(1X,E12.5),'Pa/m in x/y direction', &
    17281758           /'     starting from dp_level_b =', F8.3, 'm', A /)
    1729 153 FORMAT (' --> Large-scale vertical motion is used in the ', &
     1759160 FORMAT (//' Large scale forcing and nudging:'/ &
     1760              ' -------------------------------'/)
     1761161 FORMAT (' --> No large scale forcing from external is used (default) ')
     1762162 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
     1763163 FORMAT ('     - large scale advection tendencies ')
     1764164 FORMAT ('     - large scale subsidence velocity w_subs ')
     1765165 FORMAT ('     - large scale subsidence tendencies ')
     1766167 FORMAT ('     - and geostrophic wind components ug and vg')
     1767168 FORMAT (' --> Large-scale vertical motion is used in the ', &
    17301768                  'prognostic equation(s) for')
    1731 154 FORMAT ('     the scalar(s) only')
    1732 155 FORMAT (' --> Nudging is used - initial profiles for u, v, pt and q ',&
    1733                   'correspond to the')
    1734 156 FORMAT ('     first profiles in NUDGING_DATA')
    1735 157 FORMAT (' --> Large scale forcing from external file (LSF_DATA) is used: ')
    1736 158 FORMAT ('     prescribed surfaces fluxes and geostrophic wind components')
     1769169 FORMAT ('     the scalar(s) only')
     1770170 FORMAT (' --> Nudging is used')
     1771171 FORMAT (' --> No nudging is used (default) ')
     1772180 FORMAT ('     - prescribed surface values for temperature')
     1773181 FORMAT ('     - prescribed surface fluxes for humidity')
     1774182 FORMAT ('     - prescribed surface values for temperature')
     1775183 FORMAT ('     - prescribed surface fluxes for humidity')
    17371776200 FORMAT (//' Run time and time step information:'/ &
    17381777             ' ----------------------------------'/)
     
    18321871320 FORMAT ('       Predefined constant momentumflux:  u: ',F9.6,' m**2/s**2'/ &
    18331872            '                                          v: ',F9.6,' m**2/s**2')
     1873321 FORMAT (//' Initial profiles:'/ &
     1874              ' ----------------')
    18341875325 FORMAT (//' List output:'/ &
    18351876             ' -----------'//  &
     
    19892030428 FORMAT (/'    Initial profiles (u, v, pt, q) are taken from file '/ &
    19902031             '    NUDGING_DATA')
    1991 429 FORMAT (/'    Geostrophic wind profiles (ug, vg) are are taken from file '/ &
    1992              '    LSF_DATA')
    19932032430 FORMAT (//' Cloud physics quantities / methods:'/ &
    19942033              ' ----------------------------------'/)
     
    20172056454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
    20182057455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
    2019 460 FORMAT (/'    Profiles for large scale vertical velocity are '/ &
    2020              '    taken from file LSF_DATA')
    20212058470 FORMAT (//' Actions during the simulation:'/ &
    20222059              ' -----------------------------'/)
  • palm/trunk/SOURCE/ls_forcing.f90

    r1354 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Usage of large scale forcing for pt and q enabled:
     23! Added new subroutine ls_advec for horizontal large scale advection and large
     24! scale subsidence,
     25! error message in init_ls_forcing specified,
     26! variable t renamed nt
    2327!
    2428! Former revisions:
     
    6165
    6266    PRIVATE
    63     PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert 
     67    PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec
    6468    SAVE
    6569
     70    INTERFACE ls_advec
     71       MODULE PROCEDURE ls_advec
     72       MODULE PROCEDURE ls_advec_ij
     73    END INTERFACE ls_advec
    6674
    6775 CONTAINS
     
    7078
    7179       USE arrays_3d,                                                          &
    72            ONLY:  p_surf, pt_surf, q_surf, qsws_surf, shf_surf, time_surf,     &
    73                   time_vert, ug_vert, vg_vert, wsubs_vert, zu
     80           ONLY:  p_surf, pt_lsa, pt_subs, pt_surf, q_lsa, q_subs, q_surf,     &
     81                  qsws_surf, shf_surf, time_surf, time_vert, ug_vert, vg_vert, &
     82                  wsubs_vert, zu
    7483
    7584       USE control_parameters,                                                 &
     
    7786
    7887       USE indices,                                                            &
    79            ONLY:  nzb, nzt
     88           ONLY:  ngp_sums_ls, nzb, nz, nzt
    8089
    8190       USE kinds
    8291
     92       USE statistics,                                                         &
     93           ONLY:  sums_ls_l
     94
    8395
    8496       IMPLICIT NONE
    8597
    86        CHARACTER(100) ::  chmess              !:
    87        CHARACTER(1)   ::  hash                !:
    88 
    89        INTEGER(iwp) ::  ierrn                 !:
    90        INTEGER(iwp) ::  finput = 90           !:
    91        INTEGER(iwp) ::  k                     !:
    92        INTEGER(iwp) ::  t                     !:
    93 
    94        REAL(wp) ::  fac                       !:
    95        REAL(wp) ::  highheight                !:
    96        REAL(wp) ::  highug_vert               !:
    97        REAL(wp) ::  highvg_vert               !:
    98        REAL(wp) ::  highwsubs_vert            !:
    99        REAL(wp) ::  lowheight                 !:
    100        REAL(wp) ::  lowug_vert                !:
    101        REAL(wp) ::  lowvg_vert                !:
    102        REAL(wp) ::  lowwsubs_vert             !:
    103        REAL(wp) ::  r_dummy                   !:
    104 
    105        ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),         &
    106                  qsws_surf(0:nlsf), shf_surf(0:nlsf), time_vert(0:nlsf),  &
    107                  time_surf(0:nlsf), ug_vert(nzb:nzt+1,0:nlsf),            &
    108                  vg_vert(nzb:nzt+1,0:nlsf), wsubs_vert(nzb:nzt+1,0:nlsf) )
    109 
    110        p_surf     = 0.0_wp;  pt_surf  = 0.0_wp;  q_surf    = 0.0_wp
    111        qsws_surf  = 0.0_wp;  shf_surf = 0.0_wp;  time_vert = 0.0_wp
    112        time_surf  = 0.0_wp;  ug_vert  = 0.0_wp;  vg_vert   = 0.0_wp
    113        wsubs_vert = 0.0_wp
    114 
     98       CHARACTER(100) ::  chmess      !:
     99       CHARACTER(1)   ::  hash        !:
     100
     101       INTEGER(iwp) ::  ierrn         !:
     102       INTEGER(iwp) ::  finput = 90   !:
     103       INTEGER(iwp) ::  k             !:
     104       INTEGER(iwp) ::  nt             !:
     105
     106       REAL(wp) ::  fac               !:
     107       REAL(wp) ::  highheight        !:
     108       REAL(wp) ::  highug_vert       !:
     109       REAL(wp) ::  highvg_vert       !:
     110       REAL(wp) ::  highwsubs_vert    !:
     111       REAL(wp) ::  lowheight         !:
     112       REAL(wp) ::  lowug_vert        !:
     113       REAL(wp) ::  lowvg_vert        !:
     114       REAL(wp) ::  lowwsubs_vert     !:
     115       REAL(wp) ::  highpt_lsa        !:
     116       REAL(wp) ::  lowpt_lsa         !:
     117       REAL(wp) ::  highq_lsa         !:
     118       REAL(wp) ::  lowq_lsa          !:
     119       REAL(wp) ::  highpt_subs       !:
     120       REAL(wp) ::  lowpt_subs        !:
     121       REAL(wp) ::  highq_subs        !:
     122       REAL(wp) ::  lowq_subs         !:
     123       REAL(wp) ::  r_dummy           !:
     124
     125       ALLOCATE( p_surf(0:nlsf), pt_lsa(nzb:nzt+1,0:nlsf),                     &
     126                 pt_subs(nzb:nzt+1,0:nlsf), pt_surf(0:nlsf),                   &
     127                 q_lsa(nzb:nzt+1,0:nlsf), q_subs(nzb:nzt+1,0:nlsf),            &
     128                 q_surf(0:nlsf), qsws_surf(0:nlsf), shf_surf(0:nlsf),          &
     129                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
     130                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
     131                 wsubs_vert(nzb:nzt+1,0:nlsf) )
     132
     133       p_surf = 0.0_wp; pt_lsa = 0.0; pt_subs = 0.0; pt_surf = 0.0_wp
     134       q_lsa = 0.0; q_subs = 0.0; q_surf = 0.0_wp; qsws_surf = 0.0_wp
     135       shf_surf = 0.0_wp; time_vert = 0.0_wp; time_surf = 0.0_wp;
     136       ug_vert = 0.0_wp; vg_vert = 0.0_wp; wsubs_vert = 0.0_wp
     137
     138!
     139!--    Array for storing large scale forcing and nudging tendencies at each
     140!--    timestep for data output
     141       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
     142       sums_ls_l = 0.0_wp
     143
     144       ngp_sums_ls = (nz+2)*6
    115145
    116146       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
     
    136166!
    137167!--    Surface values are read in
    138        t     = 0
     168       nt     = 0
    139169       ierrn = 0
    140170
    141        DO WHILE ( time_surf(t) < end_time )
    142           t = t + 1
    143           READ ( finput, *, IOSTAT = ierrn ) time_surf(t), shf_surf(t), &
    144                                              qsws_surf(t), pt_surf(t),  &
    145                                              q_surf(t), p_surf(t)
     171       DO WHILE ( time_surf(nt) < end_time )
     172          nt = nt + 1
     173          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
     174                                             qsws_surf(nt), pt_surf(nt),       &
     175                                             q_surf(nt), p_surf(nt)
    146176
    147177          IF ( ierrn < 0 )  THEN
     
    156186       IF ( time_surf(1) > end_time )  THEN
    157187          WRITE ( message_string, * ) 'No time dependent surface variables in ',&
    158                                      '&LSF_DATA for end of run found - ',  &
     188                                     '&LSF_DATA for end of run found - ',      &
    159189                                     'lsf_surf is set to FALSE'
    160190          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
     
    171201!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
    172202
    173        t = 0
    174        DO WHILE ( time_vert(t) < end_time )
    175           t = t + 1
     203       nt = 0
     204       DO WHILE ( time_vert(nt) < end_time )
     205          nt = nt + 1
    176206          hash = "#"
    177207          ierrn = 1 ! not zero
     
    180210!--       from there onwards the profiles will be read
    181211          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
    182              READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(t)
     212             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
    183213             IF ( ierrn < 0 )  THEN
    184214                WRITE( message_string, * ) 'No time dependent vertical profiles',&
     
    188218          ENDDO
    189219
    190           IF ( t == 1 .AND. time_vert(t) > end_time ) EXIT
    191 
    192           READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, &
    193                                            lowwsubs_vert
     220          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
     221
     222          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
     223                                           lowwsubs_vert, lowpt_lsa,           &
     224                                           lowq_lsa, lowpt_subs, lowq_subs
    194225          IF ( ierrn /= 0 )  THEN
    195226             message_string = 'errors in file LSF_DATA'
     
    197228          ENDIF
    198229
    199           READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, highvg_vert, &
    200                                            highwsubs_vert
     230          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
     231                                           highvg_vert, highwsubs_vert,        &
     232                                           highpt_lsa, highq_lsa,              &
     233                                           highpt_subs, highq_subs
    201234     
    202235          IF ( ierrn /= 0 )  THEN
     
    212245                lowvg_vert    = highvg_vert
    213246                lowwsubs_vert = highwsubs_vert
     247                lowpt_lsa     = highpt_lsa
     248                lowq_lsa      = highq_lsa
     249                lowpt_subs    = highpt_subs
     250                lowq_subs     = highq_subs
    214251
    215252                ierrn = 0
    216                 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
    217                                                  highvg_vert, highwsubs_vert
     253                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
     254                                                 highvg_vert, highwsubs_vert,  &
     255                                                 highpt_lsa, highq_lsa,        &
     256                                                 highpt_subs, highq_subs
    218257
    219258                IF ( ierrn /= 0 )  THEN
    220                    message_string = 'errors in file LSF_DATA'
    221                    CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     259                   WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm ', &
     260                        'is higher than the maximum height in LSF_DATA which ',&
     261                        'is ', lowheight, 'm. Interpolation on PALM ',         &
     262                        'grid is not possible.'
     263                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
    222264                ENDIF
    223265
     
    228270             fac = (highheight-zu(k))/(highheight - lowheight)
    229271
    230              ug_vert(k,t)    = fac*lowug_vert    + (1-fac)*highug_vert
    231              vg_vert(k,t)    = fac*lowvg_vert    + (1-fac)*highvg_vert
    232              wsubs_vert(k,t) = fac*lowwsubs_vert + (1-fac)*highwsubs_vert
     272             ug_vert(k,nt)    = fac*lowug_vert    + (1.0_wp-fac)*highug_vert
     273             vg_vert(k,nt)    = fac*lowvg_vert    + (1.0_wp-fac)*highvg_vert
     274             wsubs_vert(k,nt) = fac*lowwsubs_vert + (1.0_wp-fac)*highwsubs_vert
     275
     276             pt_lsa(k,nt)     = fac*lowpt_lsa     + (1.0_wp-fac)*highpt_lsa
     277             q_lsa(k,nt)      = fac*lowq_lsa      + (1.0_wp-fac)*highq_lsa
     278             pt_subs(k,nt)    = fac*lowpt_subs    + (1.0_wp-fac)*highpt_subs
     279             q_subs(k,nt)     = fac*lowq_subs     + (1.0_wp-fac)*highq_subs
    233280
    234281          ENDDO
     
    241288 
    242289       IF ( time_vert(1) > end_time )  THEN
    243           WRITE ( message_string, * ) 'Time dependent large scale profile ',&
    244                              'forcing from&LSF_DATA sets in after end of ' ,&
     290          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
     291                             'forcing from&LSF_DATA sets in after end of ' ,   &
    245292                             'simulation - lsf_vert is set to FALSE'
    246293          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
     
    266313       USE kinds
    267314
    268 
    269315       IMPLICIT NONE
    270316
    271        INTEGER(iwp) ::  t                     !:
     317       INTEGER(iwp) ::  nt                     !:
    272318
    273319       REAL(wp)             :: fac            !:
     
    276322!
    277323!--    Interpolation in time of LSF_DATA at the surface
    278        t = 1
    279        DO WHILE ( time > time_surf(t) )
    280           t = t + 1
     324       nt = 1
     325       DO WHILE ( time > time_surf(nt) )
     326          nt = nt + 1
    281327       ENDDO
    282        IF ( time /= time_surf(t) )  THEN
    283          t = t - 1
    284        ENDIF
    285 
    286        fac = ( time -time_surf(t) ) / ( time_surf(t+1) - time_surf(t) )
     328       IF ( time /= time_surf(nt) )  THEN
     329         nt = nt - 1
     330       ENDIF
     331
     332       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
    287333
    288334       IF ( ibc_pt_b == 0 )  THEN
     
    290336!--       In case of Dirichlet boundary condition shf must not
    291337!--       be set - it is calculated via MOST in prandtl_fluxes
    292           pt_surface = pt_surf(t) + fac * ( pt_surf(t+1) - pt_surf(t) )
     338          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
    293339
    294340       ELSEIF ( ibc_pt_b == 1 )  THEN
     
    296342!--       In case of Neumann boundary condition pt_surface is needed for
    297343!--       calculation of reference density
    298           shf        = shf_surf(t) + fac * ( shf_surf(t+1) - shf_surf(t) )
    299           pt_surface = pt_surf(t) + fac * ( pt_surf(t+1) - pt_surf(t) )
     344          shf        = shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) )
     345          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
    300346
    301347       ENDIF
     
    305351!--       In case of Dirichlet boundary condition qsws must not
    306352!--       be set - it is calculated via MOST in prandtl_fluxes
    307           q_surface = q_surf(t) + fac * ( q_surf(t+1) - q_surf(t) )
     353          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
    308354
    309355       ELSEIF ( ibc_q_b == 1 )  THEN
    310356
    311           qsws = qsws_surf(t) + fac * ( qsws_surf(t+1) - qsws_surf(t) )
    312 
    313        ENDIF
    314 
    315        surface_pressure = p_surf(t) + fac * ( p_surf(t+1) - p_surf(t) )
     357          qsws = qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) )
     358
     359       ENDIF
     360
     361       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
    316362
    317363    END SUBROUTINE ls_forcing_surf
     
    328374       USE kinds
    329375
    330 
    331376       IMPLICIT NONE
    332377
    333        INTEGER(iwp) ::  t                     !:
     378       INTEGER(iwp) ::  nt                     !:
    334379
    335380       REAL(wp)             ::  fac           !:
     
    338383!
    339384!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
    340        t = 1
    341        DO WHILE ( time > time_vert(t) )
    342           t = t + 1
     385       nt = 1
     386       DO WHILE ( time > time_vert(nt) )
     387          nt = nt + 1
    343388       ENDDO
    344        IF ( time /= time_vert(t) )  THEN
    345          t = t - 1
    346        ENDIF
    347 
    348        fac = ( time-time_vert(t) ) / ( time_vert(t+1)-time_vert(t) )
    349 
    350        ug     = ug_vert(:,t) + fac * ( ug_vert(:,t+1) - ug_vert(:,t) )
    351        vg     = vg_vert(:,t) + fac * ( vg_vert(:,t+1) - vg_vert(:,t) )
     389       IF ( time /= time_vert(nt) )  THEN
     390         nt = nt - 1
     391       ENDIF
     392
     393       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
     394
     395       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
     396       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
    352397
    353398       IF ( large_scale_subsidence )  THEN
    354           w_subs = wsubs_vert(:,t) + fac * ( wsubs_vert(:,t+1) - wsubs_vert(:,t) )
     399          w_subs = wsubs_vert(:,nt)                                            &
     400                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
    355401       ENDIF
    356402
     
    358404
    359405
     406!------------------------------------------------------------------------------!
     407! Call for all grid points
     408!------------------------------------------------------------------------------!
     409    SUBROUTINE ls_advec ( time, prog_var )
     410
     411       USE arrays_3d,                                                          &
     412           ONLY:  pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert       
     413       
     414       USE control_parameters,                                                 &
     415           ONLY:  large_scale_subsidence, use_subsidence_tendencies
     416       
     417       USE indices
     418       
     419       USE kinds
     420
     421       IMPLICIT NONE
     422
     423       CHARACTER (LEN=*) ::  prog_var   !:
     424
     425       REAL(wp), INTENT(in)  :: time    !:
     426       REAL(wp) :: fac                  !: 
     427
     428       INTEGER(iwp) ::  i               !:
     429       INTEGER(iwp) ::  j               !:
     430       INTEGER(iwp) ::  k               !:
     431       INTEGER(iwp) ::  nt               !:
     432
     433!
     434!--    Interpolation in time of LSF_DATA
     435       nt = 1
     436       DO WHILE ( time > time_vert(nt) )
     437          nt = nt + 1
     438       ENDDO
     439       IF ( time /= time_vert(nt) )  THEN
     440         nt = nt - 1
     441       ENDIF
     442
     443       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
     444
     445!
     446!--    Add horizontal large scale advection tendencies of pt and q
     447       SELECT CASE ( prog_var )
     448
     449          CASE ( 'pt' )
     450
     451             DO  i = nxl, nxr
     452                DO  j = nys, nyn
     453                   DO  k = nzb_u_inner(j,i)+1, nzt
     454                      tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt)                 &
     455                                    + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
     456                   ENDDO
     457                ENDDO
     458             ENDDO
     459
     460          CASE ( 'q' )
     461
     462             DO  i = nxl, nxr
     463                DO  j = nys, nyn
     464                   DO  k = nzb_u_inner(j,i)+1, nzt
     465                      tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt)                  &
     466                                    + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
     467                   ENDDO
     468                ENDDO
     469             ENDDO
     470
     471       END SELECT
     472
     473!
     474!--    Subsidence of pt and q with prescribed subsidence tendencies
     475       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
     476
     477          SELECT CASE ( prog_var )
     478
     479             CASE ( 'pt' )
     480
     481                DO  i = nxl, nxr
     482                   DO  j = nys, nyn
     483                      DO  k = nzb_u_inner(j,i)+1, nzt
     484                         tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac *     &
     485                                       ( pt_subs(k,nt+1) - pt_subs(k,nt) )
     486                      ENDDO
     487                   ENDDO
     488                ENDDO
     489 
     490             CASE ( 'q' )
     491
     492                DO  i = nxl, nxr
     493                   DO  j = nys, nyn
     494                      DO  k = nzb_u_inner(j,i)+1, nzt
     495                         tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac *      &
     496                                       ( q_subs(k,nt+1) - q_subs(k,nt) )
     497                      ENDDO
     498                   ENDDO
     499                ENDDO
     500
     501          END SELECT
     502
     503       ENDIF
     504
     505    END SUBROUTINE ls_advec
     506
     507
     508!------------------------------------------------------------------------------!
     509! Call for grid point i,j
     510!------------------------------------------------------------------------------!
     511    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
     512
     513       USE arrays_3d,                                                          &
     514           ONLY:  pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert       
     515       
     516       USE control_parameters,                                                 &
     517           ONLY:  large_scale_subsidence, use_subsidence_tendencies
     518       
     519       USE indices
     520       
     521       USE kinds
     522
     523       IMPLICIT NONE
     524
     525       CHARACTER (LEN=*) ::  prog_var   !:
     526
     527       REAL(wp), INTENT(in)  :: time    !:
     528       REAL(wp) :: fac                  !:
     529
     530       INTEGER(iwp) ::  i               !:
     531       INTEGER(iwp) ::  j               !:
     532       INTEGER(iwp) ::  k               !:
     533       INTEGER(iwp) ::  nt               !:
     534
     535!
     536!--    Interpolation in time of LSF_DATA
     537       nt = 1
     538       DO WHILE ( time > time_vert(nt) )
     539          nt = nt + 1
     540       ENDDO
     541       IF ( time /= time_vert(nt) )  THEN
     542         nt = nt - 1
     543       ENDIF
     544
     545       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
     546
     547!
     548!--    Add horizontal large scale advection tendencies of pt and q
     549       SELECT CASE ( prog_var )
     550
     551          CASE ( 'pt' )
     552
     553             DO  k = nzb_u_inner(j,i)+1, nzt
     554                tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt)                       &
     555                              + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) )
     556             ENDDO
     557
     558          CASE ( 'q' )
     559
     560             DO  k = nzb_u_inner(j,i)+1, nzt
     561                tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt)                        &
     562                              + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) )
     563             ENDDO
     564
     565       END SELECT
     566
     567!
     568!--    Subsidence of pt and q with prescribed profiles
     569       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
     570
     571          SELECT CASE ( prog_var )
     572
     573             CASE ( 'pt' )
     574
     575                DO  k = nzb_u_inner(j,i)+1, nzt
     576                   tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac *           &
     577                                 ( pt_subs(k,nt+1) - pt_subs(k,nt) )
     578                ENDDO
     579 
     580             CASE ( 'q' )
     581
     582                DO  k = nzb_u_inner(j,i)+1, nzt
     583                   tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac *            &
     584                                 ( q_subs(k,nt+1) - q_subs(k,nt) )
     585                ENDDO
     586
     587          END SELECT
     588
     589       ENDIF
     590
     591    END SUBROUTINE ls_advec_ij
     592
     593
    360594 END MODULE ls_forcing_mod
  • palm/trunk/SOURCE/modules.f90

    r1362 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Usage of large scale forcing enabled:
     23! increase pr_palm from 90 to 100 to allow for more standard profiles
     24! + ngp_sums_ls, pt_lsa, pt_subs, q_lsa, q_subs, tmp_tnudge, sums_ls_l,
     25! use_subsidence_tendencies
    2326!
    2427! Former revisions:
     
    273276          nc_1d, nr_1d, ptdf_x, ptdf_y, p_surf, pt_surf, pt_1d, pt_init,       &
    274277          qsws_surf, q_1d, q_init, q_surf, qc_1d, qr_1d, rdf, rdf_sc,          &
    275           ref_state, sa_init, shf_surf, timenudge, time_surf, time_vert, ug,   &
    276           u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
     278          ref_state, sa_init, shf_surf, timenudge, time_surf, time_vert,       &
     279          tmp_tnudge, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,                &
     280          v_nzb_p1_for_vfc, w_subs, zu, zw
    277281
    278282    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
     
    281285          flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
    282286          flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
    283           mean_inflow_profiles, nrs, nrsws, nrswst, ptnudge, pt_slope_ref,     &
    284           qnudge, qs, qsws, qswst, qswst_remote, qrs, qrsws, qrswst, rif,      &
    285           saswsb, saswst, shf, tnudge, total_2d_a, total_2d_o, ts, tswst,      &
    286           ug_vert, unudge, us, usws, uswst, vnudge, vg_vert, vsws, vswst,      &
    287           wnudge, wsubs_vert, z0, z0h
     287          mean_inflow_profiles, nrs, nrsws, nrswst, ptnudge, pt_lsa,           &
     288          pt_slope_ref, pt_subs, qnudge, qs, qsws, qswst, qswst_remote, qrs,   &
     289          qrsws, qrswst, q_lsa, q_subs, rif, saswsb, saswst, shf, tnudge,      &
     290          total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us, usws, uswst, &
     291          vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, z0, z0h
    288292
    289293    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
     
    638642                use_prescribed_profile_data = .FALSE., &
    639643                use_single_reference_value = .FALSE., &
     644                use_subsidence_tendencies = .FALSE., &
    640645                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
    641646                use_ug_for_galilei_tr = .TRUE., use_upstream_for_tke = .FALSE.,&
     
    921926    USE kinds
    922927
    923     INTEGER(iwp) ::  i_left, i_right, j_north, j_south, nbgp = 3, ngp_sums, nnx,    &
    924                      nx = 0, nx_a, nx_o, nxl, nxlg, nxlu, nxr, nxrg, nx_on_file,    &
    925                      nny, ny = 0, ny_a, ny_o, nyn, nyng, nys, nysg, nysv,           &
    926                      ny_on_file, nnz, nz = 0, nzb, nzb_diff, nzb_max, nzt, nzt_diff
     928    INTEGER(iwp) ::  i_left, i_right, j_north, j_south, nbgp = 3, ngp_sums,         &
     929                     ngp_sums_ls, nnx, nx = 0, nx_a, nx_o, nxl, nxlg, nxlu, nxr,    &
     930                     nxrg, nx_on_file, nny, ny = 0, ny_a, ny_o, nyn, nyng, nys,     &
     931                     nysg, nysv, ny_on_file, nnz, nz = 0, nzb, nzb_diff, nzb_max,   &
     932                     nzt, nzt_diff
    927933
    928934
     
    13891395
    13901396    CHARACTER (LEN=40) ::  region(0:9)
    1391     INTEGER(iwp) ::  pr_palm = 90, statistic_regions = 0
     1397    INTEGER(iwp) ::  pr_palm = 100, statistic_regions = 0
    13921398    INTEGER(iwp) ::  u_max_ijk(3) = -1, v_max_ijk(3) = -1, w_max_ijk(3) = -1
    13931399    LOGICAL ::  flow_statistics_called = .FALSE.
     
    14031409                                                  sums_wssas_ws_l,                &
    14041410                                                  sums_wsqs_ws_l,                 &
    1405                                                   sums_wsqrs_ws_l
     1411                                                  sums_wsqrs_ws_l,                &
     1412                                                  sums_ls_l
    14061413                                             
    14071414    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  hom_sum, rmask, spectrum_x, &
  • palm/trunk/SOURCE/nudging.f90

    r1356 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Variable t renamed nt, variable currtnudge renamed tmp_tnudge,
     23! summation of nudging tendencies for data output added
     24! +sums_ls_l, tmp_tend
     25! Added new subroutine calc_tnudge, which calculates the current nudging time
     26! scale at each time step
    2327!
    2428! Former revisions:
     
    6569
    6670    PRIVATE
    67     PUBLIC init_nudge, nudge
     71    PUBLIC init_nudge, calc_tnudge, nudge
    6872    SAVE
    6973
     
    7882
    7983       USE arrays_3d,                                                          &
    80            ONLY:  ptnudge, qnudge, timenudge, tnudge, unudge, vnudge, wnudge,  &
    81                   zu
     84           ONLY:  ptnudge, qnudge, timenudge, tmp_tnudge, tnudge, unudge,      &
     85                  vnudge, wnudge, zu
    8286
    8387       USE control_parameters,                                                 &
     
    96100       INTEGER(iwp) ::  ierrn        !:
    97101       INTEGER(iwp) ::  k            !:
    98        INTEGER(iwp) ::  t            !:
     102       INTEGER(iwp) ::  nt            !:
    99103
    100104       CHARACTER(1) ::  hash     !:
     
    122126                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
    123127
     128       ALLOCATE( tmp_tnudge(nzb:nzt) )
     129
    124130       ALLOCATE( timenudge(0:ntnudge) )
    125 
    126131
    127132       ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
    128133       vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
    129 
    130        t = 0
     134!
     135!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
     136       tmp_tnudge = 21600.0_wp
     137
     138       nt = 0
    131139       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
    132140              FORM='FORMATTED', IOSTAT=ierrn )
     
    140148
    141149 rloop:DO
    142           t = t + 1
     150          nt = nt + 1
    143151          hash = "#"
    144152          ierrn = 1 ! not zero
     
    148156          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
    149157         
    150             READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(t)
     158            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
    151159            IF ( ierrn < 0 )  EXIT rloop
    152160
     
    202210             fac = ( highheight - zu(k) ) / ( highheight - lowheight )
    203211
    204              tnudge(k,t)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
    205              unudge(k,t)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
    206              vnudge(k,t)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
    207              wnudge(k,t)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
    208              ptnudge(k,t) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
    209              qnudge(k,t)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
     212             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
     213             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
     214             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
     215             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
     216             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
     217             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
    210218          ENDDO
    211219
     
    225233    END SUBROUTINE init_nudge
    226234
     235
     236    SUBROUTINE calc_tnudge ( time )
     237
     238       USE arrays_3d,                                                          &
     239           ONLY:  timenudge, tmp_tnudge, tnudge
     240
     241       USE control_parameters,                                                 &
     242           ONLY:  dt_3d
     243
     244       USE indices,                                                            &
     245           ONLY:  nzb, nzt
     246
     247       USE kinds
     248
     249       IMPLICIT NONE
     250
     251
     252       REAL(wp) ::  dtm         !:
     253       REAL(wp) ::  dtp         !:
     254       REAL(wp) ::  time        !:
     255
     256       INTEGER(iwp) ::  k   !:
     257       INTEGER(iwp) ::  nt  !:
     258
     259       nt = 1
     260       DO WHILE ( time > timenudge(nt) )
     261         nt = nt+1
     262       ENDDO
     263       IF ( time /= timenudge(1) ) THEN
     264         nt = nt-1
     265       ENDIF
     266
     267       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
     268       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
     269
     270       DO  k = nzb, nzt
     271          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
     272       ENDDO
     273
     274    END SUBROUTINE calc_tnudge
     275
    227276!------------------------------------------------------------------------------!
    228277! Call for all grid points
     
    231280
    232281       USE arrays_3d,                                                          &
    233            ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tnudge, u, unudge,  &
    234                   v, vnudge
    235 
    236        USE buoyancy_mod,                                                       &
    237            ONLY:  calc_mean_profile
     282           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
     283                  u, unudge, v, vnudge
    238284
    239285       USE control_parameters,                                                 &
    240            ONLY:  dt_3d, message_string
     286           ONLY:  dt_3d, intermediate_timestep_count, message_string
    241287
    242288       USE indices,                                                            &
    243289           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
    244290
    245        USE kinds,                                                              &
    246            ONLY:  iwp, wp
     291       USE kinds
    247292
    248293       USE statistics,                                                         &
    249            ONLY:  hom
     294           ONLY:  hom, sums_ls_l, weight_pres
    250295
    251296       IMPLICIT NONE
     
    253298       CHARACTER (LEN=*) ::  prog_var  !:
    254299
    255        REAL(wp) ::  currtnudge  !:
     300       REAL(wp) ::  tmp_tend    !:
    256301       REAL(wp) ::  dtm         !:
    257302       REAL(wp) ::  dtp         !:
     
    261306       INTEGER(iwp) ::  j  !:
    262307       INTEGER(iwp) ::  k  !:
    263        INTEGER(iwp) ::  t  !:
    264 
    265 
    266        t = 1
    267        DO WHILE ( time > timenudge(t) )
    268          t = t+1
     308       INTEGER(iwp) ::  nt  !:
     309
     310
     311       nt = 1
     312       DO WHILE ( time > timenudge(nt) )
     313         nt = nt+1
    269314       ENDDO
    270315       IF ( time /= timenudge(1) ) THEN
    271          t = t-1
     316         nt = nt-1
    272317       ENDIF
    273318
    274        dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
    275        dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )
     319       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
     320       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
    276321
    277322       SELECT CASE ( prog_var )
    278323
    279324          CASE ( 'u' )
    280 
    281              CALL calc_mean_profile( u, 1, 'nudging' )
    282325
    283326             DO  i = nxl, nxr
     
    285328                   DO  k = nzb_u_inner(j,i)+1, nzt
    286329
    287                       currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
    288                                                tnudge(k,t+1) * dtm )
    289 
    290 
    291                       tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
    292                                     ( unudge(k,t) * dtp + &
    293                                        unudge(k,t+1) * dtm ) ) / currtnudge
     330                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
     331                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     332
     333                      tend(k,j,i) = tend(k,j,i) + tmp_tend
     334
     335                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *            &
     336                                       weight_pres(intermediate_timestep_count)
    294337                   ENDDO
    295338                ENDDO
     
    298341          CASE ( 'v' )
    299342
    300              CALL calc_mean_profile( v, 2, 'nudging' )
    301 
    302343             DO  i = nxl, nxr
    303344                DO  j = nys, nyn
    304345                   DO  k = nzb_u_inner(j,i)+1, nzt
    305346
    306                       currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
    307                                                tnudge(k,t+1) * dtm )
    308 
    309 
    310                       tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
    311                                     ( vnudge(k,t) * dtp + &
    312                                        vnudge(k,t+1) * dtm ) ) / currtnudge
     347                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
     348                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     349
     350                      tend(k,j,i) = tend(k,j,i) + tmp_tend
     351
     352                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *            &
     353                                       weight_pres(intermediate_timestep_count)
    313354                   ENDDO
    314355                ENDDO
     
    317358          CASE ( 'pt' )
    318359
    319              CALL calc_mean_profile( pt, 4, 'nudging' )
    320 
    321360             DO  i = nxl, nxr
    322361                DO  j = nys, nyn
    323362                   DO  k = nzb_u_inner(j,i)+1, nzt
    324363
    325                       currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
    326                                                tnudge(k,t+1) * dtm )
    327 
    328 
    329                       tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
    330                                     ( ptnudge(k,t) * dtp + &
    331                                        ptnudge(k,t+1) * dtm ) ) / currtnudge
     364                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
     365                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     366
     367                      tend(k,j,i) = tend(k,j,i) + tmp_tend
     368
     369                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *            &
     370                                       weight_pres(intermediate_timestep_count)
    332371                   ENDDO
    333372                ENDDO
     
    336375          CASE ( 'q' )
    337376
    338              CALL calc_mean_profile( q, 41, 'nudging' )
    339 
    340377             DO  i = nxl, nxr
    341378                DO  j = nys, nyn
    342379                   DO  k = nzb_u_inner(j,i)+1, nzt
    343380
    344                       currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
    345                                                tnudge(k,t+1) * dtm )
    346 
    347 
    348                       tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
    349                                     ( qnudge(k,t) * dtp + &
    350                                        qnudge(k,t+1) * dtm ) ) / currtnudge
     381                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
     382                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     383
     384                      tend(k,j,i) = tend(k,j,i) + tmp_tend
     385
     386                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *            &
     387                                       weight_pres(intermediate_timestep_count)
    351388                   ENDDO
    352389                ENDDO
     
    369406
    370407       USE arrays_3d,                                                          &
    371            ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tnudge, u, unudge,  &
    372                   v, vnudge
    373 
    374        USE buoyancy_mod,                                                       &
    375            ONLY:  calc_mean_profile
     408           ONLY:  pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
     409                  u, unudge, v, vnudge
    376410
    377411       USE control_parameters,                                                 &
    378            ONLY:  dt_3d, message_string
     412           ONLY:  dt_3d, intermediate_timestep_count, message_string
    379413
    380414       USE indices,                                                            &
    381415           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
    382416
    383        USE kinds,                                                              &
    384            ONLY:  iwp, wp
     417       USE kinds
    385418
    386419       USE statistics,                                                         &
    387            ONLY:  hom
     420           ONLY:  hom, sums_ls_l, weight_pres
    388421
    389422       IMPLICIT NONE
     
    392425       CHARACTER (LEN=*) ::  prog_var  !:
    393426
    394        REAL(wp) ::  currtnudge  !:
     427       REAL(wp) ::  tmp_tend    !:
    395428       REAL(wp) ::  dtm         !:
    396429       REAL(wp) ::  dtp         !:
     
    400433       INTEGER(iwp) ::  j  !:
    401434       INTEGER(iwp) ::  k  !:
    402        INTEGER(iwp) ::  t  !:
    403 
    404 
    405        t = 1
    406        DO WHILE ( time > timenudge(t) )
    407          t = t+1
     435       INTEGER(iwp) ::  nt  !:
     436
     437
     438       nt = 1
     439       DO WHILE ( time > timenudge(nt) )
     440         nt = nt+1
    408441       ENDDO
    409442       IF ( time /= timenudge(1) )  THEN
    410          t = t-1
     443         nt = nt-1
    411444       ENDIF
    412445
    413        dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
    414        dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )
     446       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
     447       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
    415448
    416449       SELECT CASE ( prog_var )
     
    418451          CASE ( 'u' )
    419452
    420              CALL calc_mean_profile( u, 1, 'nudging' )
    421 
    422453             DO  k = nzb_u_inner(j,i)+1, nzt
    423454
    424                 currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
    425 
    426                 tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
    427                                ( unudge(k,t) * dtp + &
    428                                  unudge(k,t+1) * dtm ) ) / currtnudge
     455                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
     456                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     457
     458                tend(k,j,i) = tend(k,j,i) + tmp_tend
     459
     460                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
     461                                 * weight_pres(intermediate_timestep_count)
    429462             ENDDO
    430463
    431464          CASE ( 'v' )
    432465
    433              CALL calc_mean_profile( v, 2, 'nudging' )
    434 
    435466             DO  k = nzb_u_inner(j,i)+1, nzt
    436467
    437                 currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
    438 
    439                 tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
    440                                ( vnudge(k,t) * dtp + &
    441                                  vnudge(k,t+1) * dtm ) ) / currtnudge
     468                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
     469                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     470
     471                tend(k,j,i) = tend(k,j,i) + tmp_tend
     472
     473                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
     474                                 * weight_pres(intermediate_timestep_count)
    442475             ENDDO
    443476
     
    445478          CASE ( 'pt' )
    446479
    447              CALL calc_mean_profile( pt, 4, 'nudging' )
    448 
    449480             DO  k = nzb_u_inner(j,i)+1, nzt
    450481
    451                 currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
    452 
    453                 tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
    454                                ( ptnudge(k,t) * dtp + &
    455                                  ptnudge(k,t+1) * dtm ) ) / currtnudge
     482                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
     483                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     484
     485                tend(k,j,i) = tend(k,j,i) + tmp_tend
     486
     487                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
     488                                 * weight_pres(intermediate_timestep_count)
    456489             ENDDO
    457490
     
    459492          CASE ( 'q' )
    460493
    461              CALL calc_mean_profile( q, 41, 'nudging' )
    462 
    463494             DO  k = nzb_u_inner(j,i)+1, nzt
    464495
    465                 currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )
    466 
    467                 tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
    468                                ( qnudge(k,t) * dtp + &
    469                                  qnudge(k,t+1) * dtm ) ) / currtnudge
     496                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
     497                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
     498
     499                tend(k,j,i) = tend(k,j,i) + tmp_tend
     500
     501                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
     502                                 * weight_pres(intermediate_timestep_count)
    470503             ENDDO
    471504
  • palm/trunk/SOURCE/parin.f90

    r1362 r1365  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Usage of large scale forcing enabled:
     23! +use_subsidence_tendencies
    2324!
    2425! Former revisions:
     
    211212               transpose_compute_overlap, turbulence, turbulent_inflow,        &
    212213               ug_surface, ug_vertical_gradient, ug_vertical_gradient_level,   &
    213                use_surface_fluxes, use_cmax, use_top_fluxes,                   &
    214                use_ug_for_galilei_tr, use_upstream_for_tke, uv_heights,        &
    215                u_bulk, u_profile, vg_surface, vg_vertical_gradient,            &
     214               use_subsidence_tendencies, use_surface_fluxes, use_cmax,        &
     215               use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke,    &
     216               uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient,&
    216217               vg_vertical_gradient_level, v_bulk, v_profile,                  &
    217218               wall_adjustment, wall_heatflux, wall_humidityflux,              &
     
    295296             top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
    296297             transpose_compute_overlap, turbulence, turbulent_inflow, &
    297              ug_surface, ug_vertical_gradient, &
     298             use_subsidence_tendencies, ug_surface, ug_vertical_gradient, &
    298299             ug_vertical_gradient_level, use_surface_fluxes, use_cmax, &
    299300             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
  • palm/trunk/SOURCE/prognostic_equations.f90

    r1362 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Calls of ls_advec for large scale advection added,
     23! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
     24! new argument ls_index added to the calls of subsidence
     25! +ls_index
    2326!
    2427! Former revisions:
     
    143146               dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
    144147               icloud_scheme, inflow_l, intermediate_timestep_count,           &
    145                intermediate_timestep_count_max, large_scale_subsidence,        &
    146                neutral, nudging, ocean, outflow_l, outflow_s, passive_scalar,  &
    147                plant_canopy, precipitation, prho_reference, prho_reference,    &
    148                prho_reference, pt_reference, pt_reference, pt_reference,       &
    149                radiation, scalar_advec, scalar_advec, simulated_time,          &
    150                sloping_surface, timestep_scheme, tsc, use_upstream_for_tke,    &
     148               intermediate_timestep_count_max, large_scale_forcing,           &
     149               large_scale_subsidence, neutral, nudging, ocean, outflow_l,     &
     150               outflow_s, passive_scalar, plant_canopy, precipitation,         &
     151               prho_reference, prho_reference, prho_reference, pt_reference,   &
     152               pt_reference, pt_reference, radiation, scalar_advec,            &
     153               scalar_advec, simulated_time, sloping_surface, timestep_scheme, &
     154               tsc, use_subsidence_tendencies, use_upstream_for_tke,           &
    151155               use_upstream_for_tke, use_upstream_for_tke, wall_heatflux,      &
    152156               wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
     
    222226
    223227    USE kinds
     228
     229    USE ls_forcing_mod,                                                        &
     230        ONLY:  ls_advec
    224231
    225232    USE microphysics_mod,                                                      &
     
    307314          i_omp_start = i
    308315       ENDIF
    309  
     316
    310317       DO  j = nys, nyn
    311318!
     
    537544
    538545!
     546!--          Large scale advection
     547             IF ( large_scale_forcing )  THEN
     548                CALL ls_advec( i, j, simulated_time, 'pt' )
     549             ENDIF     
     550
     551!
    539552!--          If required, compute effect of large-scale subsidence/ascent
    540              IF ( large_scale_subsidence )  THEN
    541                 CALL subsidence( i, j, tend, pt, pt_init )
     553             IF ( large_scale_subsidence  .AND.                                &
     554                  .NOT. use_subsidence_tendencies )  THEN
     555                CALL subsidence( i, j, tend, pt, pt_init, 2 )
    542556             ENDIF
    543557
     
    662676
    663677!
     678!--          Large scale advection
     679             IF ( large_scale_forcing )  THEN
     680                CALL ls_advec( i, j, simulated_time, 'q' )
     681             ENDIF
     682
     683!
    664684!--          If required compute influence of large-scale subsidence/ascent
    665              IF ( large_scale_subsidence )  THEN
    666                 CALL subsidence( i, j, tend, q, q_init )
     685             IF ( large_scale_subsidence  .AND.                                &
     686                  .NOT. use_subsidence_tendencies )  THEN
     687                CALL subsidence( i, j, tend, q, q_init, 3 )
    667688             ENDIF
    668689
     
    11771198
    11781199!
     1200!--    Large scale advection
     1201       IF ( large_scale_forcing )  THEN
     1202          CALL ls_advec( simulated_time, 'pt' )
     1203       ENDIF
     1204
     1205!
    11791206!--    If required compute influence of large-scale subsidence/ascent
    1180        IF ( large_scale_subsidence )  THEN
    1181           CALL subsidence( tend, pt, pt_init )
     1207       IF ( large_scale_subsidence  .AND.                                      &
     1208            .NOT. use_subsidence_tendencies )  THEN
     1209          CALL subsidence( tend, pt, pt_init, 2 )
    11821210       ENDIF
    11831211
     
    13651393!--    Sink or source of scalar concentration due to canopy elements
    13661394       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
    1367        
     1395
     1396!
     1397!--    Large scale advection
     1398       IF ( large_scale_forcing )  THEN
     1399          CALL ls_advec( simulated_time, 'q' )
     1400       ENDIF
     1401
    13681402!
    13691403!--    If required compute influence of large-scale subsidence/ascent
    1370        IF ( large_scale_subsidence )  THEN
    1371          CALL subsidence( tend, q, q_init )
     1404       IF ( large_scale_subsidence  .AND.                                      &
     1405            .NOT. use_subsidence_tendencies )  THEN
     1406         CALL subsidence( tend, q, q_init, 3 )
    13721407       ENDIF
    13731408
     
    19782013
    19792014!
     2015!--    Large scale advection
     2016       IF ( large_scale_forcing )  THEN
     2017          CALL ls_advec( simulated_time, 'pt' )
     2018       ENDIF
     2019
     2020!
    19802021!--    If required compute influence of large-scale subsidence/ascent
    1981        IF ( large_scale_subsidence )  THEN
    1982           CALL subsidence( tend, pt, pt_init )
     2022       IF ( large_scale_subsidence  .AND.                                      &
     2023            .NOT. use_subsidence_tendencies )  THEN
     2024          CALL subsidence( tend, pt, pt_init, 2 )
    19832025       ENDIF
    19842026
     
    21422184
    21432185!
     2186!--    Large scale advection
     2187       IF ( large_scale_forcing )  THEN
     2188          CALL ls_advec( simulated_time, 'q' )
     2189       ENDIF
     2190
     2191!
    21442192!--    If required compute influence of large-scale subsidence/ascent
    2145        IF ( large_scale_subsidence )  THEN
    2146          CALL subsidence( tend, q, q_init )
     2193       IF ( large_scale_subsidence  .AND.                                      &
     2194            .NOT. use_subsidence_tendencies )  THEN
     2195         CALL subsidence( tend, q, q_init, 3 )
    21472196       ENDIF
    21482197
  • palm/trunk/SOURCE/subsidence.f90

    r1354 r1365  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Summation of subsidence tendencies for data output added
     23! +ls_index, sums_ls_l, tmp_tend
    2324!
    2425! Former revisions:
     
    9091       REAL(wp)     ::  ws_surface !:
    9192
    92        IF ( .NOT. ALLOCATED( w_subs )) THEN
     93       IF ( .NOT. ALLOCATED( w_subs ))  THEN
    9394          ALLOCATE( w_subs(nzb:nzt+1) )
    9495          w_subs = 0.0_wp
    9596       ENDIF
    9697
    97       IF ( ocean )  THEN
     98       IF ( ocean )  THEN
    9899          message_string = 'Applying large scale vertical motion is not ' // &
    99100                           'allowed for ocean runs'
     
    111112      subs_vertical_gradient_level_i(1) = 0
    112113      DO  k = 1, nzt+1
    113          IF ( i < 11 ) THEN
     114         IF ( i < 11 )  THEN
    114115            IF ( subs_vertical_gradient_level(i) < zu(k)  .AND. &
    115116                 subs_vertical_gradient_level(i) >= 0.0_wp )  THEN
     
    140141
    141142
    142     SUBROUTINE subsidence( tendency, var, var_init )
     143    SUBROUTINE subsidence( tendency, var, var_init, ls_index )
    143144
    144145       USE arrays_3d,                                                          &
     
    146147
    147148       USE control_parameters,                                                 &
    148            ONLY:  dt_3d
     149           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing
    149150
    150151       USE indices,                                                            &
     
    153154
    154155       USE kinds
     156
     157       USE statistics,                                                         &
     158           ONLY:  sums_ls_l, weight_pres
    155159
    156160       IMPLICIT NONE
     
    159163       INTEGER(iwp) ::  j !:
    160164       INTEGER(iwp) ::  k !:
    161 
     165       INTEGER(iwp) ::  ls_index !:
     166
     167       REAL(wp)     ::  tmp_tend !:
    162168       REAL(wp)     ::  tmp_grad !:
    163169   
     
    174180          DO  j = nys, nyn
    175181             DO  k = nzb_s_inner(j,i)+1, nzt
    176                 IF ( w_subs(k) < 0.0_wp ) THEN    ! large-scale subsidence
    177                    tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
    178                                      ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
    179                 ELSE                           ! large-scale ascent
    180                    tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
    181                                      ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
     182                IF ( w_subs(k) < 0.0_wp )  THEN    ! large-scale subsidence
     183                   tmp_tend = - w_subs(k) *                                    &
     184                              ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
     185                ELSE                               ! large-scale ascent
     186                   tmp_tend = - w_subs(k) *                                    &
     187                              ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
     188                ENDIF
     189
     190                tendency(k,j,i) = tendency(k,j,i) + tmp_tend
     191
     192                IF ( large_scale_forcing )  THEN
     193                   sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend    &
     194                                      * weight_pres(intermediate_timestep_count)
    182195                ENDIF
    183196             ENDDO
     
    190203
    191204       DO  k = nzb, nzt
    192           IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
     205          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
    193206             var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
    194207                               ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
     
    198211!--   At the upper boundary, the initial profile is shifted with aid of
    199212!--   the gradient tmp_grad. (This is ok if the gradients are linear.)
    200       IF ( w_subs(nzt) < 0.0_wp ) THEN
     213      IF ( w_subs(nzt) < 0.0_wp )  THEN
    201214         tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
    202215         var_mod(nzt+1) = var_init(nzt+1) -  &
     
    206219
    207220      DO  k = nzt+1, nzb+1, -1
    208          IF ( w_subs(k) >= 0.0_wp ) THEN  ! large-scale ascent
     221         IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
    209222            var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
    210223                               ( var_init(k) - var_init(k-1) ) * ddzu(k)
     
    214227!--   At the lower boundary shifting is not necessary because the
    215228!--   subsidence velocity w_subs(nzb) vanishes.
    216 
    217 
    218       IF ( w_subs(nzb+1) >= 0.0_wp ) THEN
     229      IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
    219230         var_mod(nzb) = var_init(nzb)
    220231      ENDIF
     
    225236 END SUBROUTINE subsidence
    226237
    227  SUBROUTINE subsidence_ij( i, j, tendency, var, var_init )
     238 SUBROUTINE subsidence_ij( i, j, tendency, var, var_init, ls_index )
    228239
    229240       USE arrays_3d,                                                          &
     
    231242
    232243       USE control_parameters,                                                 &
    233            ONLY:  dt_3d
     244           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing
    234245
    235246       USE indices,                                                            &
     
    237248
    238249       USE kinds
     250
     251       USE statistics,                                                         &
     252           ONLY:  sums_ls_l, weight_pres
    239253
    240254       IMPLICIT NONE
     
    243257       INTEGER(iwp) ::  j !:
    244258       INTEGER(iwp) ::  k !:
    245 
     259       INTEGER(iwp) ::  ls_index !:
     260
     261       REAL(wp)     ::  tmp_tend !:
    246262       REAL(wp)     ::  tmp_grad !:
    247263   
     
    256272!--    Influence of w_subsidence on the current tendency term
    257273       DO  k = nzb_s_inner(j,i)+1, nzt
    258           IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
    259              tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
    260                                ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
    261           ELSE                             ! large-scale ascent
    262              tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
    263                                ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
     274          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
     275             tmp_tend = - w_subs(k) * ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
     276          ELSE                                 ! large-scale ascent
     277             tmp_tend = - w_subs(k) * ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
     278          ENDIF
     279
     280          tendency(k,j,i) = tendency(k,j,i) + tmp_tend
     281
     282          IF ( large_scale_forcing )  THEN
     283             sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend          &
     284                                     * weight_pres(intermediate_timestep_count)
    264285          ENDIF
    265286       ENDDO
     
    269290!--    Shifting of the initial profile is especially necessary with Rayleigh
    270291!--    damping switched on
    271        IF ( i == nxl .AND. j == nys ) THEN ! shifting only once per PE
     292       IF ( i == nxl .AND. j == nys )  THEN ! shifting only once per PE
    272293
    273294          DO  k = nzb, nzt
    274              IF ( w_subs(k) < 0.0_wp ) THEN      ! large-scale subsidence
     295             IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
    275296                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
    276297                                  ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
     
    280301!--       At the upper boundary, the initial profile is shifted with aid of
    281302!--       the gradient tmp_grad. (This is ok if the gradients are linear.)
    282           IF ( w_subs(nzt) < 0.0_wp ) THEN
     303          IF ( w_subs(nzt) < 0.0_wp )  THEN
    283304             tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
    284305             var_mod(nzt+1) = var_init(nzt+1) -  &
     
    288309
    289310          DO  k = nzt+1, nzb+1, -1
    290              IF ( w_subs(k) >= 0.0_wp ) THEN  ! large-scale ascent
     311             IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
    291312                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
    292313                                   ( var_init(k) - var_init(k-1) ) * ddzu(k)
     
    296317!--       At the lower boundary shifting is not necessary because the
    297318!--       subsidence velocity w_subs(nzb) vanishes.
    298 
    299 
    300           IF ( w_subs(nzb+1) >= 0.0_wp ) THEN
     319          IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
    301320             var_mod(nzb) = var_init(nzb)
    302321          ENDIF
  • palm/trunk/SOURCE/time_integration.f90

    r1343 r1365  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 !
     22! Reset sums_ls_l to zero at each timestep
     23! +sums_ls_l
     24! Calculation of reference state (previously in subroutine calc_mean_profile)
     25
    2426! Former revisions:
    2527! -----------------
     
    126128
    127129    USE arrays_3d,                                                             &
    128         ONLY:  diss, e_p, nr_p, prho, pt, pt_p, ql, ql_c, ql_v, ql_vp, qr_p,   &
    129                q_p, rho, sa_p, tend, u, u_p, v, vpt, v_p, w_p
    130 
    131     USE buoyancy_mod,                                                          &
     130        ONLY:  diss, e_p, nr_p, prho, pt, pt_p, q, ql, ql_c, ql_v, ql_vp, qr_p,&
     131               q_p, ref_state, rho, sa_p, tend, u, u_p, v, vpt, v_p, w_p
     132
     133    USE calc_mean_profile_mod,                                                 &
    132134        ONLY:  calc_mean_profile
    133135
     
    149151               intermediate_timestep_count_max, large_scale_forcing,           &
    150152               loop_optimization, lsf_surf, lsf_vert, masks, mid,              &
    151                netcdf_data_format, neutral, nr_timesteps_this_run, ocean,      &
    152                on_device, passive_scalar, prandtl_layer, precipitation,        &
     153               netcdf_data_format, neutral, nr_timesteps_this_run, nudging,    &
     154               ocean, on_device, passive_scalar, prandtl_layer, precipitation, &
    153155               prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
    154156               run_coupled, simulated_time, simulated_time_chr,                &
     
    181183        ONLY:  ls_forcing_surf, ls_forcing_vert
    182184
     185    USE nudge_mod,                                                             &
     186        ONLY:  calc_tnudge
     187
    183188    USE particle_attributes,                                                   &
    184189        ONLY:  particle_advection, particle_advection_start, wang_kernel
     
    194199
    195200    USE statistics,                                                            &
    196         ONLY:  flow_statistics_called, hom, pr_palm
     201        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
    197202
    198203    USE user_actions_mod,                                                      &
     
    253258!--    Determine ug, vg and w_subs in dependence on data from external file
    254259!--    LSF_DATA
    255        IF ( large_scale_forcing .AND. lsf_vert ) THEN
     260       IF ( large_scale_forcing .AND. lsf_vert )  THEN
    256261           CALL ls_forcing_vert ( simulated_time )
     262           sums_ls_l = 0.0_wp
    257263       ENDIF
    258264
     
    283289!--          buoyancy terms (WARNING: only the respective last call of
    284290!--          calc_mean_profile defines the reference state!)
    285              IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4, 'time_int' )
    286              IF ( ocean         )  CALL calc_mean_profile( rho, 64, 'time_int' )
    287              IF ( humidity      )  CALL calc_mean_profile( vpt, 44, 'time_int' )
     291             IF ( .NOT. neutral )  THEN
     292                CALL calc_mean_profile( pt, 4 )
     293                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
     294             ENDIF
     295             IF ( ocean )  THEN
     296                CALL calc_mean_profile( rho, 64 )
     297                ref_state(:)  = hom(:,1,64,0)
     298             ENDIF
     299             IF ( humidity )  THEN
     300                CALL calc_mean_profile( vpt, 44 )
     301                ref_state(:)  = hom(:,1,44,0)
     302             ENDIF
     303
    288304          ENDIF
    289305
     
    291307          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
    292308               intermediate_timestep_count == 1 )  CALL ws_statistics
     309!
     310!--       In case of nudging calculate current nudging time scale and horizontal
     311!--       means of u,v,pt and q
     312          IF ( nudging )  THEN
     313             CALL calc_tnudge( simulated_time )
     314             CALL calc_mean_profile( u, 1 )
     315             CALL calc_mean_profile( v, 2 )
     316             CALL calc_mean_profile( pt, 4 )
     317             CALL calc_mean_profile( q, 41 )
     318          ENDIF
    293319
    294320!
Note: See TracChangeset for help on using the changeset viewer.