Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 1361 for palm/trunk/SOURCE


Ignore:
Timestamp:
Apr 16, 2014 3:17:48 PM (11 years ago)
Author:
hoffmann
Message:

improved version of two-moment cloud physics

Location:
palm/trunk/SOURCE
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r1360 r1361  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# cpulog added to microphysics
    2323#
    2424# Former revisions:
    … …  
    341341ls_forcing.o: modules.o cpulog.o mod_kinds.o
    342342message.o: modules.o mod_kinds.o
    343 microphysics.o: modules.o mod_kinds.o
     343microphysics.o: modules.o cpulog.o mod_kinds.o
    344344modules.o: modules.f90 mod_kinds.o
    345345mod_kinds.o: mod_kinds.f90

  • palm/trunk/SOURCE/advec_s_bc.f90

    r1354 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! nr and qr added
    2323!
    2424! Former revisions:
    … …  
    926926       ENDIF
    927927
     928    ELSEIF ( sk_char == 'qr' )  THEN
     929
     930!
     931!--    Rain water content boundary condition at the bottom boundary:
     932!--    Dirichlet (fixed surface rain water content).
     933       DO  i = nxl, nxr
     934          DO  j = nys, nyn
     935             sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     936             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     937             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     938          ENDDO
     939       ENDDO
     940
     941!
     942!--    Rain water content boundary condition at the top boundary: Dirichlet
     943       DO  i = nxl, nxr
     944          DO  j = nys, nyn
     945             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     946             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     947          ENDDO
     948       ENDDO
     949
     950    ELSEIF ( sk_char == 'nr' )  THEN
     951
     952!
     953!--    Rain drop concentration boundary condition at the bottom boundary:
     954!--    Dirichlet (fixed surface rain drop concentration).
     955       DO  i = nxl, nxr
     956          DO  j = nys, nyn
     957             sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     958             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     959             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     960          ENDDO
     961       ENDDO
     962
     963!
     964!--    Rain drop concentration boundary condition at the top boundary: Dirichlet
     965       DO  i = nxl, nxr
     966          DO  j = nys, nyn
     967             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     968             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     969          ENDDO
     970       ENDDO
     971
    928972    ELSEIF ( sk_char == 'e' )  THEN
    929973

  • palm/trunk/SOURCE/advec_ws.f90

    r1354 r1361  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! accelerator and vector version for qr and nr added
    2323!
    2424! Former revisions:
    … …  
    22332233       USE statistics,                                                        &
    22342234           ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,           &
    2235                   weight_substep
     2235                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
    22362236
    22372237       IMPLICIT NONE
    … …  
    26562656                        *   weight_substep(intermediate_timestep_count)
    26572657                    ENDDO
     2658                 CASE ( 'qr' )
     2659                    DO  k = nzb, nzt
     2660                       sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
     2661                        + ( flux_t(k) + diss_t(k) )                          &
     2662                        *   weight_substep(intermediate_timestep_count)
     2663                    ENDDO
     2664                 CASE ( 'nr' )
     2665                    DO  k = nzb, nzt
     2666                       sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
     2667                        + ( flux_t(k) + diss_t(k) )                          &
     2668                        *   weight_substep(intermediate_timestep_count)
     2669                    ENDDO
    26582670
    26592671              END SELECT
    … …  
    26902702!        USE statistics,                                                       &
    26912703!            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
    2692 !                   weight_substep
     2704!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
    26932705
    26942706       IMPLICIT NONE
    … …  
    29752987!                   CASE ( 'q' )
    29762988!                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
    2977 !                      + ( flux_t + diss_t )                                &
     2989!                      + ( flux_t + diss_t )                                 &
     2990!                      *   weight_substep(intermediate_timestep_count)
     2991!                   CASE ( 'qr' )
     2992!                      sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
     2993!                      + ( flux_t + diss_t )                                 &
     2994!                      *   weight_substep(intermediate_timestep_count)
     2995!                   CASE ( 'nr' )
     2996!                      sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
     2997!                      + ( flux_t + diss_t )                                 &
    29782998!                      *   weight_substep(intermediate_timestep_count)
    29792999!

  • palm/trunk/SOURCE/boundary_conds.f90

    r1354 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bottom and top boundary conditions of rain water content (qr) and
     23! rain drop concentration (nr) changed to Dirichlet
    2324!
    2425! Former revisions:
    … …  
    271272       q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
    272273
    273        IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                    &
    274             precipitation )  THEN
     274       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
    275275!             
    276 !--       Surface conditions rain water (Neumann)
     276!--       Surface conditions rain water (Dirichlet)
    277277          DO  i = nxlg, nxrg
    278278             DO  j = nysg, nyng
    279                 qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
    280                 nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
    281              ENDDO
    282           ENDDO
    283 !
    284 !--       Top boundary condition for rain water (Neumann)
    285           qr_p(nzt+1,:,:) = qr_p(nzt,:,:)
    286           nr_p(nzt+1,:,:) = nr_p(nzt,:,:)
     279                qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
     280                nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
     281             ENDDO
     282          ENDDO
     283!
     284!--       Top boundary condition for rain water (Dirichlet)
     285          qr_p(nzt+1,:,:) = 0.0_wp
     286          nr_p(nzt+1,:,:) = 0.0_wp
    287287           
    288288       ENDIF

  • palm/trunk/SOURCE/check_parameters.f90

    r1360 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
    23 !
     22! PA0363 removed
     23! PA0362 changed
     24!
    2425! Former revisions:
    2526! -----------------
    … …  
    862863    ENDIF
    863864
    864     IF ( loop_optimization /= 'cache' .AND.  cloud_physics  .AND.            &
    865          icloud_scheme == 0 ) THEN
     865    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
     866                 loop_optimization == 'vector' )                               &
     867         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
    866868       message_string = 'cloud_scheme = seifert_beheng requires ' // &
    867                         'loop_optimization = cache'
     869                        'loop_optimization = "cache" or "vector"'
    868870       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
    869871    ENDIF
    870 
    871 !    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    872 !         .NOT. precipitation  .AND.  .NOT. drizzle ) THEN
    873 !       message_string = 'cloud_scheme = seifert_beheng requires ' // &
    874 !                        'precipitation = .TRUE. or drizzle = .TRUE.'
    875 !       CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )
    876 !    ENDIF
    877872
    878873!

  • palm/trunk/SOURCE/init_3d_model.f90

    r1360 r1361  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! tend_* removed
     23! Bugfix: w_subs is not allocated anymore if it is already allocated
    2324!
    2425! Former revisions:
    … …  
    366367                ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    367368#endif
    368 !
    369 !--             3D-tendency arrays
    370                 ALLOCATE( tend_pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    371                           tend_q(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    372369
    373370                IF ( precipitation )  THEN
    … …  
    375372!--                1D-arrays
    376373                   ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
    377 !
    378 !
    379 !--                3D-tendency arrays
    380                    ALLOCATE( tend_nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    381                              tend_qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     374
    382375!
    383376!--                2D-rain water content and rain drop concentration arrays
    … …  
    475468!
    476469!-- 1D-array for large scale subsidence velocity
    477     ALLOCATE ( w_subs(nzb:nzt+1) )
    478     w_subs = 0.0_wp
    479 
     470    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
     471       ALLOCATE ( w_subs(nzb:nzt+1) )
     472       w_subs = 0.0_wp
     473    ENDIF
    480474
    481475!

  • palm/trunk/SOURCE/init_cloud_physics.f90

    r1354 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! sed_qc_const is now calculated here (2-moment microphysics)
    2323!
    2424! Former revisions:
    … …  
    7373    USE cloud_parameters,                                                      &
    7474        ONLY:  bfactor, cp, c_sedimentation, dpirho_l, dt_precipitation,       &
    75                hyrho, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute,              &
     75               hyrho, k_st, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute,  &
    7676               molecular_weight_of_solute, molecular_weight_of_water, pirho_l, &
    77                pt_d_t, rho_l, r_d, r_v, schmidt, schmidt_p_1d3, t_d_pt,        &
    78                vanthoff, w_precipitation
     77               pt_d_t, rho_l, r_d, r_v, sed_qc_const, schmidt, schmidt_p_1d3,  &
     78               sigma_gc, t_d_pt, vanthoff, w_precipitation
    7979       
    8080    USE constants,                                                             &
    … …  
    108108
    109109    pirho_l  = pi * rho_l / 6.0_wp
    110     dpirho_l = 1.0_wp / pirho_l
     110    dpirho_l = 1.0_wp / pirho_l
     111!
     112!-- constant fot the sedimentation of cloud water (2-moment cloud physics)
     113    sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                   &
     114                          )**( 2.0_wp / 3.0_wp ) *                             &
     115                   EXP( 5.0_wp * LOG( sigma_gc )**2 )
    111116!
    112117!-- Calculate timestep according to precipitation

  • palm/trunk/SOURCE/microphysics.f90

    r1354 r1361  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix in sedimentation_rain: Index corrected.
     23! Vectorized version of adjust_cloud added.
     24! Little reformatting of the code.
    2325!
    2426! Former revisions:
    … …  
    126128    SUBROUTINE microphysics_control
    127129
    128        USE arrays_3d
    129        USE cloud_parameters
    130        USE control_parameters
    131        USE grid_variables
    132        USE indices
     130       USE arrays_3d,                                                          &
     131           ONLY:  hyp, nr, pt, pt_init, q, qc, qr, zu
     132
     133       USE cloud_parameters,                                                   &
     134           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
     135
     136       USE control_parameters,                                                 &
     137           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
     138                  g, intermediate_timestep_count,                              &
     139                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
     140                  rho_surface,surface_pressure
     141
     142       USE indices,                                                            &
     143           ONLY:  nzb, nzt
     144
    133145       USE kinds
    134        USE statistics
     146
     147       USE statistics,                                                         &
     148           ONLY:  weight_pres
    135149
    136150       IMPLICIT NONE
    … …  
    139153       INTEGER(iwp) ::  j                 !:
    140154       INTEGER(iwp) ::  k                 !:
    141  
    142        DO  i = nxl, nxr
    143           DO  j = nys, nyn
    144              DO  k = nzb_s_inner(j,i)+1, nzt
    145 
    146              ENDDO
    147           ENDDO
    148        ENDDO
    149 
    150     END SUBROUTINE microphysics_control
    151 
    152     SUBROUTINE adjust_cloud
    153 
    154        USE arrays_3d
    155        USE cloud_parameters
    156        USE indices
    157        USE kinds
    158 
    159        IMPLICIT NONE
    160 
    161        INTEGER(iwp) ::  i                 !:
    162        INTEGER(iwp) ::  j                 !:
    163        INTEGER(iwp) ::  k                 !:
    164 
    165  
    166        DO  i = nxl, nxr
    167           DO  j = nys, nyn
    168              DO  k = nzb_s_inner(j,i)+1, nzt
    169 
    170              ENDDO
    171           ENDDO
    172        ENDDO
    173 
    174     END SUBROUTINE adjust_cloud
    175 
    176 
    177     SUBROUTINE autoconversion
    178 
    179        USE arrays_3d
    180        USE cloud_parameters
    181        USE control_parameters
    182        USE grid_variables
    183        USE indices
    184        USE kinds
    185 
    186        IMPLICIT NONE
    187 
    188        INTEGER(iwp) ::  i                 !:
    189        INTEGER(iwp) ::  j                 !:
    190        INTEGER(iwp) ::  k                 !:
    191 
    192        DO  i = nxl, nxr
    193           DO  j = nys, nyn
    194              DO  k = nzb_s_inner(j,i)+1, nzt
    195 
    196              ENDDO
    197           ENDDO
    198        ENDDO
    199 
    200     END SUBROUTINE autoconversion
    201 
    202 
    203     SUBROUTINE accretion
    204 
    205        USE arrays_3d
    206        USE cloud_parameters
    207        USE control_parameters
    208        USE indices
    209        USE kinds
    210 
    211        IMPLICIT NONE
    212 
    213        INTEGER(iwp) ::  i                 !:
    214        INTEGER(iwp) ::  j                 !:
    215        INTEGER(iwp) ::  k                 !:
    216 
    217        DO  i = nxl, nxr
    218           DO  j = nys, nyn
    219              DO  k = nzb_s_inner(j,i)+1, nzt
    220 
    221              ENDDO
    222           ENDDO
    223        ENDDO
    224 
    225     END SUBROUTINE accretion
    226 
    227 
    228     SUBROUTINE selfcollection_breakup
    229 
    230        USE arrays_3d
    231        USE cloud_parameters
    232        USE control_parameters
    233        USE indices
    234        USE kinds
    235 
    236        IMPLICIT NONE
    237 
    238        INTEGER(iwp) ::  i                 !:
    239        INTEGER(iwp) ::  j                 !:
    240        INTEGER(iwp) ::  k                 !:
    241 
    242  
    243        DO  i = nxl, nxr
    244           DO  j = nys, nyn
    245              DO  k = nzb_s_inner(j,i)+1, nzt
    246 
    247              ENDDO
    248           ENDDO
    249        ENDDO
    250 
    251     END SUBROUTINE selfcollection_breakup
    252 
    253 
    254     SUBROUTINE evaporation_rain
    255 
    256        USE arrays_3d
    257        USE cloud_parameters
    258        USE constants
    259        USE control_parameters
    260        USE indices
    261        USE kinds
    262 
    263        IMPLICIT NONE
    264 
    265        INTEGER(iwp) ::  i                 !:
    266        INTEGER(iwp) ::  j                 !:
    267        INTEGER(iwp) ::  k                 !:
    268  
    269        DO  i = nxl, nxr
    270           DO  j = nys, nyn
    271              DO  k = nzb_s_inner(j,i)+1, nzt
    272 
    273              ENDDO
    274           ENDDO
    275        ENDDO
    276 
    277     END SUBROUTINE evaporation_rain
    278 
    279 
    280     SUBROUTINE sedimentation_cloud
    281 
    282        USE arrays_3d
    283        USE cloud_parameters
    284        USE constants
    285        USE control_parameters
    286        USE indices
    287        USE kinds
    288 
    289        IMPLICIT NONE
    290 
    291        INTEGER(iwp) ::  i                 !:
    292        INTEGER(iwp) ::  j                 !:
    293        INTEGER(iwp) ::  k                 !:
    294  
    295        DO  i = nxl, nxr
    296           DO  j = nys, nyn
    297              DO  k = nzb_s_inner(j,i)+1, nzt
    298 
    299              ENDDO
    300           ENDDO
    301        ENDDO
    302 
    303     END SUBROUTINE sedimentation_cloud
    304 
    305 
    306     SUBROUTINE sedimentation_rain
    307 
    308        USE arrays_3d
    309        USE cloud_parameters
    310        USE constants
    311        USE control_parameters
    312        USE indices
    313        USE kinds
    314        USE statistics
    315 
    316        IMPLICIT NONE
    317 
    318        INTEGER(iwp) ::  i                 !:
    319        INTEGER(iwp) ::  j                 !:
    320        INTEGER(iwp) ::  k                 !:
    321  
    322        DO  i = nxl, nxr
    323           DO  j = nys, nyn
    324              DO  k = nzb_s_inner(j,i)+1, nzt
    325 
    326              ENDDO
    327           ENDDO
    328        ENDDO
    329 
    330     END SUBROUTINE sedimentation_rain
    331 
    332 
    333 !------------------------------------------------------------------------------!
    334 ! Call for grid point i,j
    335 !------------------------------------------------------------------------------!
    336 
    337     SUBROUTINE microphysics_control_ij( i, j )
    338 
    339        USE arrays_3d,                                                          &
    340            ONLY:  hyp, nc_1d,  nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,     &
    341                   qc_1d, qr, qr_1d, tend_nr, tend_pt, tend_q, tend_qr, zu
    342 
    343        USE cloud_parameters,                                                   &
    344            ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
    345 
    346        USE control_parameters,                                                 &
    347            ONLY:  drizzle, dt_3d, dt_micro, g, intermediate_timestep_count,    &
    348                   large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
    349                   rho_surface,surface_pressure
    350 
    351        USE indices,                                                            &
    352            ONLY:  nzb, nzt
    353 
    354        USE kinds
    355 
    356        USE statistics,                                                         &
    357            ONLY:  weight_pres
    358 
    359        IMPLICIT NONE
    360 
    361        INTEGER(iwp) ::  i                 !:
    362        INTEGER(iwp) ::  j                 !:
    363        INTEGER(iwp) ::  k                 !:
    364155
    365156       REAL(wp)     ::  t_surface         !:
    366157
    367        IF ( large_scale_forcing .AND. lsf_surf ) THEN
     158       IF ( large_scale_forcing  .AND. lsf_surf ) THEN
    368159!
    369160!--       Calculate:
    … …  
    374165          DO  k = nzb, nzt+1
    375166             hyp(k)    = surface_pressure * 100.0_wp * &
    376                          ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
     167                         ( ( t_surface - g / cp * zu(k) ) /                    &
     168                         t_surface )**(1.0_wp / 0.286_wp)
    377169             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
    378170             t_d_pt(k) = 1.0_wp / pt_d_t(k)
    … …  
    384176       ENDIF
    385177
    386 
    387        dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
    388 !
    389 !--    Adjust unrealistic values
    390        IF ( precipitation )  CALL adjust_cloud( i,j )
    391 !
    392 !--    Use 1-d arrays
    393        q_1d(:)  = q(:,j,i)
    394        pt_1d(:) = pt(:,j,i)
    395        qc_1d(:) = qc(:,j,i)
    396        nc_1d(:) = nc_const
    397        IF ( precipitation )  THEN
    398           qr_1d(:) = qr(:,j,i)
    399           nr_1d(:) = nr(:,j,i)
     178!
     179!--    Compute length of time step
     180       IF ( call_microphysics_at_all_substeps )  THEN
     181          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     182       ELSE
     183          dt_micro = dt_3d
    400184       ENDIF
     185
    401186!
    402187!--    Compute cloud physics
    403188       IF ( precipitation )  THEN
    404           CALL autoconversion( i,j )
    405           CALL accretion( i,j )
    406           CALL selfcollection_breakup( i,j )
    407           CALL evaporation_rain( i,j )
    408           CALL sedimentation_rain( i,j )
     189          CALL adjust_cloud
     190          CALL autoconversion
     191          CALL accretion
     192          CALL selfcollection_breakup
     193          CALL evaporation_rain
     194          CALL sedimentation_rain
    409195       ENDIF
    410196
    411        IF ( drizzle )  CALL sedimentation_cloud( i,j )
    412 !
    413 !--    Derive tendencies
    414        tend_q(:,j,i)  = ( q_1d(:) - q(:,j,i) ) / dt_micro
    415        tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
    416        IF ( precipitation )  THEN
    417           tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
    418           tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
    419        ENDIF
    420 
    421     END SUBROUTINE microphysics_control_ij
    422 
    423     SUBROUTINE adjust_cloud_ij( i, j )
     197       IF ( drizzle )  CALL sedimentation_cloud
     198
     199    END SUBROUTINE microphysics_control
     200
     201    SUBROUTINE adjust_cloud
    424202
    425203       USE arrays_3d,                                                          &
    … …  
    429207           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
    430208
     209       USE cpulog,                                                             &
     210           ONLY:  cpu_log, log_point_s
     211
    431212       USE indices,                                                            &
    432            ONLY:  nzb, nzb_s_inner, nzt
     213           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    433214
    434215       USE kinds
    … …  
    439220       INTEGER(iwp) ::  j                 !:
    440221       INTEGER(iwp) ::  k                 !:
    441 !
    442 !--    Adjust number of raindrops to avoid nonlinear effects in
    443 !--    sedimentation and evaporation of rain drops due to too small or
    444 !--    too big weights of rain drops (Stevens and Seifert, 2008).
    445 !--    The same procedure is applied to cloud droplets if they are determined
    446 !--    prognostically. 
    447        DO  k = nzb_s_inner(j,i)+1, nzt
    448 
    449           IF ( qr(k,j,i) <= eps_sb )  THEN
    450              qr(k,j,i) = 0.0_wp
    451              nr(k,j,i) = 0.0_wp
    452           ELSE
    453 !
    454 !--          Adjust number of raindrops to avoid nonlinear effects in
    455 !--          sedimentation and evaporation of rain drops due to too small or
    456 !--          too big weights of rain drops (Stevens and Seifert, 2008).
    457              IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    458                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
    459              ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
    460                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
    461              ENDIF
    462 
    463           ENDIF
    464 
     222
     223       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
     224
     225       DO  i = nxl, nxr
     226          DO  j = nys, nyn
     227             DO  k = nzb_s_inner(j,i)+1, nzt
     228                IF ( qr(k,j,i) <= eps_sb )  THEN
     229                   qr(k,j,i) = 0.0_wp
     230                   nr(k,j,i) = 0.0_wp
     231                ELSE
     232!
     233!--                Adjust number of raindrops to avoid nonlinear effects in
     234!--                sedimentation and evaporation of rain drops due to too small
     235!--                or too big weights of rain drops (Stevens and Seifert, 2008).
     236                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
     237                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
     238                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
     239                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
     240                   ENDIF
     241
     242                ENDIF
     243             ENDDO
     244          ENDDO
    465245       ENDDO
    466246
    467     END SUBROUTINE adjust_cloud_ij
    468 
    469 
    470     SUBROUTINE autoconversion_ij( i, j )
     247       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
     248
     249    END SUBROUTINE adjust_cloud
     250
     251
     252    SUBROUTINE autoconversion
    471253
    472254       USE arrays_3d,                                                          &
    473            ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
     255           ONLY:  diss, dzu, nr, qc, qr
    474256
    475257       USE cloud_parameters,                                                   &
    476258           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
    477                   c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
     259                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air,         &
     260                  nc_const, x0
    478261
    479262       USE control_parameters,                                                 &
    480263           ONLY:  dt_micro, rho_surface, turbulence
    481264
     265       USE cpulog,                                                             &
     266           ONLY:  cpu_log, log_point_s
     267
    482268       USE grid_variables,                                                     &
    483269           ONLY:  dx, dy
    484270
    485271       USE indices,                                                            &
    486            ONLY:  nzb, nzb_s_inner, nzt
     272           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    487273
    488274       USE kinds
    … …  
    496282       REAL(wp)     ::  alpha_cc          !:                   
    497283       REAL(wp)     ::  autocon           !:
    498        REAL(wp)     ::  epsilon           !:
     284       REAL(wp)     ::  dissipation       !:
    499285       REAL(wp)     ::  k_au              !:
    500286       REAL(wp)     ::  l_mix             !:
    … …  
    509295       REAL(wp)     ::  xc                !:
    510296
    511        k_au = k_cc / ( 20.0_wp * x0 )
    512 
    513        DO  k = nzb_s_inner(j,i)+1, nzt
    514 
    515           IF ( qc_1d(k) > eps_sb )  THEN
    516 !
    517 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    518 !--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
    519              tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
    520 !
    521 !--          Universal function for autoconversion process
    522 !--          (Seifert and Beheng, 2006):
    523              phi_au    = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
    524 !
    525 !--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
    526 !--          (Use constant nu_c = 1.0_wp instead?)
    527              nu_c      = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
    528 !
    529 !--          Mean weight of cloud droplets:
    530              xc = hyrho(k) * qc_1d(k) / nc_1d(k)
    531 !
    532 !--          Parameterized turbulence effects on autoconversion (Seifert,
    533 !--          Nuijens and Stevens, 2010)
    534              IF ( turbulence )  THEN
    535 !
    536 !--             Weight averaged radius of cloud droplets:
    537                 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
    538 
    539                 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
    540                 r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
    541                 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
    542 !
    543 !--             Mixing length (neglecting distance to ground and stratification)
    544                 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    545 !
    546 !--             Limit dissipation rate according to Seifert, Nuijens and
    547 !--             Stevens (2010)
    548                 epsilon = MIN( 0.06_wp, diss(k,j,i) )
    549 !
    550 !--             Compute Taylor-microscale Reynolds number:
    551                 re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *  &
    552                                SQRT( 15.0_wp / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp )
    553 !
    554 !--             The factor of 1.0E4 is needed to convert the dissipation rate
    555 !--             from m2 s-3 to cm2 s-3.
    556                 k_au = k_au * ( 1.0_wp +                                             &
    557                        epsilon * 1.0E4_wp * ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
    558                        ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /                 &
    559                        sigma_cc )**2 ) + beta_cc ) )
    560              ENDIF
    561 !
    562 !--          Autoconversion rate (Seifert and Beheng, 2006):
    563              autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /                &
    564                        ( nu_c + 1.0_wp )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp *   &
    565                        ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2.0_wp ) *        &
    566                        rho_surface
    567              autocon = MIN( autocon, qc_1d(k) / dt_micro )
    568 
    569              qr_1d(k) = qr_1d(k) + autocon * dt_micro
    570              qc_1d(k) = qc_1d(k) - autocon * dt_micro
    571              nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
    572 
    573           ENDIF
    574 
     297       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
     298
     299       DO  i = nxl, nxr
     300          DO  j = nys, nyn
     301             DO  k = nzb_s_inner(j,i)+1, nzt
     302
     303                IF ( qc(k,j,i) > eps_sb )  THEN
     304
     305                   k_au = k_cc / ( 20.0_wp * x0 )
     306!
     307!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
     308!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
     309                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
     310!
     311!--                Universal function for autoconversion process
     312!--                (Seifert and Beheng, 2006):
     313                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
     314                            ( 1.0_wp - tau_cloud**0.68_wp )**3
     315!
     316!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
     317!--                (Use constant nu_c = 1.0_wp instead?)
     318                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
     319!
     320!--                Mean weight of cloud droplets:
     321                   xc = hyrho(k) * qc(k,j,i) / nc_const
     322!
     323!--                Parameterized turbulence effects on autoconversion (Seifert,
     324!--                Nuijens and Stevens, 2010)
     325                   IF ( turbulence )  THEN
     326!
     327!--                   Weight averaged radius of cloud droplets:
     328                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     329
     330                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
     331                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
     332                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
     333!
     334!--                   Mixing length (neglecting distance to ground and
     335!--                   stratification)
     336                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
     337!
     338!--                   Limit dissipation rate according to Seifert, Nuijens and
     339!--                   Stevens (2010)
     340                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
     341!
     342!--                   Compute Taylor-microscale Reynolds number:
     343                      re_lambda = 6.0_wp / 11.0_wp *                           &
     344                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
     345                                  SQRT( 15.0_wp / kin_vis_air ) *              &
     346                                  dissipation**( 1.0_wp / 6.0_wp )
     347!
     348!--                   The factor of 1.0E4 is needed to convert the dissipation
     349!--                   rate from m2 s-3 to cm2 s-3.
     350                      k_au = k_au * ( 1.0_wp +                                 &
     351                                      dissipation * 1.0E4_wp *                 &
     352                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
     353                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &     
     354                                                                      r_cc ) / &
     355                                                        sigma_cc )**2          &
     356                                                      ) + beta_cc              &
     357                                      )                                        &
     358                                    )
     359                   ENDIF
     360!
     361!--                Autoconversion rate (Seifert and Beheng, 2006):
     362                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
     363                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
     364                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
     365                             rho_surface
     366                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
     367
     368                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
     369                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
     370                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
     371
     372                ENDIF
     373
     374             ENDDO
     375          ENDDO
    575376       ENDDO
    576377
    577     END SUBROUTINE autoconversion_ij
    578 
    579 
    580     SUBROUTINE accretion_ij( i, j )
     378       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
     379
     380    END SUBROUTINE autoconversion
     381
     382
     383    SUBROUTINE accretion
    581384
    582385       USE arrays_3d,                                                          &
    583            ONLY:  diss, qc_1d, qr_1d
     386           ONLY:  diss, qc, qr
    584387
    585388       USE cloud_parameters,                                                   &
    … …  
    589392           ONLY:  dt_micro, rho_surface, turbulence
    590393
     394       USE cpulog,                                                             &
     395           ONLY:  cpu_log, log_point_s
     396
    591397       USE indices,                                                            &
    592            ONLY:  nzb, nzb_s_inner, nzt
     398           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    593399
    594400       USE kinds
    … …  
    606412       REAL(wp)     ::  xc                !:
    607413
    608        DO  k = nzb_s_inner(j,i)+1, nzt
    609           IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
    610 !
    611 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    612              tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
    613 !
    614 !--          Universal function for accretion process
    615 !--          (Seifert and Beheng, 2001):
    616              phi_ac = tau_cloud / ( tau_cloud + 5.0E-5_wp )
    617              phi_ac = ( phi_ac**2.0_wp )**2.0_wp
    618 !
    619 !--          Parameterized turbulence effects on autoconversion (Seifert,
    620 !--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
    621 !--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
    622              IF ( turbulence )  THEN
    623                 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                             &
    624                                  MIN( 600.0_wp, diss(k,j,i) * 1.0E4_wp )**0.25_wp )
    625              ELSE
    626                 k_cr = k_cr0                       
    627              ENDIF
    628 !
    629 !--          Accretion rate (Seifert and Beheng, 2006):
    630              accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                 &
    631                     SQRT( rho_surface * hyrho(k) )
    632              accr = MIN( accr, qc_1d(k) / dt_micro )
    633 
    634              qr_1d(k) = qr_1d(k) + accr * dt_micro
    635              qc_1d(k) = qc_1d(k) - accr * dt_micro
    636 
    637           ENDIF
    638 
     414       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
     415
     416       DO  i = nxl, nxr
     417          DO  j = nys, nyn
     418             DO  k = nzb_s_inner(j,i)+1, nzt
     419
     420                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
     421!
     422!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
     423                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
     424!
     425!--                Universal function for accretion process (Seifert and
     426!--                Beheng, 2001):
     427                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
     428!
     429!--                Parameterized turbulence effects on autoconversion (Seifert,
     430!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
     431!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
     432                   IF ( turbulence )  THEN
     433                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
     434                                       MIN( 600.0_wp,                          &
     435                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
     436                                     )
     437                   ELSE
     438                      k_cr = k_cr0                       
     439                   ENDIF
     440!
     441!--                Accretion rate (Seifert and Beheng, 2006):
     442                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
     443                          SQRT( rho_surface * hyrho(k) )
     444                   accr = MIN( accr, qc(k,j,i) / dt_micro )
     445
     446                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro
     447                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
     448
     449                ENDIF
     450
     451             ENDDO
     452          ENDDO
    639453       ENDDO
    640454
    641     END SUBROUTINE accretion_ij
    642 
    643 
    644     SUBROUTINE selfcollection_breakup_ij( i, j )
     455       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
     456
     457    END SUBROUTINE accretion
     458
     459
     460    SUBROUTINE selfcollection_breakup
    645461
    646462       USE arrays_3d,                                                          &
    647            ONLY:  nr_1d, qr_1d
     463           ONLY:  nr, qr
    648464
    649465       USE cloud_parameters,                                                   &
    … …  
    653469           ONLY:  dt_micro, rho_surface
    654470
     471       USE cpulog,                                                             &
     472           ONLY:  cpu_log, log_point_s
     473
    655474       USE indices,                                                            &
    656            ONLY:  nzb, nzb_s_inner, nzt
     475           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    657476
    658477       USE kinds
    … …  
    669488       REAL(wp)     ::  selfcoll          !:
    670489
    671        DO  k = nzb_s_inner(j,i)+1, nzt
    672           IF ( qr_1d(k) > eps_sb )  THEN
    673 !
    674 !--          Selfcollection rate (Seifert and Beheng, 2001):
    675              selfcoll = k_rr * nr_1d(k) * qr_1d(k) *         &
    676                         SQRT( hyrho(k) * rho_surface )
    677 !
    678 !--          Weight averaged diameter of rain drops:
    679              dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
    680 !
    681 !--          Collisional breakup rate (Seifert, 2008):
    682              IF ( dr >= 0.3E-3_wp )  THEN
    683                 phi_br  = k_br * ( dr - 1.1E-3_wp )
    684                 breakup = selfcoll * ( phi_br + 1.0_wp )
    685              ELSE
    686                 breakup = 0.0_wp
    687              ENDIF
    688 
    689              selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
    690              nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
    691 
    692           ENDIF         
     490       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
     491
     492       DO  i = nxl, nxr
     493          DO  j = nys, nyn
     494             DO  k = nzb_s_inner(j,i)+1, nzt
     495                IF ( qr(k,j,i) > eps_sb )  THEN
     496!
     497!--                Selfcollection rate (Seifert and Beheng, 2001):
     498                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
     499                              SQRT( hyrho(k) * rho_surface )
     500!
     501!--                Weight averaged diameter of rain drops:
     502                   dr = ( hyrho(k) * qr(k,j,i) /                               &
     503                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
     504!
     505!--                Collisional breakup rate (Seifert, 2008):
     506                   IF ( dr >= 0.3E-3_wp )  THEN
     507                      phi_br  = k_br * ( dr - 1.1E-3_wp )
     508                      breakup = selfcoll * ( phi_br + 1.0_wp )
     509                   ELSE
     510                      breakup = 0.0_wp
     511                   ENDIF
     512
     513                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
     514                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
     515
     516                ENDIF         
     517             ENDDO
     518          ENDDO
    693519       ENDDO
    694520
    695     END SUBROUTINE selfcollection_breakup_ij
    696 
    697 
    698     SUBROUTINE evaporation_rain_ij( i, j )
     521       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
     522
     523    END SUBROUTINE selfcollection_breakup
     524
     525
     526    SUBROUTINE evaporation_rain
     527
    699528!
    700529!--    Evaporation of precipitable water. Condensation is neglected for
    … …  
    702531
    703532       USE arrays_3d,                                                          &
    704            ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
     533           ONLY:  hyp, nr, pt, q,  qc, qr
    705534
    706535       USE cloud_parameters,                                                   &
    … …  
    716545           ONLY:  dt_micro
    717546
     547       USE cpulog,                                                             &
     548           ONLY:  cpu_log, log_point_s
     549
    718550       USE indices,                                                            &
    719            ONLY:  nzb, nzb_s_inner, nzt
     551           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    720552
    721553       USE kinds
    … …  
    745577       REAL(wp)     ::  xr                !:
    746578
    747        DO  k = nzb_s_inner(j,i)+1, nzt
    748           IF ( qr_1d(k) > eps_sb )  THEN
    749 !
    750 !--          Actual liquid water temperature:
    751              t_l = t_d_pt(k) * pt_1d(k)
    752 !
    753 !--          Saturation vapor pressure at t_l:
    754              e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp ) )
    755 !
    756 !--          Computation of saturation humidity:
    757              q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
    758              alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
    759              q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
    760 !
    761 !--          Supersaturation:
    762              sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp )
    763 !
    764 !--          Actual temperature:
    765              temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
    766    
    767              g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /   &
    768                       ( thermal_conductivity_l * temp ) + r_v * temp / &
    769                       ( diff_coeff_l * e_s ) )
    770 !
    771 !--          Mean weight of rain drops
    772              xr = hyrho(k) * qr_1d(k) / nr_1d(k)
    773 !
    774 !--          Weight averaged diameter of rain drops:
    775              dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
    776 !
    777 !--          Compute ventilation factor and intercept parameter
    778 !--          (Seifert and Beheng, 2006; Seifert, 2008):
    779              IF ( ventilation_effect )  THEN
    780 !
    781 !--             Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
    782 !--             Stevens and Seifert, 2008):
    783                 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
    784 !
    785 !--             Slope parameter of gamma distribution (Seifert, 2008):
    786                 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *       &
    787                              ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
    788 
    789                 mu_r_2   = mu_r + 2.0_wp
    790                 mu_r_5d2 = mu_r + 2.5_wp
    791                 f_vent = a_vent * gamm( mu_r_2 ) *                            &
    792                          lambda_r**( -mu_r_2 ) +                              &
    793                          b_vent * schmidt_p_1d3 *                             &
    794                          SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
    795                          lambda_r**( -mu_r_5d2 ) *                            &
    796                          ( 1.0_wp - 0.5_wp * ( b_term / a_term ) *            &
    797                          ( lambda_r /                                         &
    798                          (       c_term + lambda_r ) )**mu_r_5d2 -            &
    799                                  0.125_wp * ( b_term / a_term )**2.0_wp *     &
    800                          ( lambda_r /                                         &
    801                          ( 2.0_wp * c_term + lambda_r ) )**mu_r_5d2 -         &
    802                                  0.0625_wp * ( b_term / a_term )**3.0_wp *    &
    803                          ( lambda_r /                                         &
    804                          ( 3.0_wp * c_term + lambda_r ) )**mu_r_5d2 -         &
    805                                  0.0390625_wp * ( b_term / a_term )**4.0_wp * &
    806                          ( lambda_r /                                         &
    807                          ( 4.0_wp * c_term + lambda_r ) )**mu_r_5d2 )
    808                 nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /             &
    809                          gamm( mu_r + 1.0_wp )
    810              ELSE
    811                 f_vent = 1.0_wp
    812                 nr_0   = nr_1d(k) * dr
    813              ENDIF
    814 !
    815 !--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
    816              evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /    &
    817                     hyrho(k)
    818 
    819              evap    = MAX( evap, -qr_1d(k) / dt_micro )
    820              evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
    821                             -nr_1d(k) / dt_micro )
    822 
    823              qr_1d(k) = qr_1d(k) + evap * dt_micro
    824              nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
    825           ENDIF         
    826 
     579       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
     580
     581       DO  i = nxl, nxr
     582          DO  j = nys, nyn
     583             DO  k = nzb_s_inner(j,i)+1, nzt
     584                IF ( qr(k,j,i) > eps_sb )  THEN
     585!
     586!--                Actual liquid water temperature:
     587                   t_l = t_d_pt(k) * pt(k,j,i)
     588!
     589!--                Saturation vapor pressure at t_l:
     590                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
     591                                          ( t_l - 35.86_wp )                   &
     592                                        )
     593!
     594!--                Computation of saturation humidity:
     595                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     596                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     597                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
     598                           ( 1.0_wp + alpha * q_s )
     599!
     600!--                Supersaturation:
     601                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
     602!
     603!--                Evaporation needs only to be calculated in subsaturated regions
     604                   IF ( sat < 0.0_wp )  THEN
     605!
     606!--                   Actual temperature:
     607                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
     608             
     609                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
     610                                          l_v / ( thermal_conductivity_l * temp ) &
     611                                          + r_v * temp / ( diff_coeff_l * e_s )   &
     612                                        )
     613!
     614!--                   Mean weight of rain drops
     615                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
     616!
     617!--                   Weight averaged diameter of rain drops:
     618                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
     619!
     620!--                   Compute ventilation factor and intercept parameter
     621!--                   (Seifert and Beheng, 2006; Seifert, 2008):
     622                      IF ( ventilation_effect )  THEN
     623!
     624!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
     625!--                      2005; Stevens and Seifert, 2008):
     626                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
     627                                                          ( dr - 1.4E-3_wp ) ) )
     628!
     629!--                      Slope parameter of gamma distribution (Seifert, 2008):
     630                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
     631                                      ( mu_r + 1.0_wp )                        &
     632                                    )**( 1.0_wp / 3.0_wp ) / dr
     633
     634                         mu_r_2   = mu_r + 2.0_wp
     635                         mu_r_5d2 = mu_r + 2.5_wp
     636
     637                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
     638                                  lambda_r**( -mu_r_2 ) + b_vent *              &
     639                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
     640                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
     641                                  ( 1.0_wp -                                    &
     642                                    0.5_wp * ( b_term / a_term ) *              &
     643                                    ( lambda_r / ( c_term + lambda_r )          &
     644                                    )**mu_r_5d2 -                               &
     645                                    0.125_wp * ( b_term / a_term )**2 *         &
     646                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
     647                                    )**mu_r_5d2 -                               &
     648                                    0.0625_wp * ( b_term / a_term )**3 *        &
     649                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
     650                                    )**mu_r_5d2 -                               &
     651                                    0.0390625_wp * ( b_term / a_term )**4 *     &
     652                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
     653                                    )**mu_r_5d2                                 &
     654                                  )
     655
     656                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
     657                                  gamm( mu_r + 1.0_wp )
     658                      ELSE
     659                         f_vent = 1.0_wp
     660                         nr_0   = nr(k,j,i) * dr
     661                      ENDIF
     662   !
     663   !--                Evaporation rate of rain water content (Seifert and
     664   !--                Beheng, 2006):
     665                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
     666                                hyrho(k)
     667                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
     668                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
     669                                     -nr(k,j,i) / dt_micro )
     670
     671                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
     672                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
     673
     674                   ENDIF
     675                ENDIF         
     676
     677             ENDDO
     678          ENDDO
    827679       ENDDO
    828680
    829     END SUBROUTINE evaporation_rain_ij
    830 
    831 
    832     SUBROUTINE sedimentation_cloud_ij( i, j )
     681       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
     682
     683    END SUBROUTINE evaporation_rain
     684
     685
     686    SUBROUTINE sedimentation_cloud
    833687
    834688       USE arrays_3d,                                                          &
    835            ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
     689           ONLY:  ddzu, dzu, pt, q, qc
    836690
    837691       USE cloud_parameters,                                                   &
    838            ONLY:  eps_sb, hyrho, k_st, l_d_cp, prr, pt_d_t, rho_l, sigma_gc
     692           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
    839693
    840694       USE constants,                                                          &
    … …  
    844698           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
    845699
     700       USE cpulog,                                                             &
     701           ONLY:  cpu_log, log_point_s
     702
    846703       USE indices,                                                            &
    847            ONLY:  nzb, nzb_s_inner, nzt
     704           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    848705
    849706       USE kinds
    … …  
    855712       INTEGER(iwp) ::  k                 !:
    856713
    857        REAL(wp)     ::  sed_qc_const      !:
    858 
    859 
    860        REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc
    861 
    862 !
    863 !--    Sedimentation of cloud droplets (Heus et al., 2010):
    864        sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ))**( 2.0_wp / 3.0_wp ) *   &
    865                      EXP( 5.0_wp * LOG( sigma_gc )**2 )
    866 
     714       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !:
     715
     716       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
     717
     718!
     719!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
    867720       sed_qc(nzt+1) = 0.0_wp
    868721
    869        DO  k = nzt, nzb_s_inner(j,i)+1, -1
    870           IF ( qc_1d(k) > eps_sb )  THEN
    871              sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * &
    872                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
    873           ELSE
    874              sed_qc(k) = 0.0_wp
    875           ENDIF
    876 
    877           sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /     &
    878                                       dt_micro + sed_qc(k+1) )
    879 
    880           q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /  &
    881                                 hyrho(k) * dt_micro
    882           qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
    883                                 hyrho(k) * dt_micro
    884           pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
    885                                 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
    886 
     722       DO  i = nxl, nxr
     723          DO  j = nys, nyn
     724             DO  k = nzt, nzb_s_inner(j,i)+1, -1
     725
     726                IF ( qc(k,j,i) > eps_sb )  THEN
     727                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
     728                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
     729                ELSE
     730                   sed_qc(k) = 0.0_wp
     731                ENDIF
     732
     733                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
     734                                            dt_micro + sed_qc(k+1)             &
     735                               )
     736
     737                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
     738                                        ddzu(k+1) / hyrho(k) * dt_micro
     739                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
     740                                        ddzu(k+1) / hyrho(k) * dt_micro
     741                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
     742                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
     743                                        pt_d_t(k) * dt_micro
     744
     745             ENDDO
     746          ENDDO
    887747       ENDDO
    888748
    889     END SUBROUTINE sedimentation_cloud_ij
    890 
    891 
    892     SUBROUTINE sedimentation_rain_ij( i, j )
     749       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
     750
     751    END SUBROUTINE sedimentation_cloud
     752
     753
     754    SUBROUTINE sedimentation_rain
    893755
    894756       USE arrays_3d,                                                          &
    895            ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
     757           ONLY:  ddzu, dzu, nr, pt, q, qr
    896758
    897759       USE cloud_parameters,                                                   &
    … …  
    901763
    902764       USE control_parameters,                                                 &
    903            ONLY:  dt_do2d_xy, dt_micro, dt_3d, intermediate_timestep_count,    &
     765           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
     766                  dt_3d, intermediate_timestep_count,                          &
    904767                  intermediate_timestep_count_max,                             &
    905768                  precipitation_amount_interval, time_do2d_xy
    906769
     770       USE cpulog,                                                             &
     771           ONLY:  cpu_log, log_point_s
     772
    907773       USE indices,                                                            &
    908            ONLY:  nzb, nzb_s_inner, nzt
     774           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
    909775
    910776       USE kinds
    … …  
    931797       REAL(wp)     ::  z_run                      !:
    932798
    933       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr      !:
    934       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr      !:
    935       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr      !:
    936       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr      !:
    937       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !:
    938       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !:
    939       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr    !:
    940       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr    !:
    941       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr      !:
    942       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr      !:
    943 
    944 
     799       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !:
     800       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !:
     801       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !:
     802       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !:
     803       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !:
     804       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !:
     805       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !:
     806       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !:
     807       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !:
     808       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !:
     809
     810       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
    945811!
    946812!--    Computation of sedimentation flux. Implementation according to Stevens
    947 !--    and Seifert (2008).
     813!--    and Seifert (2008). Code is based on UCLA-LES.
     814       IF ( intermediate_timestep_count == 1 )  prr(:,:,:) = 0.0_wp
     815!
     816!--    Compute velocities
     817       DO  i = nxl, nxr
     818          DO  j = nys, nyn
     819             DO  k = nzb_s_inner(j,i)+1, nzt
     820                IF ( qr(k,j,i) > eps_sb )  THEN
     821!
     822!--                Weight averaged diameter of rain drops:
     823                   dr = ( hyrho(k) * qr(k,j,i) /                               &
     824                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
     825!
     826!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
     827!--                Stevens and Seifert, 2008):
     828                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
     829                                                     ( dr - 1.4E-3_wp ) ) )
     830!
     831!--                Slope parameter of gamma distribution (Seifert, 2008):
     832                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
     833                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
     834
     835                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
     836                                               a_term - b_term * ( 1.0_wp +    &
     837                                                  c_term /                     &
     838                                                  lambda_r )**( -1.0_wp *      &
     839                                                  ( mu_r + 1.0_wp ) )          &
     840                                              )                                &
     841                                )
     842
     843                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
     844                                               a_term - b_term * ( 1.0_wp +    &
     845                                                  c_term /                     &
     846                                                  lambda_r )**( -1.0_wp *      &
     847                                                  ( mu_r + 4.0_wp ) )          &
     848                                             )                                 &
     849                                )
     850                ELSE
     851                   w_nr(k) = 0.0_wp
     852                   w_qr(k) = 0.0_wp
     853                ENDIF
     854             ENDDO
     855!
     856!--          Adjust boundary values
     857             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
     858             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
     859             w_nr(nzt+1) = 0.0_wp
     860             w_qr(nzt+1) = 0.0_wp
     861!
     862!--          Compute Courant number
     863             DO  k = nzb_s_inner(j,i)+1, nzt
     864                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
     865                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
     866                          dt_micro * ddzu(k)
     867                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
     868                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
     869                          dt_micro * ddzu(k)
     870             ENDDO     
     871!
     872!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
     873             IF ( limiter_sedimentation )  THEN
     874
     875                DO k = nzb_s_inner(j,i)+1, nzt
     876                   d_mean = 0.5_wp * ( qr(k+1,j,i) + qr(k-1,j,i) )
     877                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
     878                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
     879
     880                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
     881                                                              2.0_wp * d_max,  &
     882                                                              ABS( d_mean ) )
     883
     884                   d_mean = 0.5_wp * ( nr(k+1,j,i) + nr(k-1,j,i) )
     885                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
     886                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
     887
     888                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
     889                                                              2.0_wp * d_max,  &
     890                                                              ABS( d_mean ) )
     891                ENDDO
     892
     893             ELSE
     894
     895                nr_slope = 0.0_wp
     896                qr_slope = 0.0_wp
     897
     898             ENDIF
     899
     900             sed_nr(nzt+1) = 0.0_wp
     901             sed_qr(nzt+1) = 0.0_wp
     902!
     903!--          Compute sedimentation flux
     904             DO  k = nzt, nzb_s_inner(j,i)+1, -1
     905!
     906!--             Sum up all rain drop number densities which contribute to the flux
     907!--             through k-1/2
     908                flux  = 0.0_wp
     909                z_run = 0.0_wp ! height above z(k)
     910                k_run = k
     911                c_run = MIN( 1.0_wp, c_nr(k) )
     912                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
     913                   flux  = flux + hyrho(k_run) *                               &
     914                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
     915                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
     916                   z_run = z_run + dzu(k_run)
     917                   k_run = k_run + 1
     918                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
     919                ENDDO
     920!
     921!--             It is not allowed to sediment more rain drop number density than
     922!--             available
     923                flux = MIN( flux,                                              &
     924                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
     925                            dt_micro                                           &
     926                          )
     927
     928                sed_nr(k) = flux / dt_micro
     929                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
     930                                        ddzu(k+1) / hyrho(k) * dt_micro
     931!
     932!--             Sum up all rain water content which contributes to the flux
     933!--             through k-1/2
     934                flux  = 0.0_wp
     935                z_run = 0.0_wp ! height above z(k)
     936                k_run = k
     937                c_run = MIN( 1.0_wp, c_qr(k) )
     938
     939                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
     940
     941                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
     942                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
     943                                  0.5_wp ) * c_run * dzu(k_run)
     944                   z_run = z_run + dzu(k_run)
     945                   k_run = k_run + 1
     946                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
     947
     948                ENDDO
     949!
     950!--             It is not allowed to sediment more rain water content than
     951!--             available
     952                flux = MIN( flux,                                              &
     953                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
     954                            dt_micro                                           &
     955                          )
     956
     957                sed_qr(k) = flux / dt_micro
     958
     959                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
     960                                        ddzu(k+1) / hyrho(k) * dt_micro
     961                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
     962                                        ddzu(k+1) / hyrho(k) * dt_micro
     963                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
     964                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
     965                                        pt_d_t(k) * dt_micro
     966!
     967!--             Compute the rain rate
     968                IF ( call_microphysics_at_all_substeps )  THEN
     969                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *            &
     970                                weight_substep(intermediate_timestep_count)
     971                ELSE
     972                   prr(k,j,i) = sed_qr(k) / hyrho(k)
     973                ENDIF
     974
     975             ENDDO
     976          ENDDO
     977       ENDDO
     978
     979!
     980!--    Precipitation amount
     981       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
     982            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
     983                   precipitation_amount_interval )  THEN
     984          DO  i = nxl, nxr
     985             DO  j = nys, nyn
     986                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
     987                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
     988                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
     989             ENDDO
     990          ENDDO
     991       ENDIF
     992
     993       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
     994
     995    END SUBROUTINE sedimentation_rain
     996
     997
     998!------------------------------------------------------------------------------!
     999! Call for grid point i,j
     1000!------------------------------------------------------------------------------!
     1001
     1002    SUBROUTINE microphysics_control_ij( i, j )
     1003
     1004       USE arrays_3d,                                                          &
     1005           ONLY:  hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,      &
     1006                  qc_1d, qr, qr_1d, zu
     1007
     1008       USE cloud_parameters,                                                   &
     1009           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
     1010
     1011       USE control_parameters,                                                 &
     1012           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
     1013                  g, intermediate_timestep_count, large_scale_forcing,         &
     1014                  lsf_surf, precipitation, pt_surface,                         &
     1015                  rho_surface,surface_pressure
     1016
     1017       USE indices,                                                            &
     1018           ONLY:  nzb, nzt
     1019
     1020       USE kinds
     1021
     1022       USE statistics,                                                         &
     1023           ONLY:  weight_pres
     1024
     1025       IMPLICIT NONE
     1026
     1027       INTEGER(iwp) ::  i                 !:
     1028       INTEGER(iwp) ::  j                 !:
     1029       INTEGER(iwp) ::  k                 !:
     1030
     1031       REAL(wp)     ::  t_surface         !:
     1032
     1033       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
     1034!
     1035!--       Calculate:
     1036!--       pt / t : ratio of potential and actual temperature (pt_d_t)
     1037!--       t / pt : ratio of actual and potential temperature (t_d_pt)
     1038!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
     1039          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
     1040          DO  k = nzb, nzt+1
     1041             hyp(k)    = surface_pressure * 100.0_wp * &
     1042                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
     1043             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
     1044             t_d_pt(k) = 1.0_wp / pt_d_t(k)
     1045             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     1046          ENDDO
     1047!
     1048!--       Compute reference density
     1049          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
     1050       ENDIF
     1051
     1052!
     1053!--    Compute length of time step
     1054       IF ( call_microphysics_at_all_substeps )  THEN
     1055          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     1056       ELSE
     1057          dt_micro = dt_3d
     1058       ENDIF
     1059
     1060!
     1061!--    Use 1d arrays
     1062       q_1d(:)  = q(:,j,i)
     1063       pt_1d(:) = pt(:,j,i)
     1064       qc_1d(:) = qc(:,j,i)
     1065       nc_1d(:) = nc_const
     1066       IF ( precipitation )  THEN
     1067          qr_1d(:) = qr(:,j,i)
     1068          nr_1d(:) = nr(:,j,i)
     1069       ENDIF
     1070
     1071!
     1072!--    Compute cloud physics
     1073       IF ( precipitation )  THEN
     1074          CALL adjust_cloud( i,j )
     1075          CALL autoconversion( i,j )
     1076          CALL accretion( i,j )
     1077          CALL selfcollection_breakup( i,j )
     1078          CALL evaporation_rain( i,j )
     1079          CALL sedimentation_rain( i,j )
     1080       ENDIF
     1081
     1082       IF ( drizzle )  CALL sedimentation_cloud( i,j )
     1083
     1084!
     1085!--    Store results on the 3d arrays
     1086       q(:,j,i)  = q_1d(:)
     1087       pt(:,j,i) = pt_1d(:)
     1088       IF ( precipitation )  THEN
     1089          qr(:,j,i) = qr_1d(:)
     1090          nr(:,j,i) = nr_1d(:)
     1091       ENDIF
     1092
     1093    END SUBROUTINE microphysics_control_ij
     1094
     1095    SUBROUTINE adjust_cloud_ij( i, j )
     1096
     1097       USE arrays_3d,                                                          &
     1098           ONLY:  qr_1d, nr_1d
     1099
     1100       USE cloud_parameters,                                                   &
     1101           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
     1102
     1103       USE indices,                                                            &
     1104           ONLY:  nzb, nzb_s_inner, nzt
     1105
     1106       USE kinds
     1107
     1108       IMPLICIT NONE
     1109
     1110       INTEGER(iwp) ::  i                 !:
     1111       INTEGER(iwp) ::  j                 !:
     1112       INTEGER(iwp) ::  k                 !:
     1113!
     1114!--    Adjust number of raindrops to avoid nonlinear effects in
     1115!--    sedimentation and evaporation of rain drops due to too small or
     1116!--    too big weights of rain drops (Stevens and Seifert, 2008).
     1117!--    The same procedure is applied to cloud droplets if they are determined
     1118!--    prognostically. 
     1119       DO  k = nzb_s_inner(j,i)+1, nzt
     1120
     1121          IF ( qr_1d(k) <= eps_sb )  THEN
     1122             qr_1d(k) = 0.0_wp
     1123             nr_1d(k) = 0.0_wp
     1124          ELSE
     1125!
     1126!--          Adjust number of raindrops to avoid nonlinear effects in
     1127!--          sedimentation and evaporation of rain drops due to too small or
     1128!--          too big weights of rain drops (Stevens and Seifert, 2008).
     1129             IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
     1130                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
     1131             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
     1132                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
     1133             ENDIF
     1134
     1135          ENDIF
     1136
     1137       ENDDO
     1138
     1139    END SUBROUTINE adjust_cloud_ij
     1140
     1141
     1142    SUBROUTINE autoconversion_ij( i, j )
     1143
     1144       USE arrays_3d,                                                          &
     1145           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
     1146
     1147       USE cloud_parameters,                                                   &
     1148           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
     1149                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
     1150
     1151       USE control_parameters,                                                 &
     1152           ONLY:  dt_micro, rho_surface, turbulence
     1153
     1154       USE grid_variables,                                                     &
     1155           ONLY:  dx, dy
     1156
     1157       USE indices,                                                            &
     1158           ONLY:  nzb, nzb_s_inner, nzt
     1159
     1160       USE kinds
     1161
     1162       IMPLICIT NONE
     1163
     1164       INTEGER(iwp) ::  i                 !:
     1165       INTEGER(iwp) ::  j                 !:
     1166       INTEGER(iwp) ::  k                 !:
     1167
     1168       REAL(wp)     ::  alpha_cc          !:                   
     1169       REAL(wp)     ::  autocon           !:
     1170       REAL(wp)     ::  dissipation       !:
     1171       REAL(wp)     ::  k_au              !:
     1172       REAL(wp)     ::  l_mix             !:
     1173       REAL(wp)     ::  nu_c              !:
     1174       REAL(wp)     ::  phi_au            !:
     1175       REAL(wp)     ::  r_cc              !:
     1176       REAL(wp)     ::  rc                !:
     1177       REAL(wp)     ::  re_lambda         !:
     1178       REAL(wp)     ::  selfcoll          !:
     1179       REAL(wp)     ::  sigma_cc          !:
     1180       REAL(wp)     ::  tau_cloud         !:
     1181       REAL(wp)     ::  xc                !:
     1182
     1183       DO  k = nzb_s_inner(j,i)+1, nzt
     1184
     1185          IF ( qc_1d(k) > eps_sb )  THEN
     1186
     1187             k_au = k_cc / ( 20.0_wp * x0 )
     1188!
     1189!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
     1190!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
     1191             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
     1192!
     1193!--          Universal function for autoconversion process
     1194!--          (Seifert and Beheng, 2006):
     1195             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
     1196!
     1197!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
     1198!--          (Use constant nu_c = 1.0_wp instead?)
     1199             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
     1200!
     1201!--          Mean weight of cloud droplets:
     1202             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
     1203!
     1204!--          Parameterized turbulence effects on autoconversion (Seifert,
     1205!--          Nuijens and Stevens, 2010)
     1206             IF ( turbulence )  THEN
     1207!
     1208!--             Weight averaged radius of cloud droplets:
     1209                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     1210
     1211                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
     1212                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
     1213                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
     1214!
     1215!--             Mixing length (neglecting distance to ground and stratification)
     1216                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
     1217!
     1218!--             Limit dissipation rate according to Seifert, Nuijens and
     1219!--             Stevens (2010)
     1220                dissipation = MIN( 0.06_wp, diss(k,j,i) )
     1221!
     1222!--             Compute Taylor-microscale Reynolds number:
     1223                re_lambda = 6.0_wp / 11.0_wp *                                 &
     1224                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
     1225                            SQRT( 15.0_wp / kin_vis_air ) *                    &
     1226                            dissipation**( 1.0_wp / 6.0_wp )
     1227!
     1228!--             The factor of 1.0E4 is needed to convert the dissipation rate
     1229!--             from m2 s-3 to cm2 s-3.
     1230                k_au = k_au * ( 1.0_wp +                                       &
     1231                                dissipation * 1.0E4_wp *                       &
     1232                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
     1233                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
     1234                                                  sigma_cc )**2                &
     1235                                                ) + beta_cc                    &
     1236                                )                                              &
     1237                              )
     1238             ENDIF
     1239!
     1240!--          Autoconversion rate (Seifert and Beheng, 2006):
     1241             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
     1242                       ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 *            &
     1243                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
     1244                       rho_surface
     1245             autocon = MIN( autocon, qc_1d(k) / dt_micro )
     1246
     1247             qr_1d(k) = qr_1d(k) + autocon * dt_micro
     1248             qc_1d(k) = qc_1d(k) - autocon * dt_micro
     1249             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
     1250
     1251          ENDIF
     1252
     1253       ENDDO
     1254
     1255    END SUBROUTINE autoconversion_ij
     1256
     1257
     1258    SUBROUTINE accretion_ij( i, j )
     1259
     1260       USE arrays_3d,                                                          &
     1261           ONLY:  diss, qc_1d, qr_1d
     1262
     1263       USE cloud_parameters,                                                   &
     1264           ONLY:  eps_sb, hyrho, k_cr0
     1265
     1266       USE control_parameters,                                                 &
     1267           ONLY:  dt_micro, rho_surface, turbulence
     1268
     1269       USE indices,                                                            &
     1270           ONLY:  nzb, nzb_s_inner, nzt
     1271
     1272       USE kinds
     1273
     1274       IMPLICIT NONE
     1275
     1276       INTEGER(iwp) ::  i                 !:
     1277       INTEGER(iwp) ::  j                 !:
     1278       INTEGER(iwp) ::  k                 !:
     1279
     1280       REAL(wp)     ::  accr              !:
     1281       REAL(wp)     ::  k_cr              !:
     1282       REAL(wp)     ::  phi_ac            !:
     1283       REAL(wp)     ::  tau_cloud         !:
     1284       REAL(wp)     ::  xc                !:
     1285
     1286       DO  k = nzb_s_inner(j,i)+1, nzt
     1287          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
     1288!
     1289!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
     1290             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
     1291!
     1292!--          Universal function for accretion process
     1293!--          (Seifert and Beheng, 2001):
     1294             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
     1295!
     1296!--          Parameterized turbulence effects on autoconversion (Seifert,
     1297!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
     1298!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
     1299             IF ( turbulence )  THEN
     1300                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
     1301                                 MIN( 600.0_wp,                                &
     1302                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
     1303                               )
     1304             ELSE
     1305                k_cr = k_cr0                       
     1306             ENDIF
     1307!
     1308!--          Accretion rate (Seifert and Beheng, 2006):
     1309             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
     1310             accr = MIN( accr, qc_1d(k) / dt_micro )
     1311
     1312             qr_1d(k) = qr_1d(k) + accr * dt_micro
     1313             qc_1d(k) = qc_1d(k) - accr * dt_micro
     1314
     1315          ENDIF
     1316
     1317       ENDDO
     1318
     1319    END SUBROUTINE accretion_ij
     1320
     1321
     1322    SUBROUTINE selfcollection_breakup_ij( i, j )
     1323
     1324       USE arrays_3d,                                                          &
     1325           ONLY:  nr_1d, qr_1d
     1326
     1327       USE cloud_parameters,                                                   &
     1328           ONLY:  dpirho_l, eps_sb, hyrho, k_br, k_rr
     1329
     1330       USE control_parameters,                                                 &
     1331           ONLY:  dt_micro, rho_surface
     1332
     1333       USE indices,                                                            &
     1334           ONLY:  nzb, nzb_s_inner, nzt
     1335
     1336       USE kinds
     1337   
     1338       IMPLICIT NONE
     1339
     1340       INTEGER(iwp) ::  i                 !:
     1341       INTEGER(iwp) ::  j                 !:
     1342       INTEGER(iwp) ::  k                 !:
     1343
     1344       REAL(wp)     ::  breakup           !:
     1345       REAL(wp)     ::  dr                !:
     1346       REAL(wp)     ::  phi_br            !:
     1347       REAL(wp)     ::  selfcoll          !:
     1348
     1349       DO  k = nzb_s_inner(j,i)+1, nzt
     1350          IF ( qr_1d(k) > eps_sb )  THEN
     1351!
     1352!--          Selfcollection rate (Seifert and Beheng, 2001):
     1353             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
     1354!
     1355!--          Weight averaged diameter of rain drops:
     1356             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
     1357!
     1358!--          Collisional breakup rate (Seifert, 2008):
     1359             IF ( dr >= 0.3E-3_wp )  THEN
     1360                phi_br  = k_br * ( dr - 1.1E-3_wp )
     1361                breakup = selfcoll * ( phi_br + 1.0_wp )
     1362             ELSE
     1363                breakup = 0.0_wp
     1364             ENDIF
     1365
     1366             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
     1367             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
     1368
     1369          ENDIF         
     1370       ENDDO
     1371
     1372    END SUBROUTINE selfcollection_breakup_ij
     1373
     1374
     1375    SUBROUTINE evaporation_rain_ij( i, j )
     1376!
     1377!--    Evaporation of precipitable water. Condensation is neglected for
     1378!--    precipitable water.
     1379
     1380       USE arrays_3d,                                                          &
     1381           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
     1382
     1383       USE cloud_parameters,                                                   &
     1384           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
     1385                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
     1386                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
     1387                  t_d_pt, ventilation_effect
     1388
     1389       USE constants,                                                          &
     1390           ONLY:  pi
     1391
     1392       USE control_parameters,                                                 &
     1393           ONLY:  dt_micro
     1394
     1395       USE indices,                                                            &
     1396           ONLY:  nzb, nzb_s_inner, nzt
     1397
     1398       USE kinds
     1399
     1400       IMPLICIT NONE
     1401
     1402       INTEGER(iwp) ::  i                 !:
     1403       INTEGER(iwp) ::  j                 !:
     1404       INTEGER(iwp) ::  k                 !:
     1405
     1406       REAL(wp)     ::  alpha             !:
     1407       REAL(wp)     ::  dr                !:
     1408       REAL(wp)     ::  e_s               !:
     1409       REAL(wp)     ::  evap              !:
     1410       REAL(wp)     ::  evap_nr           !:
     1411       REAL(wp)     ::  f_vent            !:
     1412       REAL(wp)     ::  g_evap            !:
     1413       REAL(wp)     ::  lambda_r          !:
     1414       REAL(wp)     ::  mu_r              !:
     1415       REAL(wp)     ::  mu_r_2            !:
     1416       REAL(wp)     ::  mu_r_5d2          !:
     1417       REAL(wp)     ::  nr_0              !:
     1418       REAL(wp)     ::  q_s               !:
     1419       REAL(wp)     ::  sat               !:
     1420       REAL(wp)     ::  t_l               !:
     1421       REAL(wp)     ::  temp              !:
     1422       REAL(wp)     ::  xr                !:
     1423
     1424       DO  k = nzb_s_inner(j,i)+1, nzt
     1425          IF ( qr_1d(k) > eps_sb )  THEN
     1426!
     1427!--          Actual liquid water temperature:
     1428             t_l = t_d_pt(k) * pt_1d(k)
     1429!
     1430!--          Saturation vapor pressure at t_l:
     1431             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
     1432                                    ( t_l - 35.86_wp )                         &
     1433                                  )
     1434!
     1435!--          Computation of saturation humidity:
     1436             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     1437             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     1438             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
     1439!
     1440!--          Supersaturation:
     1441             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
     1442!
     1443!--          Evaporation needs only to be calculated in subsaturated regions
     1444             IF ( sat < 0.0_wp )  THEN
     1445!
     1446!--             Actual temperature:
     1447                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
     1448       
     1449                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
     1450                                    ( thermal_conductivity_l * temp ) +        &
     1451                                    r_v * temp / ( diff_coeff_l * e_s )        &
     1452                                  )
     1453!
     1454!--             Mean weight of rain drops
     1455                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
     1456!
     1457!--             Weight averaged diameter of rain drops:
     1458                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
     1459!
     1460!--             Compute ventilation factor and intercept parameter
     1461!--             (Seifert and Beheng, 2006; Seifert, 2008):
     1462                IF ( ventilation_effect )  THEN
     1463!
     1464!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
     1465!--                Stevens and Seifert, 2008):
     1466                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
     1467!
     1468!--                Slope parameter of gamma distribution (Seifert, 2008):
     1469                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
     1470                                ( mu_r + 1.0_wp )                              &
     1471                              )**( 1.0_wp / 3.0_wp ) / dr
     1472
     1473                   mu_r_2   = mu_r + 2.0_wp
     1474                   mu_r_5d2 = mu_r + 2.5_wp
     1475
     1476                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
     1477                            b_vent * schmidt_p_1d3 *                           &
     1478                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
     1479                            lambda_r**( -mu_r_5d2 ) *                          &
     1480                            ( 1.0_wp -                                         &
     1481                              0.5_wp * ( b_term / a_term ) *                   &
     1482                              ( lambda_r / ( c_term + lambda_r )               &
     1483                              )**mu_r_5d2 -                                    &
     1484                              0.125_wp * ( b_term / a_term )**2 *              &
     1485                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
     1486                              )**mu_r_5d2 -                                    &
     1487                              0.0625_wp * ( b_term / a_term )**3 *             &
     1488                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
     1489                              )**mu_r_5d2 -                                    &
     1490                              0.0390625_wp * ( b_term / a_term )**4 *          &
     1491                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
     1492                              )**mu_r_5d2                                      &
     1493                            )
     1494
     1495                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
     1496                            gamm( mu_r + 1.0_wp )
     1497                ELSE
     1498                   f_vent = 1.0_wp
     1499                   nr_0   = nr_1d(k) * dr
     1500                ENDIF
     1501!
     1502!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
     1503                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
     1504                evap    = MAX( evap, -qr_1d(k) / dt_micro )
     1505                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
     1506                               -nr_1d(k) / dt_micro )
     1507
     1508                qr_1d(k) = qr_1d(k) + evap    * dt_micro
     1509                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
     1510
     1511             ENDIF
     1512          ENDIF         
     1513
     1514       ENDDO
     1515
     1516    END SUBROUTINE evaporation_rain_ij
     1517
     1518
     1519    SUBROUTINE sedimentation_cloud_ij( i, j )
     1520
     1521       USE arrays_3d,                                                          &
     1522           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
     1523
     1524       USE cloud_parameters,                                                   &
     1525           ONLY:  eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
     1526
     1527       USE constants,                                                          &
     1528           ONLY:  pi
     1529
     1530       USE control_parameters,                                                 &
     1531           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
     1532
     1533       USE indices,                                                            &
     1534           ONLY:  nzb, nzb_s_inner, nzt
     1535
     1536       USE kinds
     1537       
     1538       IMPLICIT NONE
     1539
     1540       INTEGER(iwp) ::  i                 !:
     1541       INTEGER(iwp) ::  j                 !:
     1542       INTEGER(iwp) ::  k                 !:
     1543
     1544       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !:
     1545
     1546!
     1547!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
     1548       sed_qc(nzt+1) = 0.0_wp
     1549
     1550       DO  k = nzt, nzb_s_inner(j,i)+1, -1
     1551          IF ( qc_1d(k) > eps_sb )  THEN
     1552             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
     1553                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
     1554          ELSE
     1555             sed_qc(k) = 0.0_wp
     1556          ENDIF
     1557
     1558          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
     1559                                      dt_micro + sed_qc(k+1)                   &
     1560                         )
     1561
     1562          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     1563                                hyrho(k) * dt_micro
     1564          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     1565                                hyrho(k) * dt_micro
     1566          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     1567                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
     1568
     1569       ENDDO
     1570
     1571    END SUBROUTINE sedimentation_cloud_ij
     1572
     1573
     1574    SUBROUTINE sedimentation_rain_ij( i, j )
     1575
     1576       USE arrays_3d,                                                          &
     1577           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
     1578
     1579       USE cloud_parameters,                                                   &
     1580           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
     1581                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
     1582                  pt_d_t, stp
     1583
     1584       USE control_parameters,                                                 &
     1585           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
     1586                  dt_3d, intermediate_timestep_count,                          &
     1587                  intermediate_timestep_count_max,                             &
     1588                  precipitation_amount_interval, time_do2d_xy
     1589
     1590       USE indices,                                                            &
     1591           ONLY:  nzb, nzb_s_inner, nzt
     1592
     1593       USE kinds
     1594
     1595       USE statistics,                                                         &
     1596           ONLY:  weight_substep
     1597       
     1598       IMPLICIT NONE
     1599
     1600       INTEGER(iwp) ::  i                          !:
     1601       INTEGER(iwp) ::  j                          !:
     1602       INTEGER(iwp) ::  k                          !:
     1603       INTEGER(iwp) ::  k_run                      !:
     1604
     1605       REAL(wp)     ::  c_run                      !:
     1606       REAL(wp)     ::  d_max                      !:
     1607       REAL(wp)     ::  d_mean                     !:
     1608       REAL(wp)     ::  d_min                      !:
     1609       REAL(wp)     ::  dr                         !:
     1610       REAL(wp)     ::  dt_sedi                    !:
     1611       REAL(wp)     ::  flux                       !:
     1612       REAL(wp)     ::  lambda_r                   !:
     1613       REAL(wp)     ::  mu_r                       !:
     1614       REAL(wp)     ::  z_run                      !:
     1615
     1616       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !:
     1617       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !:
     1618       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !:
     1619       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_qr     !:
     1620       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !:
     1621       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !:
     1622       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !:
     1623       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !:
     1624       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !:
     1625       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !:
     1626
     1627
     1628!
     1629!--    Computation of sedimentation flux. Implementation according to Stevens
     1630!--    and Seifert (2008). Code is based on UCLA-LES.
    9481631       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
    9491632!
    … …  
    9601643!
    9611644!--          Slope parameter of gamma distribution (Seifert, 2008):
    962              lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *          &
    963                         ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
    964 
    965              w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp +    &
    966                        c_term / lambda_r )**( -1.0_wp * ( mu_r + 1.0_wp ) ) ) )
    967              w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp +    &
    968                        c_term / lambda_r )**( -1.0_wp * ( mu_r + 4.0_wp ) ) ) )
     1645             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
     1646                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
     1647
     1648             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
     1649                                         a_term - b_term * ( 1.0_wp +          &
     1650                                            c_term / lambda_r )**( -1.0_wp *   &
     1651                                            ( mu_r + 1.0_wp ) )                &
     1652                                        )                                      &
     1653                          )
     1654             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
     1655                                         a_term - b_term * ( 1.0_wp +          &
     1656                                            c_term / lambda_r )**( -1.0_wp *   &
     1657                                            ( mu_r + 4.0_wp ) )                &
     1658                                       )                                       &
     1659                          )
    9691660          ELSE
    9701661             w_nr(k) = 0.0_wp
    … …  
    9811672!--    Compute Courant number
    9821673       DO  k = nzb_s_inner(j,i)+1, nzt
    983           c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
     1674          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
    9841675                    dt_micro * ddzu(k)
    985           c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
     1676          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
    9861677                    dt_micro * ddzu(k)
    9871678       ENDDO     
    … …  
    9951686             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
    9961687
    997              qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
    998                                                      ABS( d_mean ) )
     1688             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
     1689                                                        2.0_wp * d_max,        &
     1690                                                        ABS( d_mean ) )
    9991691
    10001692             d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) )
    … …  
    10021694             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
    10031695
    1004              nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
    1005                                                      ABS( d_mean ) )
     1696             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
     1697                                                        2.0_wp * d_max,        &
     1698                                                        ABS( d_mean ) )
    10061699          ENDDO
    10071700
    … …  
    10261719          c_run = MIN( 1.0_wp, c_nr(k) )
    10271720          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
    1028              flux  = flux + hyrho(k_run) *                                    &
    1029                      ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *     &
     1721             flux  = flux + hyrho(k_run) *                                     &
     1722                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
    10301723                     0.5_wp ) * c_run * dzu(k_run)
    10311724             z_run = z_run + dzu(k_run)
    … …  
    10361729!--       It is not allowed to sediment more rain drop number density than
    10371730!--       available
    1038           flux = MIN( flux,                                                   &
     1731          flux = MIN( flux,                                                    &
    10391732                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
    10401733
    10411734          sed_nr(k) = flux / dt_micro
    1042           nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /    &
    1043                                  hyrho(k) * dt_micro
     1735          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
     1736                                    hyrho(k) * dt_micro
    10441737!
    10451738!--       Sum up all rain water content which contributes to the flux
    … …  
    10501743          c_run = MIN( 1.0_wp, c_qr(k) )
    10511744
    1052           DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt-1 )
    1053 
    1054              flux  = flux + hyrho(k_run) *                                    &
    1055                      ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *    &
     1745          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
     1746
     1747             flux  = flux + hyrho(k_run) *                                     &
     1748                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
    10561749                     0.5_wp ) * c_run * dzu(k_run)
    10571750             z_run = z_run + dzu(k_run)
    … …  
    10621755!
    10631756!--       It is not allowed to sediment more rain water content than available
    1064           flux = MIN( flux,                                                   &
     1757          flux = MIN( flux,                                                    &
    10651758                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
    10661759
    10671760          sed_qr(k) = flux / dt_micro
    10681761
    1069           qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
     1762          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
    10701763                                hyrho(k) * dt_micro
    1071           q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
     1764          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
    10721765                                hyrho(k) * dt_micro
    1073           pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
     1766          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
    10741767                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
    10751768!
    10761769!--       Compute the rain rate
    1077           prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
    1078                        weight_substep(intermediate_timestep_count)
     1770          IF ( call_microphysics_at_all_substeps )  THEN
     1771             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                  &
     1772                          weight_substep(intermediate_timestep_count)
     1773          ELSE
     1774             prr(k,j,i) = sed_qr(k) / hyrho(k)
     1775          ENDIF
     1776
    10791777       ENDDO
    10801778
    10811779!
    10821780!--    Precipitation amount
    1083        IF ( intermediate_timestep_count == intermediate_timestep_count_max    &
    1084             .AND.  ( dt_do2d_xy - time_do2d_xy ) <                            &
    1085             precipitation_amount_interval )  THEN
    1086 
    1087           precipitation_amount(j,i) = precipitation_amount(j,i) +   &
    1088                                       prr(nzb_s_inner(j,i)+1,j,i) *      &
     1781       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
     1782            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
     1783                   precipitation_amount_interval )  THEN
     1784
     1785          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
     1786                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
    10891787                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
    10901788       ENDIF
    … …  
    10921790    END SUBROUTINE sedimentation_rain_ij
    10931791
    1094 
     1792!------------------------------------------------------------------------------!
     1793! Call for all optimizations
     1794!------------------------------------------------------------------------------!
    10951795!
    10961796!-- This function computes the gamma function (Press et al., 1992).

  • palm/trunk/SOURCE/modules.f90

    r1360 r1361  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! tend_* removed
     23! call_microphysics_at_all_substeps added
     24! default of drizzle set to true
    2325!
    2426! Former revisions:
    … …  
    287289          diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
    288290          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
    289           lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
    290           tend_nr, tend_q, tend_qr, tric, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,   &
    291           v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
     291          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric,       &
     292          u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,       &
     293          w_m_n, w_m_r, w_m_s
    292294
    293295#if defined( __nopointer )
    … …  
    368370                 b_vent = 0.308_wp,    & !: coef. for ventilation effect
    369371                 beta_cc = 3.09E-4_wp, & !: coef. in turb. parametrization (cm-2 s3)
    370                  bfactor,           &
     372                 bfactor,              & !: solution effect on diffusional growth
    371373                 c_1 = 4.82E-6_wp,     & !: coef. in turb. parametrization (m)
    372374                 c_2 = 4.8E-6_wp,      & !: coef. in turb. parametrization (m)
    … …  
    384386                cp = 1005.0_wp,       & !: heat capacity of dry air (J kg-1 K-1)
    385387                diff_coeff_l = 0.23E-4_wp, & !: diffusivity of water vapor (m2 s-1)
    386                 effective_coll_efficiency, & !:
     388                effective_coll_efficiency, & !: effective collision efficiency
    387389                eps_ros = 1.0E-4_wp,  & !: accuracy of Rosenbrock method
    388390                eps_sb = 1.0E-20_wp,  & !: threshold in two-moments scheme
    … …  
    408410                schmidt = 0.71_wp,    & !: Schmidt number
    409411                schmidt_p_1d3,        & !: schmidt**( 1.0 / 3.0 )
     412                sed_qc_const,         & !: const. for sedimentation of cloud water
    410413                sigma_gc = 1.3_wp,    & !: log-normal geometric standard deviation
    411414                stp = 2.5066282746310005_wp, & !: parameter in gamma function
    … …  
    592595                bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., &
    593596                bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE.,&
     597                call_microphysics_at_all_substeps = .FALSE., &
    594598                call_psolver_at_all_substeps = .TRUE., &
    595599                cloud_droplets = .FALSE., cloud_physics = .FALSE., &
    … …  
    604608                do_sum = .FALSE., &
    605609                dp_external = .FALSE., dp_smooth = .FALSE., &
    606                 drizzle = .FALSE., dt_fixed = .FALSE., &
     610                drizzle = .TRUE., dt_fixed = .FALSE., &
    607611                dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
    608612                first_call_lpm = .TRUE., &

  • palm/trunk/SOURCE/parin.f90

    r1360 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! +call_microphysics_at_all_substeps
    2323!
    2424! Former revisions:
    … …  
    148148               bottom_salinityflux, building_height, building_length_x,        &
    149149               building_length_y, building_wall_left, building_wall_south,     &
    150                call_psolver_at_all_substeps, canopy_mode, canyon_height,       &
     150               call_microphysics_at_all_substeps, call_psolver_at_all_substeps,&
     151               canopy_mode, canyon_height,                                     &
    151152               canyon_width_x, canyon_width_y, canyon_wall_left,               &
    152153               canyon_wall_south, cfl_factor,                                  &
    … …  
    248249             bottom_salinityflux, building_height, building_length_x, &
    249250             building_length_y, building_wall_left, building_wall_south, &
    250              call_psolver_at_all_substeps, canopy_mode, canyon_height, &
     251             call_psolver_at_all_substeps, call_microphysics_at_all_substeps, &
     252             canopy_mode, canyon_height, &
    251253             canyon_width_x, canyon_width_y, canyon_wall_left, &
    252254             canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, &

  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r1341 r1361  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
     23! drop concentration (nrsws) added
    2324!
    2425! Former revisions:
    … …  
    6465
    6566    USE arrays_3d,                                                             &
    66         ONLY:  e, pt, q, qs, qsws, rif, shf, ts, u, us, usws, v, vpt, vsws,    &
    67                zu, zw, z0, z0h
     67        ONLY:  e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf,   &
     68               ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
    6869
    6970    USE control_parameters,                                                    &
    70         ONLY:  constant_heatflux, constant_waterflux, coupling_mode, g,        &
    71                humidity, ibc_e_b, kappa, large_scale_forcing, lsf_surf,        &
    72                passive_scalar, pt_surface, q_surface, rif_max, rif_min,        &
    73                run_coupled, surface_pressure
     71        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
     72               coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa,      &
     73               large_scale_forcing, lsf_surf, passive_scalar, precipitation,   &
     74               pt_surface, q_surface, rif_max, rif_min, run_coupled,           &
     75               surface_pressure
    7476
    7577    USE indices,                                                               &
    … …  
    9698!
    9799!-- Data information for accelerators
    98     !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) &
    99     !$acc      present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h )
     100    !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
     101    !$acc      present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v )     &
     102    !$acc      present( vpt, vsws, zu, zw, z0, z0h )
    100103!
    101104!-- Compute theta*
    … …  
    383386          ENDDO
    384387       ENDIF
     388
     389       IF ( cloud_physics  .AND.  icloud_scheme == 0  &
     390            .AND.  precipitation )  THEN
     391
     392          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
     393          !$acc kernels loop independent
     394          DO  i = nxlg, nxrg
     395             !$acc loop independent
     396             DO  j = nysg, nyng
     397
     398                k   = nzb_s_inner(j,i)
     399                z_p = zu(k+1) - zw(k)
     400
     401                IF ( rif(j,i) >= 0.0 )  THEN
     402!
     403!--                Stable stratification
     404                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / (          &
     405                                  LOG( z_p / z0h(j,i) ) +                      &
     406                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
     407                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / (          &
     408                                  LOG( z_p / z0h(j,i) ) +                      &
     409                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
     410
     411                ELSE
     412!
     413!--                Unstable stratification
     414                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
     415                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
     416 
     417                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) /   (        &
     418                             LOG( z_p / z0h(j,i) ) -                           &
     419                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
     420                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) /   (        &
     421                             LOG( z_p / z0h(j,i) ) -                           &
     422                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
     423
     424                ENDIF
     425
     426             ENDDO
     427          ENDDO
     428
     429       ENDIF
     430
    385431    ENDIF
    386432
    … …  
    396442       CALL exchange_horiz_2d( qsws )
    397443       !$acc update device( qsws )
     444       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     445            precipitation )  THEN
     446          !$acc update host( qrsws, nrsws )
     447          CALL exchange_horiz_2d( qrsws )
     448          CALL exchange_horiz_2d( nrsws )
     449          !$acc update device( qrsws, nrsws )
     450       ENDIF
    398451    ENDIF
    399452
    … …  
    425478
    426479!
     480!-- Compute (turbulent) fluxes of rain water content and rain drop concentartion
     481    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
     482       !$OMP PARALLEL DO
     483       !$acc kernels loop independent
     484       DO  i = nxlg, nxrg
     485          !$acc loop independent
     486          DO  j = nysg, nyng
     487             qrsws(j,i) = -qrs(j,i) * us(j,i)
     488             nrsws(j,i) = -nrs(j,i) * us(j,i)
     489          ENDDO
     490       ENDDO
     491    ENDIF
     492
     493!
    427494!-- Bottom boundary condition for the TKE
    428495    IF ( ibc_e_b == 2 )  THEN

  • palm/trunk/SOURCE/prognostic_equations.f90

    r1354 r1361  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Two-moment microphysics moved to the start of prognostic equations. This makes
     23! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
     24! Additionally, it is allowed to call the microphysics just once during the time
     25! step (not at each sub-time step).
     26!
     27! Two-moment cloud physics added for vector and accelerator optimization.
    2328!
    2429! Former revisions:
    … …  
    126131               nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
    127132               q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc, rho,    &
    128                sa, sa_init, sa_p, saswsb, saswst, shf, tend, tend_nr,          &
    129                tend_pt, tend_q, tend_qr, te_m, tnr_m, tpt_m, tq_m, tqr_m,      &
    130                tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, v, vg, vpt, v_p,    &
    131                w, w_p
     133               sa, sa_init, sa_p, saswsb, saswst, shf, tend, te_m, tnr_m,      &
     134               tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, &
     135               v, vg, vpt, v_p, w, w_p
    132136       
    133137    USE control_parameters,                                                    &
    134         ONLY:  cloud_physics, constant_diffusion, cthf, dp_external,           &
     138        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
     139               constant_diffusion, cthf, dp_external,                          &
    135140               dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
    136141               icloud_scheme, inflow_l, intermediate_timestep_count,           &
    … …  
    301306 
    302307       DO  j = nys, nyn
     308!
     309!--       If required, calculate cloud microphysical impacts (two-moment scheme)
     310          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
     311               ( intermediate_timestep_count == 1  .OR.                        &
     312                 call_microphysics_at_all_substeps )                           &
     313             )  THEN
     314             CALL microphysics_control( i, j )
     315          ENDIF
    303316!
    304317!--       Tendency terms for u-velocity component
    … …  
    481494             ENDIF
    482495          ENDIF
    483 !
    484 !--       If required, calculate tendencies for total water content, liquid water
    485 !--       potential temperature, rain water content and rain drop concentration
    486           IF ( cloud_physics  .AND.  icloud_scheme == 0 )  CALL microphysics_control( i, j )
     496
    487497!
    488498!--       If required, compute prognostic equation for potential temperature
    … …  
    511521
    512522!
    513 !--          Using microphysical tendencies (latent heat)
    514              IF ( cloud_physics )  THEN
    515                 IF ( icloud_scheme == 0 )  THEN
    516                    tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i)
    517                 ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
    518                    CALL impact_of_latent_heat( i, j )
    519                 ENDIF
     523!--          If required compute impact of latent heat due to precipitation
     524             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
     525                  precipitation )  THEN
     526                CALL impact_of_latent_heat( i, j )
    520527             ENDIF
    521528
    … …  
    641648     
    642649!
    643 !--          Using microphysical tendencies
    644              IF ( cloud_physics )  THEN
    645                 IF ( icloud_scheme == 0 )  THEN
    646                    tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i)
    647                 ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
    648                    CALL calc_precipitation( i, j )
    649                 ENDIF
     650!--          If required compute decrease of total water content due to
     651!--          precipitation
     652             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
     653                  precipitation )  THEN
     654                CALL calc_precipitation( i, j )             
    650655             ENDIF
    651656!
    … …  
    712717                ENDIF
    713718                CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
    714 !
    715 !--             Using microphysical tendencies (autoconversion, accretion,
    716 !--             evaporation; if required: sedimentation)
    717                 tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
    718719
    719720                CALL user_actions( i, j, 'qr-tendency' )
    … …  
    757758                ENDIF
    758759                CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
    759 !
    760 !--             Using microphysical tendencies (autoconversion, accretion,
    761 !--             selfcollection, breakup, evaporation;
    762 !--             if required: sedimentation)
    763                 tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)
    764760
    765761                CALL user_actions( i, j, 'nr-tendency' )
    … …  
    883879
    884880!
     881!-- If required, calculate cloud microphysical impacts (two-moment scheme)
     882    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
     883         ( intermediate_timestep_count == 1  .OR.                              &
     884           call_microphysics_at_all_substeps )                                 &
     885       )  THEN
     886       CALL cpu_log( log_point(51), 'microphysics', 'start' )
     887       CALL microphysics_control
     888       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
     889    ENDIF
     890
     891!
    885892!-- u-velocity component
    886893    CALL cpu_log( log_point(5), 'u-equation', 'start' )
    … …  
    11561163!
    11571164!--    If required compute impact of latent heat due to precipitation
    1158        IF ( precipitation )  THEN
     1165       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
    11591166          CALL impact_of_latent_heat
    11601167       ENDIF
    … …  
    13481355!--    If required compute decrease of total water content due to
    13491356!--    precipitation
    1350        IF ( precipitation )  THEN
     1357       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
    13511358          CALL calc_precipitation
    13521359       ENDIF
    … …  
    14061413
    14071414       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     1415
     1416!
     1417!--    If required, calculate prognostic equations for rain water content
     1418!--    and rain drop concentration
     1419       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
     1420
     1421          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
     1422
     1423!
     1424!--       Calculate prognostic equation for rain water content
     1425          sbt = tsc(2)
     1426          IF ( scalar_advec == 'bc-scheme' )  THEN
     1427
     1428             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1429!
     1430!--             Bott-Chlond scheme always uses Euler time step. Thus:
     1431                sbt = 1.0_wp
     1432             ENDIF
     1433             tend = 0.0_wp
     1434             CALL advec_s_bc( qr, 'qr' )
     1435
     1436          ENDIF
     1437
     1438!
     1439!--       qr-tendency terms with no communication
     1440          IF ( scalar_advec /= 'bc-scheme' )  THEN
     1441             tend = 0.0_wp
     1442             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1443                IF ( ws_scheme_sca )  THEN
     1444                   CALL advec_s_ws( qr, 'qr' )
     1445                ELSE
     1446                   CALL advec_s_pw( qr )
     1447                ENDIF
     1448             ELSE
     1449                CALL advec_s_up( qr )
     1450             ENDIF
     1451          ENDIF
     1452
     1453          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
     1454
     1455          CALL user_actions( 'qr-tendency' )
     1456
     1457!
     1458!--       Prognostic equation for rain water content
     1459          DO  i = nxl, nxr
     1460             DO  j = nys, nyn
     1461                DO  k = nzb_s_inner(j,i)+1, nzt
     1462                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
     1463                                                     tsc(3) * tqr_m(k,j,i) )   &
     1464                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
     1465                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     1466                ENDDO
     1467             ENDDO
     1468          ENDDO
     1469
     1470!
     1471!--       Calculate tendencies for the next Runge-Kutta step
     1472          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1473             IF ( intermediate_timestep_count == 1 )  THEN
     1474                DO  i = nxl, nxr
     1475                   DO  j = nys, nyn
     1476                      DO  k = nzb_s_inner(j,i)+1, nzt
     1477                         tqr_m(k,j,i) = tend(k,j,i)
     1478                      ENDDO
     1479                   ENDDO
     1480                ENDDO
     1481             ELSEIF ( intermediate_timestep_count < &
     1482                      intermediate_timestep_count_max )  THEN
     1483                DO  i = nxl, nxr
     1484                   DO  j = nys, nyn
     1485                      DO  k = nzb_s_inner(j,i)+1, nzt
     1486                         tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
     1487                                                                   tqr_m(k,j,i)
     1488                      ENDDO
     1489                   ENDDO
     1490                ENDDO
     1491             ENDIF
     1492          ENDIF
     1493
     1494          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
     1495          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
     1496
     1497!
     1498!--       Calculate prognostic equation for rain drop concentration
     1499          sbt = tsc(2)
     1500          IF ( scalar_advec == 'bc-scheme' )  THEN
     1501
     1502             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1503!
     1504!--             Bott-Chlond scheme always uses Euler time step. Thus:
     1505                sbt = 1.0_wp
     1506             ENDIF
     1507             tend = 0.0_wp
     1508             CALL advec_s_bc( nr, 'nr' )
     1509
     1510          ENDIF
     1511
     1512!
     1513!--       nr-tendency terms with no communication
     1514          IF ( scalar_advec /= 'bc-scheme' )  THEN
     1515             tend = 0.0_wp
     1516             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1517                IF ( ws_scheme_sca )  THEN
     1518                   CALL advec_s_ws( nr, 'nr' )
     1519                ELSE
     1520                   CALL advec_s_pw( nr )
     1521                ENDIF
     1522             ELSE
     1523                CALL advec_s_up( nr )
     1524             ENDIF
     1525          ENDIF
     1526
     1527          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
     1528
     1529          CALL user_actions( 'nr-tendency' )
     1530
     1531!
     1532!--       Prognostic equation for rain drop concentration
     1533          DO  i = nxl, nxr
     1534             DO  j = nys, nyn
     1535                DO  k = nzb_s_inner(j,i)+1, nzt
     1536                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
     1537                                                     tsc(3) * tnr_m(k,j,i) )   &
     1538                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
     1539                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     1540                ENDDO
     1541             ENDDO
     1542          ENDDO
     1543
     1544!
     1545!--       Calculate tendencies for the next Runge-Kutta step
     1546          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1547             IF ( intermediate_timestep_count == 1 )  THEN
     1548                DO  i = nxl, nxr
     1549                   DO  j = nys, nyn
     1550                      DO  k = nzb_s_inner(j,i)+1, nzt
     1551                         tnr_m(k,j,i) = tend(k,j,i)
     1552                      ENDDO
     1553                   ENDDO
     1554                ENDDO
     1555             ELSEIF ( intermediate_timestep_count < &
     1556                      intermediate_timestep_count_max )  THEN
     1557                DO  i = nxl, nxr
     1558                   DO  j = nys, nyn
     1559                      DO  k = nzb_s_inner(j,i)+1, nzt
     1560                         tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
     1561                                                                   tnr_m(k,j,i)
     1562                      ENDDO
     1563                   ENDDO
     1564                ENDDO
     1565             ENDIF
     1566          ENDIF
     1567
     1568          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
     1569
     1570       ENDIF
    14081571
    14091572    ENDIF
    … …  
    15101673    ENDIF
    15111674
    1512 
    15131675 END SUBROUTINE prognostic_equations_vector
    15141676
    … …  
    15391701          runge_step = 2
    15401702       ENDIF
     1703    ENDIF
     1704
     1705!
     1706!-- If required, calculate cloud microphysical impacts (two-moment scheme)
     1707    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
     1708         ( intermediate_timestep_count == 1  .OR.                              &
     1709           call_microphysics_at_all_substeps )                                 &
     1710       )  THEN
     1711       CALL cpu_log( log_point(51), 'microphysics', 'start' )
     1712       CALL microphysics_control
     1713       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
    15411714    ENDIF
    15421715
    … …  
    17911964!
    17921965!--    If required compute impact of latent heat due to precipitation
    1793        IF ( precipitation )  THEN
     1966       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
    17941967          CALL impact_of_latent_heat
    17951968       ENDIF
    … …  
    19572130!--    If required compute decrease of total water content due to
    19582131!--    precipitation
    1959        IF ( precipitation )  THEN
     2132       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
    19602133          CALL calc_precipitation
    19612134       ENDIF
    … …  
    19992172
    20002173       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     2174
     2175!
     2176!--    If required, calculate prognostic equations for rain water content
     2177!--    and rain drop concentration
     2178       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
     2179
     2180          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
     2181!
     2182!--       qr-tendency terms with communication
     2183          sbt = tsc(2)
     2184          IF ( scalar_advec == 'bc-scheme' )  THEN
     2185
     2186             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     2187!
     2188!--             Bott-Chlond scheme always uses Euler time step. Thus:
     2189                sbt = 1.0_wp
     2190             ENDIF
     2191             tend = 0.0_wp
     2192             CALL advec_s_bc( qr, 'qr' )
     2193
     2194          ENDIF
     2195
     2196!
     2197!--       qr-tendency terms with no communication
     2198          IF ( scalar_advec /= 'bc-scheme' )  THEN
     2199             tend = 0.0_wp
     2200             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2201                IF ( ws_scheme_sca )  THEN
     2202                   CALL advec_s_ws( qr, 'qr' )
     2203                ELSE
     2204                   CALL advec_s_pw( qr )
     2205                ENDIF
     2206             ELSE
     2207                CALL advec_s_up( qr )
     2208             ENDIF
     2209          ENDIF
     2210
     2211          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
     2212
     2213          CALL user_actions( 'qr-tendency' )
     2214
     2215!
     2216!--       Prognostic equation for rain water content
     2217          DO  i = i_left, i_right
     2218             DO  j = j_south, j_north
     2219                DO  k = nzb_s_inner(j,i)+1, nzt
     2220                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
     2221                                                       tsc(3) * tqr_m(k,j,i) ) &
     2222                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
     2223                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     2224!
     2225!--                Tendencies for the next Runge-Kutta step
     2226                   IF ( runge_step == 1 )  THEN
     2227                      tqr_m(k,j,i) = tend(k,j,i)
     2228                   ELSEIF ( runge_step == 2 )  THEN
     2229                      tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
     2230                                                                tqr_m(k,j,i)
     2231                   ENDIF
     2232                ENDDO
     2233             ENDDO
     2234          ENDDO
     2235
     2236          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
     2237          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
     2238
     2239!
     2240!--       nr-tendency terms with communication
     2241          sbt = tsc(2)
     2242          IF ( scalar_advec == 'bc-scheme' )  THEN
     2243
     2244             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     2245!
     2246!--             Bott-Chlond scheme always uses Euler time step. Thus:
     2247                sbt = 1.0_wp
     2248             ENDIF
     2249             tend = 0.0_wp
     2250             CALL advec_s_bc( nr, 'nr' )
     2251
     2252          ENDIF
     2253
     2254!
     2255!--       nr-tendency terms with no communication
     2256          IF ( scalar_advec /= 'bc-scheme' )  THEN
     2257             tend = 0.0_wp
     2258             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2259                IF ( ws_scheme_sca )  THEN
     2260                   CALL advec_s_ws( nr, 'nr' )
     2261                ELSE
     2262                   CALL advec_s_pw( nr )
     2263                ENDIF
     2264             ELSE
     2265                CALL advec_s_up( nr )
     2266             ENDIF
     2267          ENDIF
     2268
     2269          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
     2270
     2271          CALL user_actions( 'nr-tendency' )
     2272
     2273!
     2274!--       Prognostic equation for rain drop concentration
     2275          DO  i = i_left, i_right
     2276             DO  j = j_south, j_north
     2277                DO  k = nzb_s_inner(j,i)+1, nzt
     2278                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
     2279                                                       tsc(3) * tnr_m(k,j,i) ) &
     2280                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
     2281                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     2282!
     2283!--                Tendencies for the next Runge-Kutta step
     2284                   IF ( runge_step == 1 )  THEN
     2285                      tnr_m(k,j,i) = tend(k,j,i)
     2286                   ELSEIF ( runge_step == 2 )  THEN
     2287                      tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
     2288                                                                tnr_m(k,j,i)
     2289                   ENDIF
     2290                ENDDO
     2291             ENDDO
     2292          ENDDO
     2293
     2294          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
     2295
     2296       ENDIF
    20012297
    20022298    ENDIF
    … …  
    20952391    ENDIF
    20962392
    2097 
    20982393 END SUBROUTINE prognostic_equations_acc
    20992394
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |